home *** CD-ROM | disk | FTP | other *** search
/ Usenet 1994 October / usenetsourcesnewsgroupsinfomagicoctober1994disk2.iso / misc / volume21 / cloops / part02 / kernel.c < prev   
C/C++ Source or Header  |  1991-07-25  |  12KB  |  498 lines
  1. /*
  2.  * This file is part of the Livermore Loops transliteration into C.
  3.  * Copyright (C) 1991 by Martin Fouts
  4.  *
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 1, or (at your option)
  8.  * any later version.
  9.  *
  10.  * This program is distributed in the hope that it will be useful,
  11.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13.  * GNU General Public License for more details.
  14.  *
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include "types.h"
  21. #include "externs.h"
  22.  
  23. #define isign(i,j) ((j < 0) ? -i : i)
  24. #define mod2n(i,j) (i % j) + (j / 2) - isign((j/2),i)
  25.  
  26. #define amax1(a,b) ((a < b) ? b : a)
  27. #define amin1(a,b) ((a > b) ? b : a)
  28.  
  29. #define TEST(x) if ((DoTest & (1<<(x-1))) == (1<<(x-1)))
  30.  
  31. extern long int DoTest;
  32.  
  33. Void test();
  34. double sqrt();
  35. double exp();
  36.  
  37. Void kernel()
  38. {
  39.   Int i, l, k, il, ipntp, ipnt, lw, j, nl1, nl2, kx, ky, ip, i1, j1;
  40.   Int i2, j2, nr, nz;
  41.   Int ii, lb, j4, ink;
  42.   Int kn, jn, kb5i;
  43.   Float xtemp;
  44.  
  45.   test((Int)0);
  46.   TEST(1) {
  47.     for (l = 0; l < lp; l++)
  48.     for (k = 0; k < n; k++)
  49.       x[k]= q + y[k]*(r*z[k+10] + t*z[k+11]);
  50. }
  51.   test((Int)1);
  52.   TEST(2) {
  53.   for (l = 0; l < lp; l++) {
  54.     il= n;
  55.     ipntp= 0;
  56.   l_222:   ipnt= ipntp;
  57.     ipntp= ipntp+il;
  58.     il= il/2;
  59.     i= ipntp;
  60.     for (k = ipnt+1; k < ipntp; k += 2) {
  61.       i= i+1;
  62.       x[i]= x[k] - v[k]*x[k-1] - v[k+1]*x[k+1];
  63.     }
  64.     if( il>1) goto l_222;
  65.   }
  66. }
  67.   test((Int)2);
  68.   TEST(3) {
  69.   for (l = 0; l < lp; l++) {
  70.     q = 0.0;
  71.     for (k = 0; k < n; k++) {
  72.       q += z[k]*x[k];
  73.     }
  74.   }
  75. }
  76.   test((Int)3);
  77.   TEST(4) {
  78.   m= (1001-7)/2;
  79.   for (l = 0; l < lp; l++) {
  80.     for (k = 6; k < 1001; k += m) {
  81.       lw= k-6;
  82.       for (j = 4; j < n; j += 5) {
  83.     x[k-1]= x[k-1] - x[lw]*y[j];
  84.     lw= lw+1;
  85.       }
  86.       x[k-1]= y[4]*x[k-1];
  87.     }
  88.   }
  89. }
  90.   test((Int)4);
  91.   TEST(5) {
  92.   for (l = 0; l < lp; l++)
  93.     for (i = 1; i < n; i++)
  94.       x[i]= z[i]*(y[i] - x[i-1]);
  95. }
  96.   test((Int)5);
  97.   TEST(6) {
  98.   for (l = 0; l < lp; l++)
  99.     for (i = 1; i < n; i++)
  100.       for (k = 0; k < i; k++)
  101.     w[i] += b[i][k] * w[i-k-1];
  102. }
  103.   test((Int)6);
  104.   TEST(7) {
  105.   for (l = 0; l < lp; l++)
  106.     for (k = 0; k < n; k++)
  107.       x[k]=     u[k  ] + r*( z[k  ] + r*y[k  ]) +
  108.     t*( u[k+3] + r*( u[k+2] + r*u[k+1]) +
  109.        t*( u[k+6] + r*( u[k+5] + r*u[k+4])));
  110. }
  111.   test((Int)7);
  112.   TEST(8) {
  113.   for (l = 0; l < lp; l++) {
  114.     nl1 = 0;
  115.     nl2 = 1;
  116.     for (kx = 1; kx < 3; kx++) {
  117.       for (ky = 1; ky < n; ky++) {
  118.     du1[ky]= u1[kx][ky+1][nl1] - u1[kx][ky-1][nl1];
  119.     du2[ky]= u2[kx][ky+1][nl1] - u2[kx][ky-1][nl1];
  120.     du3[ky]= u3[kx][ky+1][nl1] - u3[kx][ky-1][nl1];
  121.     u1[kx][ky][nl2]= u1[kx][ky][nl1]+a11*du1[ky]+a12*du2[ky]+a13*du3[ky]+
  122.       sig*(u1[kx+1][ky][nl1]-2.*u1[kx][ky][nl1]+u1[kx-1][ky][nl1]);
  123.     u2[kx][ky][nl2]= u2[kx][ky][nl1]+a21*du1[ky]+a22*du2[ky]+a23*du3[ky]+
  124.       sig*(u2[kx+1][ky][nl1]-2.*u2[kx][ky][nl1]+u2[kx-1][ky][nl1]);
  125.     u3[kx][ky][nl2]= u3[kx][ky][nl1]+a31*du1[ky]+a32*du2[ky]+a33*du3[ky]+
  126.       sig*(u3[kx+1][ky][nl1]-2.*u3[kx][ky][nl1]+u3[kx-1][ky][nl1]);
  127.       }
  128.     }
  129.   }
  130. }
  131.   test((Int)8);
  132.   TEST(9) {
  133.   for (l = 0; l < lp; l++)
  134.     for (i = 0; i < n; i++)
  135.       px[ 0][i]= dm28*px[12][i] + dm27*px[11][i] + dm26*px[10][i] +
  136.     dm25*px[9][i] + dm24*px[8][i] + dm23*px[7][i] +
  137.       dm22*px[6][i] +  c0*(px[4][i] +      px[5][i])+ px[ 2][i];
  138. }
  139.   test((Int)9);
  140.   TEST(10) {
  141.   for (l = 0; l < lp; l++) {
  142.     for (i = 0; i < n; i++) {
  143.       ar      =      cx[4][i];
  144.       br      = ar - px[4][i];
  145.       px[4][i] = ar;
  146.       cr      = br - px[5][i];
  147.       px[5][i] = br;
  148.       ar      = cr - px[6][i];
  149.       px[6][i] = cr;
  150.       br      = ar - px[7][i];
  151.       px[7][i] = ar;
  152.       cr      = br - px[8][i];
  153.       px[8][i] = br;
  154.       ar      = cr - px[9][i];
  155.       px[9][i]= cr;
  156.       br      = ar - px[10][i];
  157.       px[10][i]= ar;
  158.       cr      = br - px[11][i];
  159.       px[11][i]= br;
  160.       px[13][i]= cr - px[12][i];
  161.       px[12][i]= cr;
  162.     }
  163.   }
  164. }
  165.   test((Int)10);
  166.   TEST(11) {
  167.   for (l = 0; l < lp; l++) {
  168.     x[0]= y[0];
  169.     for (k = 1; k < n; k++) {
  170.       x[k]= x[k-1] + y[k];
  171.     }
  172.   }
  173. }
  174.   test((Int)11);
  175.   TEST(12) {
  176.   for (l = 0; l < lp; l++) {
  177.     for (k = 0; k < n; k++) {
  178.       x[k]= y[k+1] - y[k];
  179.     }
  180.   }
  181. }
  182.   test((Int)12);
  183.   TEST(13) {
  184.   for (l = 0; l < lp; l++) {
  185.     for (ip = 0; ip < n; ip++) {
  186.       i1= p[0][ip];
  187.       j1= p[1][ip];
  188.       i1=        mod2n(i1,64);
  189.       j1=        mod2n(j1,64);
  190.       p[2][ip]= p[2][ip]  + b[i1][j1];
  191.       p[3][ip]= p[3][ip]  + c[i1][j1];
  192.       p[0][ip]= p[0][ip]  + p[2][ip];
  193.       p[1][ip]= p[1][ip]  + p[3][ip];
  194.       i2= p[0][ip];
  195.       j2= p[1][ip];
  196.       i2=            mod2n(i2,64) - 1;
  197.       j2=            mod2n(j2,64) - 1;
  198.       p[0][ip]= p[0][ip]  + y[i2+32];
  199.       p[1][ip]= p[1][ip]  + z[j2+32];
  200.       i2= i2       + e[i2+32];
  201.       j2= j2       + f[j2+32];
  202.       h[i2][j2]= h[i2][j2] + 1.0;
  203.     }
  204.   }
  205. }
  206.   test((Int)13);
  207.   TEST(14) {
  208.   for (l = 0; l < lp; l++) {
  209.     for (k = 0; k < n; k++) {
  210.       ix[k]= (Int) grd[k];
  211.       xi[k]= (Float) ix[k];
  212.       ex1[k]= ex [ ix[k]-1 ];
  213.       dex1[k]= dex  [ ix[k]-1 ];
  214.     }
  215.     for (k = 0; k < n; k++) {
  216.       vx[k]= vx[k] + ex1[k] + (xx[k] - xi[k])*dex1[k];
  217.       xx[k]= xx[k] + vx[k]  + flx;
  218.       ir[k]= xx[k];
  219.       rx[k]= xx[k] - ir[k];
  220.       ir[k]= mod2n(  ir[k],512) + 1;
  221.       xx[k]= rx[k] + ir[k];
  222.     }
  223.     for (k = 0; k < n; k++) {
  224.       rh[ir[k]-1]= rh[ir[k]-1] + 1.0 - rx[k];
  225.       rh[ir[k]]= rh[ir[k]] + rx[k];
  226.     }
  227.   }
  228. }
  229.   test((Int)14);
  230.   TEST(15) {
  231.   for (l = 0; l < lp; l++) {
  232.     nr= 7;
  233.     nz= n;
  234.     ar= 0.053;
  235.     br= 0.073;
  236.     for (j = 1; j < nr; j++) {
  237.       for (k = 1; k < nz; k++) {
  238.     if (j >= (nr-1)) {
  239.       vy[k][j] = 0.0;
  240.       continue;
  241.     }
  242.     if (vh[k][j+1] > vh[k][j]) {
  243.       t = ar;
  244.     } else {
  245.       t = br;
  246.     }
  247.     if (vf[k][j] < vf[k-1][j]) {
  248.           r= amax1(vh[k-1][j], vh[k-1][j+1]);
  249.           s= vf[k-1][j];
  250.     } else {
  251.           r= amax1( vh[k][j],   vh[k][j+1]);
  252.           s= vf[k][j];
  253.     }
  254.         vy[k][j]= sqrt((double) (vg[k][j]*vg[k][j] +r*r))*t/s;
  255.     if (k >= (nz-1)) {
  256.         vs[k][j] = 0.;
  257.         continue;
  258.     }
  259.         if (vf[k][j]  < vf[k][j-1]) {
  260.           r= amax1( vg[k][j-1], vg[k+1][j-1]);
  261.           s= vf[k][j-1];
  262.           t= br;
  263.     } else {
  264.           r= amax1( vg[k][j],   vg[k+1][j]);
  265.           s= vf[k][j];
  266.           t= ar;
  267.     }
  268.         vs[k][j]= sqrt((double) (vh[k][j]*vh[k][j] +r*r))*t/s;
  269.       }
  270.     }
  271.   }
  272. }
  273.   test((Int)15);
  274.   TEST(16) {
  275.   ii= n/3;
  276.   lb= ii+ii;
  277.   k2= 0;
  278.   k3= 0;
  279.   for (l = 0; l < lp; l++) {
  280.     m= 0;
  281.     i1= m;
  282. l_410:
  283.     j2= ((n+n)*m) + 1;
  284.     for (k = 0; k < n; k++) {
  285.       k2= k2+1;
  286.       j4= j2+k+k+1;
  287.       j5= zone[j4];
  288.       if (j5-n+1 < 0) goto l_420;
  289.       if (j5-n+1 == 0) goto l_475;
  290.       goto l_450;
  291. l_415:
  292.       if (j5-n+ii+1 < 0) goto l_430;
  293.       if (j5-n+ii+1 == 0) goto l_425;
  294.       goto l_425;
  295. l_420:
  296.       if (j5-n+lb+1 < 0) goto l_435;
  297.       if (j5-n+lb+1 == 0) goto l_415;
  298.       goto l_415;
  299. l_425:
  300.       if (plan[j5] - r < 0) goto l_445;
  301.       if (plan[j5] - r == 0) goto l_480;
  302.       goto l_440;
  303. l_430:
  304.       if (plan[j5] - s < 0) goto l_445;
  305.       if (plan[j5] - s == 0) goto l_480;
  306.       goto l_440;
  307. l_435:
  308.       if (plan[j5] - t < 0) goto l_445;
  309.       if (plan[j5] - t == 0) goto l_480;
  310.       goto l_440;
  311. l_440:
  312.       if (zone[j4-1]+1 < 0) goto l_455;
  313.       if (zone[j4-1]+1 == 0) goto l_485;
  314.       goto l_470;
  315. l_445:
  316.       if (zone[j4-1]+1 < 0) goto l_470;
  317.       if (zone[j4-1]+1 == 0) goto l_485;
  318.       goto l_455;
  319. l_450:
  320.       k3= k3+1;
  321.       xtemp = d[j5] - (d[j5-1]*((t-d[j5-2])*(t-d[j5-2])
  322.                 +(s-d[j5-3])*(s-d[j5-3])
  323.                 +(r-d[j5-4])*(r-d[j5-4])));
  324.       if (xtemp < 0) goto l_445;
  325.       if (xtemp == 0) goto l_480;
  326.       goto l_440;
  327. l_455:
  328.       m= m+1;
  329.       if (m < zone[0]) goto l_465;
  330.       if (m == zone[0]) goto l_465;
  331.       goto l_460;
  332. l_460:
  333.       m= 0;
  334. l_465:
  335.       if (i1 - m < 0) goto l_410;
  336.       if (i1 - m == 0) goto l_480;
  337.       goto l_410;
  338. l_470: ;
  339.     }
  340. l_475:
  341. l_480:
  342. l_485: ;
  343.   }
  344. }
  345.   test((Int)16);
  346.   TEST(17) {
  347.   for (l = 0; l < lp; l++) {
  348.     i= n-1;
  349.     j= 0;
  350.     ink= -1;
  351.     scale= 5./3.;
  352.     xnm= 1./3.;
  353.     e6= 1.03/3.07;
  354.     goto l_61;
  355. l_60:
  356.     e6= xnm*vsp[i]+vstp[i];
  357.     vxne[i]= e6;
  358.     xnm= e6;
  359.     ve3[i]= e6;
  360.     i= i+ink;
  361.     if (i == j) goto  l_62;
  362. l_61:
  363.     e3= xnm*vlr[i] +vlin[i];
  364.     xnei= vxne[i];
  365.     vxnd[i]= e6;
  366.     xnc= scale*e3;
  367.     if ( xnm >xnc) goto  l_60;
  368.     if ( xnei>xnc) goto  l_60;
  369.     ve3[i]= e3;
  370.     e6= e3+e3-xnm;
  371.     vxne[i]= e3+e3-xnei;
  372.     xnm= e6;
  373.     i= i+ink;
  374.     if ( i != j) goto l_61;
  375. l_62:
  376.    ;
  377.   }
  378. }
  379.   test((Int)17);
  380.   TEST(18) {
  381.   for (l = 0; l < lp; l++) {
  382.     t= 0.0037;
  383.     s= 0.0041;
  384.     kn= 6;
  385.     jn= n;
  386.     for (k = 1; k < kn; k++) {
  387.       for (j = 1; j < jn; j++) {
  388.     za[j][k]= (zp[j-1][k+1]+zq[j-1][k+1]-zp[j-1][k]-zq[j-1][k])
  389.       *(zr[j][k]+zr[j-1][k])/(zm[j-1][k]+zm[j-1][k+1]);
  390.     zb[j][k]= (zp[j-1][k]+zq[j-1][k]-zp[j][k]-zq[j][k])
  391.       *(zr[j][k]+zr[j][k-1])/(zm[j][k]+zm[j-1][k]);
  392.       }
  393.     }
  394.     for (k = 1; k < kn; k++) {
  395.       for (j = 1; j < jn; j++) {
  396. #ifdef CRAY2
  397.     temp1 =  za[j][k]   *(zz[j][k]-zz[j+1][k]);
  398.     temp2 = -za[j-1][k] *(zz[j][k]-zz[j-1][k]);
  399.     temp3 = -zb[j][k]   *(zz[j][k]-zz[j][k-1]);
  400.     temp4 =  zb[j][k+1] *(zz[j][k]-zz[j][k+1]);
  401.     zu[j][k] += s*(temp1+temp2+temp3+temp4);
  402.         temp1 =za[j][k]*(zr[j][k]-zr[j+1][k]);
  403.     temp2 = -za[j-1][k] *(zr[j][k]-zr[j-1][k]);
  404.     temp3 = -zb[j][k]   *(zr[j][k]-zr[j][k-1]);
  405.     temp4 = zb[j][k+1] *(zr[j][k]-zr[j][k+1]);
  406.         zv[j][k] += s*(temp1+temp2+temp3+temp4);
  407. #else    
  408.     zu[j][k]= zu[j][k]+s*(za[j][k]*(zz[j][k]-zz[j+1][k])
  409.                   -za[j-1][k] *(zz[j][k]-zz[j-1][k])
  410.                   -zb[j][k]   *(zz[j][k]-zz[j][k-1])
  411.                   +zb[j][k+1] *(zz[j][k]-zz[j][k+1]));
  412.     zv[j][k]= zv[j][k]+s*(za[j][k]*(zr[j][k]-zr[j+1][k])
  413.                   -za[j-1][k] *(zr[j][k]-zr[j-1][k])
  414.                   -zb[j][k]   *(zr[j][k]-zr[j][k-1])
  415.                   +zb[j][k+1] *(zr[j][k]-zr[j][k+1]));
  416. #endif
  417.       }
  418.     }
  419.     for (k = 1; k < kn; k++) {
  420.       for (j = 1; j < jn; j++) {
  421.           zr[j][k]= zr[j][k]+t*zu[j][k];
  422.           zz[j][k]= zz[j][k]+t*zv[j][k];
  423.     }
  424.     }
  425.   }
  426. }
  427.   test((Int)18);
  428.   TEST(19) {
  429.   kb5i= 0;
  430.   for (l = 0; l < lp; l++) {
  431.     for (k = 0; k < n; k++) {
  432.       b5[k+kb5i]= sa[k] +stb5*sb[k];
  433.       stb5= b5[k+kb5i] -stb5;
  434.     }
  435.     for (i = 0; i < n; i++) {
  436.       k= n-i+1;
  437.       b5[k+kb5i]= sa[k] +stb5*sb[k];
  438.       stb5= b5[k+kb5i] -stb5;
  439.     }
  440.   }
  441. }
  442.   test((Int)19);
  443.   TEST(20) {
  444.   for (l = 0; l < lp; l++) {
  445.     for (k = 0; k < n; k++) {
  446.       di= y[k]-g[k]/( xx[k]+dk);
  447.       dn= 0.2;
  448.       if ( di != 0.0) dn= amax1( 0.1,amin1( z[k]/di, 0.2));
  449.       x[k]= ((w[k]+v[k]*dn)* xx[k]+u[k])/(vx[k]+v[k]*dn);
  450.       xx[k+1]= (x[k]- xx[k])*dn+ xx[k];
  451.     }
  452.   }
  453. }
  454.   test((Int)20);
  455.   TEST(21) {
  456.   for (l = 0; l < lp; l++)
  457.     for (i = 0; i < n; i++)
  458.       for (k = 0; k < n; k++)
  459.     for (j = 0; j < n; j++)
  460.       px[i][j]= px[i][j] +vy[i][k] * cx[k][j];
  461. }
  462.   test((Int)21);
  463.   TEST(22) {
  464.   expmax= 20.0;
  465.   u[n-1]= 0.99*expmax*v[n-1];
  466.   for (l = 0; l < lp; l++) {
  467.     for (k = 0; k < n; k++) {
  468.       y[k]= u[k]/v[k];
  469.       w[k]= x[k]/( exp((double) y[k]) -1.0);
  470.     }
  471.   }
  472. }
  473.   test((Int)22);
  474.   TEST(23) {
  475.   for (l = 0; l < lp; l++) {
  476.     for (j = 1; j < 6; j++) {
  477.       for (k = 1; k < n; k++) {
  478.     qa= za[k][j+1]*zr[k][j] +za[k][j-1]*zb[k][j] +
  479.       za[k+1][j]*zu[k][j] +za[k-1][j]*zv[k][j] +zz[k][j];
  480.     za[k][j]= za[k][j] +.175*(qa -za[k][j]);
  481.       }
  482.     }
  483.   }
  484. }
  485.   test((Int)23);
  486.   TEST(24) {
  487.   x[ n/2]= -1.0e+10;
  488.   for (l = 0; l < lp; l++) {
  489.     m= 1;
  490.     for (k = 1; k < n; k++) {
  491.       if ( x[k] < x[m])  m= k;
  492.     }
  493.   }
  494. }
  495.   test((Int)24);
  496.   return;
  497. }
  498.