🚀 go-pugleaf

I have come across the following problem with gcc.  when i compile the
file included here using -O or -O2 on the command line (gcc -O temp.c)
the following error occurs:

======================================================================
temp.c: In function `main':
temp.c:217: internal error--insn does not satisfy its constraints:
(insn 4667 4665 2071 (set (reg:SI 64)
       (mem:SI (reg:SI 13))) 34 {movsi+3} (nil)
   (nil))
gcc: Internal compiler error: program cc1 got fatal signal 6
======================================================================

the problem disappears when the 'for' statement on line 60 of the
following file is uncommented.  admittedly the code is poorly commented.
if there are any questions, please let me know.  i'll be glad to provide
and addn'l needed information.

======================================================================
#include <stdio.h>
#include <math.h>

#define L	32
#define LM	(L-1)
#define L3	(L*L*L)
#define FTOL	1.0e-4

main()
{
  extern double fmin3();

  int a,b,c,d,r12;
  int nxm,nxp,nym,nyp,nzm,nzp,nxyz;
  int np, npl, nmin, it, itave;		/* number of poly balls */
  int density;				/* counter-ion density */
  int nimp;				/* number of counter ions */
  int locpol[3][20]; 			/* location of poly balls */
  int pbp[L],pbm[L];			/* indices for PB cond's */
  double amf,action,feqn,nit;
  double sumexp, sumexm, gammap, gammam;
  double cp, cm, ep, ef, phiold;
  double t1, t2, dpavg, dptot;
  double phi[L3],vexc[L3];	 	/* field, poly-repulsions */
  double beta, z;
/* double kludge[75][10]; */


/*
      size = 10;
      cinp(kludge,size);
*/
  fscanf(stdin,"%d\n",&np);
  for (a=0;a<np;a++)
    fscanf(stdin,"%d %d %d\n",&locpol[0][a],&locpol[1][a],&locpol[2][a]);
  np = 3;

  for (a=0;a<L;a++) {
    pbp[a] = (a+1)%L;
    pbm[a] = (a+LM)%L;
  }
  fscanf(stdin,"%le %le\n",&beta,&z);
  fscanf(stdin,"%d\n",&density);
  nimp=density*L3;

  nit = 0.0; itave = 0;
  for (a=0;a<L3;a++) vexc[a] = 0.0;
  locpol[0][0] = locpol[1][0] = locpol[2][0] = 0;
  locpol[0][1] = locpol[1][1] = locpol[2][1] = 0;
  locpol[0][2] = locpol[1][2] = locpol[2][2] = 0;

/* the first two colloids are position along the x axis, and the third is
   placed on all the other lattice points in the 1/4th plane nearest the
   orign.  the configurations with both "extra" colloids along the x axis
   are not double counted due to the test that locpol[0][2] > locpol[0][1].
   only one "plane" is needed since in principal the sysem has a c_inf
   axis along any pair of colloids (3 points determine a plane...)
   These configurations represent all the unique 3-body configurations. */
  for ( locpol[0][1] = 1; locpol[0][1] < L/2; locpol[0][1]++ )
  /* for ( locpol[1][2] = 0; locpol[1][2] < L/2; locpol[1][2]++ ) */
  for ( locpol[0][2] = 0; locpol[0][2] < L/2; locpol[0][2]++ )
  if ( locpol[1][2] || locpol[2][2] || (locpol[0][2] > locpol[0][1]) ) {
    for (a=0;a<np;a++)
      vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 250.0;
    for (a=0;a<L3;a++)
      phi[a]=0.0;
/* set (npl,nmin)=(total # + simple ions,total # - simple ions) */
    npl=np*z + nimp;
    nmin=nimp;
/* intialize gamma, based on DH level phi's: */
    sumexp=0.0;
    sumexm=0.0;
    for (a=0;a<L3;a++) {
      sumexm += exp(-phi[a]-vexc[a]);
      sumexp += exp(phi[a]-vexc[a]);
    }
    gammap = npl/sumexp;
    gammam = nmin/sumexm;
    /* printf("initial gammap= %e %e %d\n",gammap,sumexp,L3); */
    /* printf("initial gammam= %e %e\n",gammam,sumexm); */

    action=0.0;
    for(a=0;a<np;a++) 
      action -= phi[locpol[0][a],locpol[1][a],locpol[2][a]];
    action *= z;

    for (a=0,nxyz=0;a<L3;a+=L*L) {
      nxp = pbp[a/(L*L)]*L*L;
      nxm = pbm[a/(L*L)]*L*L;
      for (b=0;b<L*L;b+=L) {
        nyp = pbp[b/L]*L;
        nym = pbm[b/L]*L;
        for (c=0;c<L;c++,nxyz++) {
          nzp = pbp[c];
          nzm = pbm[c];
	  action += -0.5*beta*phi[nxyz]*
		   (phi[nxp+b+c] + phi[a+nyp+c] + phi[a+b+nzp] +
		    phi[nxm+b+c] + phi[a+nym+c] + phi[a+b+nzm] - 6*phi[nxyz]) +
		    gammap*exp(phi[nxyz]-vexc[nxyz]) +
		    gammam*exp(-phi[nxyz]-vexc[nxyz]);
	}
      }
    }
    /* printf("initial action = %le\n",action); */

    feqn = 2*FTOL;
    nit += 1.0;
    it = 0;
    while ( fabs(feqn) > FTOL){
      it++; itave++;
      dptot=0.0;
      dpavg=2.0;

      for (a=0,nxyz=0;a<L3;a+=L*L) {
        nxp = pbp[a/(L*L)]*L*L;
        nxm = pbm[a/(L*L)]*L*L;
        for (b=0;b<L*L;b+=L) {
          nyp = pbp[b/L]*L;
          nym = pbm[b/L]*L;
          for (c=0;c<L;c++,nxyz++) {
            nzp = pbp[c];
            nzm = pbm[c];
            t1 = -beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
                 + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm]);
            for (d=0;d<np;d++)
              if ( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
                t1 -= z;
	    if ( vexc[nxyz] < 1.0 ) {
              cp = gammap;
              cm = gammam;
	    }
	    else { cp = cm = 0.0; }
            phiold = phi[nxyz];
            phi[nxyz] = fmin3(3*beta,t1,cp,cm,phiold);
            dptot = fabs(phiold-phi[nxyz])+dptot;
          }
	}
      }
      dpavg = dptot/L3;

      action=0.0;
      sumexp=0.0;
      sumexm=0.0;
      feqn=0.0;

      for (a=0,nxyz=0;a<L3;a+=L*L) {
        nxp = pbp[a/(L*L)]*L*L;
        nxm = pbm[a/(L*L)]*L*L;
        for (b=0;b<L*L;b+=L) {
          nyp = pbp[b/L]*L;
          nym = pbm[b/L]*L;
          for (c=0;c<L;c++,nxyz++) {
            nzp = pbp[c];
            nzm = pbm[c];
	    t1 = exp(phi[nxyz]-vexc[nxyz]);
	    t2 = exp(-phi[nxyz]-vexc[nxyz]);
            ep = beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
                 + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm]-6*phi[nxyz]);
            action -= 0.5*phi[nxyz]*ep;
            ef = -ep + gammap*t1 - gammam*t2;
            for (d=0;d<np;d++)
              if( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
                ef -= z;
            feqn += ef*ef;
            sumexp += t1;
            sumexm += t2;
          }
	}
      }
      feqn=feqn/L3;
      feqn=sqrt(feqn);
      action += gammap*sumexp + gammam*sumexm;
      for (a=0;a<np;a++)
        action -= z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
      if ( ! (it % 10) ) {
        gammap=npl/sumexp;
        gammam=nmin/sumexm;
	/*
	printf("action=%le accuracy=%le\n gammap=%le gammam=%le\n",
		 action,feqn,gammap,gammam);
	*/
      }
    }

/*
C  compute Helmholtz free energy (meanfield)
*/
    amf=0.0;
    for (a=0;a<L*L*L;a++)
      amf += gammap*(0.5*phi[a]-1.0)*exp(phi[a]-vexc[a])
		-gammam*(0.5*phi[a]+1.0)*exp(-phi[a]-vexc[a]);
    for (a=0;a<np;a++)
      amf += 0.5*z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
    amf += npl*log(gammap)+nmin*log(gammam);
    printf("%d %d %d %d %20.16le\n",locpol[0][1],locpol[0][2],
			       locpol[1][2],locpol[2][2],amf);
    for (a=0;a<np;a++)
      vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 0.0;
  }
      fprintf(stderr,"%le\n",itave/nit);
      /*
      write(numout,80) (((phi(nx,ny,nz),nz=0,lm),ny=0,lm),nx=0,lm)
80    format(4d16.9)    
      do 120 nx=0,lm
      write(17,21) 1.*nx,phi(nx,0,0)
      if(nx.gt.0) then
      write(18,21) 1.*nx,gammap*dexp(phi(nx,0,0))
      else
      end if
21    format(13f10.5)
120   continue
      write(17,21) (lm+1)*1.,phi(0,0,0)
      stop
      end
      */

}


/*
c uses newton-raphson scheme to find root of
c f(x)=2*a*x + b + cp*exp(x) - cm*exp(-x)
*/
double fmin3(a,b,cp,cm,phiold)
double a,b,cp,cm,phiold;
{
  int ic;
  double x0, f0, fpr0, xrt, frt;

  x0=phiold;
  f0=2*a*x0 + b + cp*exp(x0) - cm*exp(-x0);
  for (ic=0;ic<20;ic++) {
    fpr0= 2*a + cp*exp(x0) + cm*exp(-x0);
  /* approximate root: */
    xrt=x0 - f0/fpr0;
    frt=2*a*xrt + b + cp*exp(xrt) - cm*exp(-xrt);
    if( fabs(frt) < 0.01 ) return(xrt);

    x0=xrt;
    f0=frt;
  }
  fprintf(stderr,"root not found in fmin3\n");
  exit(0);
  return(1.0);
}