home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #27 / NN_1992_27.iso / spool / gnu / gcc / bug / 2816 < prev    next >
Encoding:
Text File  |  1992-11-21  |  8.0 KB  |  257 lines
  1. Newsgroups: gnu.gcc.bug
  2. Path: sparky!uunet!cs.utexas.edu!sdd.hp.com!saimiri.primate.wisc.edu!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: similar problem, different code...
  5. Message-ID: <9211212043.AA07861@life.ai.mit.edu>
  6. Sender: gnulists@ai.mit.edu
  7. Organization: GNUs Not Usenet
  8. Distribution: gnu
  9. Date: Sat, 21 Nov 1992 10:38:34 GMT
  10. Approved: bug-gcc@prep.ai.mit.edu
  11. Lines: 244
  12.  
  13. in the same vein as my earlier bug report, i ran into the following problem.
  14. when line 222 in the following code is uncommented i get the following
  15. compiler crash
  16.  
  17. ================================================================================
  18. fft.c: In function `main':
  19. fft.c:226: internal error--insn does not satisfy its constraints:
  20. (insn 4175 741 716 (set (reg:SI 68)
  21.        (mem:SI (plus:SI (reg:SI 30)
  22.                (const_int -380)))) 34 {movsi+3} (nil)
  23.    (nil))
  24. gcc: Internal compiler error: program cc1 got fatal signal 6
  25. ================================================================================
  26.  
  27. the offending (apparently) line is the call to exit(0);
  28.  
  29. ================================================================================
  30. #include <stdio.h>
  31. #include <math.h>
  32.  
  33. #define L    32
  34. #define LM    (L-1)
  35. #define L3    (L*L*L)
  36. #define FTOL    2.0e-5
  37.  
  38. main()
  39. {
  40.   extern void fourn();
  41.  
  42.   int a,b,c,d,e,f,r12;
  43.   int nxm,nxp,nym,nyp,nzm,nzp,nxyz;
  44.   int np, npl, nmin, it, itave;        /* number of poly balls */
  45.   int density;                /* counter-ion density */
  46.   int nimp;                /* number of counter ions */
  47.   int locpol[3][20];             /* location of poly balls */
  48.   int pbp[L],pbm[L];            /* indices for PB cond's */
  49.   int nn[3], dim, sign;
  50.   double amf,action,feqn,nit;
  51.   double sumexp, sumexm, gammap, gammam;
  52.   double cp, cm, ep, ef, phiold;
  53.   double t1, t2;
  54.   double phi[L3],vexc[L3];         /* field, poly-repulsions */
  55.   double beta, z;
  56.   double pil,il3;
  57.   float dphi[L3][2],psq[L3];
  58.   double eps,psqcut;
  59.  
  60.   pil = M_PI/L;
  61.   il3 = 1.0/L3;
  62.   fscanf(stdin,"%d\n",&np);
  63.   for (a=0;a<np;a++)
  64.     fscanf(stdin,"%d %d %d\n",&locpol[0][a],&locpol[1][a],&locpol[2][a]);
  65.   np = 3;
  66.  
  67.   for (a=0;a<L;a++) {
  68.     pbp[a] = (a+1)%L;
  69.     pbm[a] = (a+LM)%L;
  70.   }
  71.   fscanf(stdin,"%le %le\n",&beta,&z);
  72.   fscanf(stdin,"%d\n",&density);
  73.   nimp=density*L3;
  74.   /* nimp = 4000; */
  75.  
  76.   fscanf(stdin,"%le %le\n",&eps,&psqcut);
  77. /* initialize needed stuff for fft acceleration */
  78.   nn[0] = L;
  79.   nn[1] = L;
  80.   nn[2] = L;
  81.   dim = 3;
  82.   for (a=0;a<L;a++)
  83.     for (b=0;b<L;b++)
  84.       for (c=0;c<L;c++)
  85.     if ( a+b+c == 0 ) psq[0] = 1.0;
  86.     else psq[a+b*L+c*L*L]= 3.0/
  87.       (pow(sin(a*pil),2.0)+pow(sin(b*pil),2.0)+pow(sin(c*pil),2.0)+psqcut);
  88.  
  89.   nit = 0.0; itave = 0;
  90.   for (a=0;a<L3;a++) vexc[a] = 0.0;
  91.   locpol[0][0] = locpol[1][0] = locpol[2][0] = 0;
  92.   locpol[0][1] = locpol[1][1] = locpol[2][1] = 0;
  93.   locpol[0][2] = locpol[1][2] = locpol[2][2] = 0;
  94.  
  95. /* the first two colloids are position along the x axis, and the third is
  96.    placed on all the other lattice points in the 1/4th plane nearest the
  97.    orign.  the configurations with both "extra" colloids along the x axis
  98.    are not double counted due to the test that locpol[0][2] > locpol[0][1].
  99.    only one "plane" is needed since in principal the sysem has a c_inf
  100.    axis along any pair of colloids (3 points determine a plane...)
  101.    These configurations represent all the unique 3-body configurations. */
  102.   for ( locpol[0][1] = 3; locpol[0][1] < L/2; locpol[0][1]++ )
  103.   for ( locpol[1][2] = 0; locpol[1][2] < L/2; locpol[1][2]++ )
  104.   for ( locpol[0][2] = 0; locpol[0][2] < L/2; locpol[0][2]++ )
  105.   if ( locpol[1][2] || locpol[2][2] || (locpol[0][2] > locpol[0][1]) ) {
  106.     for (a=0;a<np;a++)
  107.       vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 250.0;
  108.     for (a=0;a<L3;a++)
  109.       phi[a]=0.0;
  110. /* set (npl,nmin)=(total # + simple ions,total # - simple ions) */
  111.     npl=np*z + nimp;
  112.     nmin=nimp;
  113. /* intialize gamma, based on DH level phi's: */
  114.     sumexp=0.0;
  115.     sumexm=0.0;
  116.     for (a=0;a<L3;a++) {
  117.       sumexm += exp(-phi[a]-vexc[a]);
  118.       sumexp += exp(phi[a]-vexc[a]);
  119.     }
  120.     gammap = npl/sumexp;
  121.     gammam = nmin/sumexm;
  122.     printf("initial gammap= %e %e %d\n",gammap,sumexp,L3);
  123.     printf("initial gammam= %e %e\n",gammam,sumexm);
  124.  
  125.     action=0.0;
  126.     for(a=0;a<np;a++) 
  127.       action -= phi[locpol[0][a]*L*L + locpol[1][a]*L + locpol[2][a]];
  128.     action *= z;
  129.  
  130.     for (a=0,nxyz=0;a<L3;a+=L*L) {
  131.       nxp = pbp[a/(L*L)]*L*L;
  132.       nxm = pbm[a/(L*L)]*L*L;
  133.       for (b=0;b<L*L;b+=L) {
  134.         nyp = pbp[b/L]*L;
  135.         nym = pbm[b/L]*L;
  136.         for (c=0;c<L;c++,nxyz++) {
  137.           nzp = pbp[c];
  138.           nzm = pbm[c];
  139.       action += -0.5*beta*phi[nxyz]*
  140.            (phi[nxp+b+c] + phi[a+nyp+c] + phi[a+b+nzp] +
  141.             phi[nxm+b+c] + phi[a+nym+c] + phi[a+b+nzm] - 6*phi[nxyz]) +
  142.             gammap*exp(phi[nxyz]-vexc[nxyz]) +
  143.             gammam*exp(-phi[nxyz]-vexc[nxyz]);
  144.     }
  145.       }
  146.     }
  147.     printf("initial action = %le\n",action);
  148.  
  149.     feqn = 2*FTOL;
  150.     nit += 1.0;
  151.     it = 0;
  152.     while ( fabs(feqn) > FTOL ){
  153.       it++; itave++;
  154.  
  155.       for (a=f=0,nxyz=0;a<L3;a+=L*L,f++) {
  156.         nxp = pbp[a/(L*L)]*L*L;
  157.         nxm = pbm[a/(L*L)]*L*L;
  158.         for (b=0;b<L*L;b+=L) {
  159.           nyp = pbp[b/L]*L;
  160.           nym = pbm[b/L]*L;
  161.           for (c=e=0;c<L;c++,nxyz++,e+=L*L) {
  162.             nzp = pbp[c];
  163.             nzm = pbm[c];
  164.             t1 = -beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
  165.                  + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm] - 6*phi[nxyz]);
  166.             for (d=0;d<np;d++)
  167.               if ( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
  168.                 t1 -= z;
  169.             t1 += gammap*exp(phi[nxyz]-vexc[nxyz]) -
  170.                gammam*exp(-phi[nxyz]-vexc[nxyz]);
  171.             dphi[f+b+e][0] = eps*t1;
  172.             dphi[f+b+e][1] = 0;
  173.           }
  174.     }
  175.       }
  176.       sign = 1;
  177.       fourn(dphi,nn,&dim,&sign);
  178. /* fourier accelerate */
  179.       for (a=0;a<L3;a++) {
  180.         dphi[a][0] *= psq[a];
  181.         dphi[a][1] *= psq[a];
  182.       }
  183.       sign = -1;
  184.       fourn(dphi,nn,&dim,&sign);
  185.  
  186.       for (nxyz=f=0;f<L;f++)
  187.         for (b=0;b<L*L;b+=L)
  188.           for (d=0;d<L3;d+=L*L)
  189.             phi[nxyz]=phi[nxyz]-il3*dphi[f+b+d][0];
  190.       /* same for these four lines
  191.       for (f=0;f<L;f++)
  192.         for (b=0;b<L;b++)
  193.           for (d=0;d<L;d++)
  194.             phi[f*L*L+b*L+d]=phi[f*L*L+b*L+d]-il3*dphi[f+b*L+d*L*L][0];
  195.       */
  196.  
  197.       action=0.0;
  198.       sumexp=0.0;
  199.       sumexm=0.0;
  200.       feqn=0.0;
  201.  
  202.       for (a=0,nxyz=0;a<L3;a+=L*L) {
  203.         nxp = pbp[a/(L*L)]*L*L;
  204.         nxm = pbm[a/(L*L)]*L*L;
  205.         for (b=0;b<L*L;b+=L) {
  206.           nyp = pbp[b/L]*L;
  207.           nym = pbm[b/L]*L;
  208.           for (c=0;c<L;c++,nxyz++) {
  209.             nzp = pbp[c];
  210.             nzm = pbm[c];
  211.         t1 = exp(phi[nxyz]-vexc[nxyz]);
  212.         t2 = exp(-phi[nxyz]-vexc[nxyz]);
  213.             ep = beta*(phi[nxp+b+c]+phi[a+nyp+c]+phi[a+b+nzp]
  214.                  + phi[nxm+b+c]+phi[a+nym+c]+phi[a+b+nzm]-6*phi[nxyz]);
  215.             action -= 0.5*phi[nxyz]*ep;
  216.             ef = -ep + gammap*t1 - gammam*t2;
  217.             for (d=0;d<np;d++)
  218.               if( nxyz == locpol[0][d]*L*L+locpol[1][d]*L+locpol[2][d] )
  219.                 ef -= z;
  220.             feqn += ef*ef;
  221.             sumexp += t1;
  222.             sumexm += t2;
  223.           }
  224.     }
  225.       }
  226.       feqn=feqn/L3;
  227.       feqn=sqrt(feqn);
  228.       action += gammap*sumexp + gammam*sumexm;
  229.       for (a=0;a<np;a++)
  230.         action -= z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
  231.       printf("action=%le accuracy=%le\n",action,feqn);
  232.       gammap=npl/sumexp;
  233.       gammam=nmin/sumexm;
  234.       printf("gammap=%le gammam=%le\n",gammap,gammam);
  235.     }
  236.  
  237. /*
  238. C  compute Helmholtz free energy (meanfield)
  239. */
  240.     amf=0.0;
  241.     for (a=0;a<L*L*L;a++)
  242.       amf += gammap*(0.5*phi[a]-1.0)*exp(phi[a]-vexc[a])
  243.         -gammam*(0.5*phi[a]+1.0)*exp(-phi[a]-vexc[a]);
  244.     for (a=0;a<np;a++)
  245.       amf += 0.5*z*phi[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]];
  246.     amf += npl*log(gammap)+nmin*log(gammam);
  247.     printf("%d %d %d %d %20.16le\n",locpol[0][1],locpol[0][2],
  248.                    locpol[1][2],locpol[2][2],amf);
  249.     for (a=0;a<np;a++)
  250.       vexc[locpol[0][a]*L*L+locpol[1][a]*L+locpol[2][a]] = 0.0;
  251.     /* commenting out the next line apparently solves the compile problem... */
  252.     exit(0);
  253.   }
  254.       fprintf(stderr,"%le\n",itave/nit);
  255. }
  256.  
  257.