home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / gnu / gcc / bug / 2814 < prev    next >
Encoding:
Text File  |  1992-11-21  |  8.1 KB  |  278 lines
  1. Newsgroups: gnu.gcc.bug
  2. Path: sparky!uunet!elroy.jpl.nasa.gov!sdd.hp.com!zaphod.mps.ohio-state.edu!cis.ohio-state.edu!ratface.chem.pitt.edu!walsh
  3. From: walsh@ratface.chem.pitt.edu (alan walsh)
  4. Subject: "bug" in gcc 2.2.2
  5. Message-ID: <9211211953.AA05885@life.ai.mit.edu>
  6. Sender: gnulists@ai.mit.edu
  7. Organization: GNUs Not Usenet
  8. Distribution: gnu
  9. Date: Sat, 21 Nov 1992 09:48:28 GMT
  10. Approved: bug-gcc@prep.ai.mit.edu
  11. Lines: 265
  12.  
  13. I have come across the following problem with gcc.  when i compile the
  14. file included here using -O or -O2 on the command line (gcc -O temp.c)
  15. the following error occurs:
  16.  
  17. ======================================================================
  18. temp.c: In function `main':
  19. temp.c:217: internal error--insn does not satisfy its constraints:
  20. (insn 4667 4665 2071 (set (reg:SI 64)
  21.        (mem:SI (reg:SI 13))) 34 {movsi+3} (nil)
  22.    (nil))
  23. gcc: Internal compiler error: program cc1 got fatal signal 6
  24. ======================================================================
  25.  
  26. the problem disappears when the 'for' statement on line 60 of the
  27. following file is uncommented.  admittedly the code is poorly commented.
  28. if there are any questions, please let me know.  i'll be glad to provide
  29. and addn'l needed information.
  30.  
  31. ======================================================================
  32. #include <stdio.h>
  33. #include <math.h>
  34.  
  35. #define L    32
  36. #define LM    (L-1)
  37. #define L3    (L*L*L)
  38. #define FTOL    1.0e-4
  39.  
  40. main()
  41. {
  42.   extern double fmin3();
  43.  
  44.   int a,b,c,d,r12;
  45.   int nxm,nxp,nym,nyp,nzm,nzp,nxyz;
  46.   int np, npl, nmin, it, itave;        /* number of poly balls */
  47.   int density;                /* counter-ion density */
  48.   int nimp;                /* number of counter ions */
  49.   int locpol[3][20];             /* location of poly balls */
  50.   int pbp[L],pbm[L];            /* indices for PB cond's */
  51.   double amf,action,feqn,nit;
  52.   double sumexp, sumexm, gammap, gammam;
  53.   double cp, cm, ep, ef, phiold;
  54.   double t1, t2, dpavg, dptot;
  55.   double phi[L3],vexc[L3];         /* field, poly-repulsions */
  56.   double beta, z;
  57. /* double kludge[75][10]; */
  58.  
  59.  
  60. /*
  61.       size = 10;
  62.       cinp(kludge,size);
  63. */
  64.   fscanf(stdin,"%d\n",&np);
  65.   for (a=0;a<np;a++)
  66.     fscanf(stdin,"%d %d %d\n",&locpol[0][a],&locpol[1][a],&locpol[2][a]);
  67.   np = 3;
  68.  
  69.   for (a=0;a<L;a++) {
  70.     pbp[a] = (a+1)%L;
  71.     pbm[a] = (a+LM)%L;
  72.   }
  73.   fscanf(stdin,"%le %le\n",&beta,&z);
  74.   fscanf(stdin,"%d\n",&density);
  75.   nimp=density*L3;
  76.  
  77.   nit = 0.0; itave = 0;
  78.   for (a=0;a<L3;a++) vexc[a] = 0.0;
  79.   locpol[0][0] = locpol[1][0] = locpol[2][0] = 0;
  80.   locpol[0][1] = locpol[1][1] = locpol[2][1] = 0;
  81.   locpol[0][2] = locpol[1][2] = locpol[2][2] = 0;
  82.  
  83. /* the first two colloids are position along the x axis, and the third is
  84.    placed on all the other lattice points in the 1/4th plane nearest the
  85.    orign.  the configurations with both "extra" colloids along the x axis
  86.    are not double counted due to the test that locpol[0][2] > locpol[0][1].
  87.    only one "plane" is needed since in principal the sysem has a c_inf
  88.    axis along any pair of colloids (3 points determine a plane...)
  89.    These configurations represent all the unique 3-body configurations. */
  90.   for ( locpol[0][1] = 1; locpol[0][1] < L/2; locpol[0][1]++ )
  91.   /* for ( locpol[1][2] = 0; locpol[1][2] < L/2; locpol[1][2]++ ) */
  92.   for ( locpol[0][2] = 0; locpol[0][2] < L/2; locpol[0][2]++ )
  93.   if ( locpol[1][2] || locpol[2][2] || (locpol[0][2] > locpol[0][1]) ) {
  94.     for (a=0;a<np;a++)
  95.       vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 250.0;
  96.     for (a=0;a<L3;a++)
  97.       phi[a]=0.0;
  98. /* set (npl,nmin)=(total # + simple ions,total # - simple ions) */
  99.     npl=np*z + nimp;
  100.     nmin=nimp;
  101. /* intialize gamma, based on DH level phi's: */
  102.     sumexp=0.0;
  103.     sumexm=0.0;
  104.     for (a=0;a<L3;a++) {
  105.       sumexm += exp(-phi[a]-vexc[a]);
  106.       sumexp += exp(phi[a]-vexc[a]);
  107.     }
  108.     gammap = npl/sumexp;
  109.     gammam = nmin/sumexm;
  110.     /* printf("initial gammap= %e %e %d\n",gammap,sumexp,L3); */
  111.     /* printf("initial gammam= %e %e\n",gammam,sumexm); */
  112.  
  113.     action=0.0;
  114.     for(a=0;a<np;a++) 
  115.       action -= phi[locpol[0][a],locpol[1][a],locpol[2][a]];
  116.     action *= z;
  117.  
  118.     for (a=0,nxyz=0;a<L3;a+=L*L) {
  119.       nxp = pbp[a/(L*L)]*L*L;
  120.       nxm = pbm[a/(L*L)]*L*L;
  121.       for (b=0;b<L*L;b+=L) {
  122.         nyp = pbp[b/L]*L;
  123.         nym = pbm[b/L]*L;
  124.         for (c=0;c<L;c++,nxyz++) {
  125.           nzp = pbp[c];
  126.           nzm = pbm[c];
  127.       action += -0.5*beta*phi[nxyz]*
  128.            (phi[nxp+b+c] + phi[a+nyp+c] + phi[a+b+nzp] +
  129.             phi[nxm+b+c] + phi[a+nym+c] + phi[a+b+nzm] - 6*phi[nxyz]) +
  130.             gammap*exp(phi[nxyz]-vexc[nxyz]) +
  131.             gammam*exp(-phi[nxyz]-vexc[nxyz]);
  132.     }
  133.       }
  134.     }
  135.     /* printf("initial action = %le\n",action); */
  136.  
  137.     feqn = 2*FTOL;
  138.     nit += 1.0;
  139.     it = 0;
  140.     while ( fabs(feqn) > FTOL){
  141.       it++; itave++;
  142.       dptot=0.0;
  143.       dpavg=2.0;
  144.  
  145.       for (a=0,nxyz=0;a<L3;a+=L*L) {
  146.         nxp = pbp[a/(L*L)]*L*L;
  147.         nxm = pbm[a/(L*L)]*L*L;
  148.         for (b=0;b<L*L;b+=L) {
  149.           nyp = pbp[b/L]*L;
  150.           nym = pbm[b/L]*L;
  151.           for (c=0;c<L;c++,nxyz++) {
  152.             nzp = pbp[c];
  153.             nzm = pbm[c];
  154.             t1 = -beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
  155.                  + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm]);
  156.             for (d=0;d<np;d++)
  157.               if ( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
  158.                 t1 -= z;
  159.         if ( vexc[nxyz] < 1.0 ) {
  160.               cp = gammap;
  161.               cm = gammam;
  162.         }
  163.         else { cp = cm = 0.0; }
  164.             phiold = phi[nxyz];
  165.             phi[nxyz] = fmin3(3*beta,t1,cp,cm,phiold);
  166.             dptot = fabs(phiold-phi[nxyz])+dptot;
  167.           }
  168.     }
  169.       }
  170.       dpavg = dptot/L3;
  171.  
  172.       action=0.0;
  173.       sumexp=0.0;
  174.       sumexm=0.0;
  175.       feqn=0.0;
  176.  
  177.       for (a=0,nxyz=0;a<L3;a+=L*L) {
  178.         nxp = pbp[a/(L*L)]*L*L;
  179.         nxm = pbm[a/(L*L)]*L*L;
  180.         for (b=0;b<L*L;b+=L) {
  181.           nyp = pbp[b/L]*L;
  182.           nym = pbm[b/L]*L;
  183.           for (c=0;c<L;c++,nxyz++) {
  184.             nzp = pbp[c];
  185.             nzm = pbm[c];
  186.         t1 = exp(phi[nxyz]-vexc[nxyz]);
  187.         t2 = exp(-phi[nxyz]-vexc[nxyz]);
  188.             ep = beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
  189.                  + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm]-6*phi[nxyz]);
  190.             action -= 0.5*phi[nxyz]*ep;
  191.             ef = -ep + gammap*t1 - gammam*t2;
  192.             for (d=0;d<np;d++)
  193.               if( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
  194.                 ef -= z;
  195.             feqn += ef*ef;
  196.             sumexp += t1;
  197.             sumexm += t2;
  198.           }
  199.     }
  200.       }
  201.       feqn=feqn/L3;
  202.       feqn=sqrt(feqn);
  203.       action += gammap*sumexp + gammam*sumexm;
  204.       for (a=0;a<np;a++)
  205.         action -= z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
  206.       if ( ! (it % 10) ) {
  207.         gammap=npl/sumexp;
  208.         gammam=nmin/sumexm;
  209.     /*
  210.     printf("action=%le accuracy=%le\n gammap=%le gammam=%le\n",
  211.          action,feqn,gammap,gammam);
  212.     */
  213.       }
  214.     }
  215.  
  216. /*
  217. C  compute Helmholtz free energy (meanfield)
  218. */
  219.     amf=0.0;
  220.     for (a=0;a<L*L*L;a++)
  221.       amf += gammap*(0.5*phi[a]-1.0)*exp(phi[a]-vexc[a])
  222.         -gammam*(0.5*phi[a]+1.0)*exp(-phi[a]-vexc[a]);
  223.     for (a=0;a<np;a++)
  224.       amf += 0.5*z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
  225.     amf += npl*log(gammap)+nmin*log(gammam);
  226.     printf("%d %d %d %d %20.16le\n",locpol[0][1],locpol[0][2],
  227.                    locpol[1][2],locpol[2][2],amf);
  228.     for (a=0;a<np;a++)
  229.       vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 0.0;
  230.   }
  231.       fprintf(stderr,"%le\n",itave/nit);
  232.       /*
  233.       write(numout,80) (((phi(nx,ny,nz),nz=0,lm),ny=0,lm),nx=0,lm)
  234. 80    format(4d16.9)    
  235.       do 120 nx=0,lm
  236.       write(17,21) 1.*nx,phi(nx,0,0)
  237.       if(nx.gt.0) then
  238.       write(18,21) 1.*nx,gammap*dexp(phi(nx,0,0))
  239.       else
  240.       end if
  241. 21    format(13f10.5)
  242. 120   continue
  243.       write(17,21) (lm+1)*1.,phi(0,0,0)
  244.       stop
  245.       end
  246.       */
  247.  
  248. }
  249.  
  250.  
  251. /*
  252. c uses newton-raphson scheme to find root of
  253. c f(x)=2*a*x + b + cp*exp(x) - cm*exp(-x)
  254. */
  255. double fmin3(a,b,cp,cm,phiold)
  256. double a,b,cp,cm,phiold;
  257. {
  258.   int ic;
  259.   double x0, f0, fpr0, xrt, frt;
  260.  
  261.   x0=phiold;
  262.   f0=2*a*x0 + b + cp*exp(x0) - cm*exp(-x0);
  263.   for (ic=0;ic<20;ic++) {
  264.     fpr0= 2*a + cp*exp(x0) + cm*exp(-x0);
  265.   /* approximate root: */
  266.     xrt=x0 - f0/fpr0;
  267.     frt=2*a*xrt + b + cp*exp(xrt) - cm*exp(-xrt);
  268.     if( fabs(frt) < 0.01 ) return(xrt);
  269.  
  270.     x0=xrt;
  271.     f0=frt;
  272.   }
  273.   fprintf(stderr,"root not found in fmin3\n");
  274.   exit(0);
  275.   return(1.0);
  276. }
  277.  
  278.