home *** CD-ROM | disk | FTP | other *** search
/ NetNews Usenet Archive 1992 #20 / NN_1992_20.iso / spool / comp / benchmar / 1381 < prev    next >
Encoding:
Text File  |  1992-09-07  |  40.5 KB  |  1,408 lines
  1. Path: sparky!uunet!spool.mu.edu!darwin.sura.net!jvnc.net!netnews.upenn.edu!msuinfo!vthnw.cvm.msu.edu!eoq
  2. From: eoq@vthnw.cvm.msu.edu (Edward Quillen)
  3. Newsgroups: comp.benchmarks
  4. Subject: C source for Livermore Loops benchmark
  5. Message-ID: <eoq.40.0@vthnw.cvm.msu.edu>
  6. Date: 6 Sep 92 01:42:11 GMT
  7. Sender: news@msuinfo.cl.msu.edu
  8. Organization: Veterinary Teaching Hospital, Michigan State Univ.
  9. Lines: 1397
  10.  
  11.  
  12.  
  13. /* Livermore Loops coded in C        Latest File Modification  20 Feb 87.
  14.  *     SUBROUTINE KERNEL   replaces the Fortran routine in LFK Test program.
  15.  ************************************************************************
  16.  *                                                                      *
  17.  *            KERNEL     executes 24 samples of "C" computation         *
  18.  *                                                                      *
  19.  ************************************************************************
  20.  *                                                                      *
  21.  *     L. L. N. L.   " C "   K E R N E L S:   M F L O P S               *
  22.  *                                                                      *
  23.  *     These kernels measure   " C "   numerical computation            *
  24.  *     rates for  a  spectrum  of  cpu-limited computational            *
  25.  *     structures or benchmarks.   Mathematical  through-put            *
  26.  *     is measured  in  units  of millions of floating-point            *
  27.  *     operations executed per second, called Megaflops/sec.            *
  28.  *                                                                      *
  29.  *     Fonzi's Law: There is not now and there never will be a language *
  30.  *                  in which it is the least bit difficult to write     *
  31.  *                  bad programs.                                       *
  32.  *                                                    F.H.MCMAHON  1972 *
  33.  ************************************************************************
  34.  *    Originally from  Greg Astfalk, AT&T, P.O.Box 900, Princeton, NJ. 08540.
  35.  *               by way of Frank McMahon (LLNL).
  36.  *    Changes made to correct many array subscripting problems,
  37.  *      make more readable (added #define's), include the original
  38.  *      FORTRAN versions of the runs as comments, and make more
  39.  *      portable by Kelly O'Hair (LLNL) and Chuck Rasbold (LLNL).
  40.  *
  41.  */
  42.  
  43. #include <stdio.h>
  44. #include <math.h>
  45.  
  46. /* Type-specifiers for function declarations */
  47. #ifdef CRAY
  48. #define CALLED_FROM_FORTRAN fortran
  49. #define FORTRAN_FUNCTION    fortran
  50. #else
  51. #define CALLED_FROM_FORTRAN
  52. #define FORTRAN_FUNCTION    extern
  53. #endif
  54.  
  55. /* These are the names to be used for the functions */
  56. #ifdef CRAY
  57. #define KERNEL kernel
  58. #define TEST   test
  59. #define SPACE  space
  60. #else
  61. #define KERNEL kernel_
  62. #define TEST   test_
  63. #define SPACE  space_
  64. #endif
  65.  
  66. /* Type-specifiers for the structs that map to common */
  67. #ifdef CRAY
  68. #define COMMON_BLOCK fortran
  69. #else
  70. #define COMMON_BLOCK extern
  71. #endif
  72.  
  73. /* Names of structs (or common blocks) */
  74. #ifdef CRAY
  75. #define ALPHA  alpha
  76. #define BETA   beta
  77. #define SPACES spaces
  78. #define SPACER spacer
  79. #define SPACE0 space0
  80. #define SPACEI spacei
  81. #define ISPACE ispace
  82. #define SPACE1 space1
  83. #define SPACE2 space2
  84. #else
  85. #define ALPHA  alpha_
  86. #define BETA   beta_
  87. #define SPACES spaces_
  88. #define SPACER spacer_
  89. #define SPACE0 space0_
  90. #define SPACEI spacei_
  91. #define ISPACE ispace_
  92. #define SPACE1 space1_
  93. #define SPACE2 space2_
  94. #endif
  95.  
  96. /* Declare the two fortran functions called */
  97. FORTRAN_FUNCTION TEST();
  98. FORTRAN_FUNCTION SPACE();
  99.  
  100. /* Define the structs or COMMON BLOCKS */
  101.  
  102. COMMON_BLOCK struct {
  103.     long Mk;
  104.     long Ik;
  105.     long Ml;
  106.     long Il;
  107.     long Mr;
  108.     long Jr;
  109.     long Npfs[47][3][8];
  110. } ALPHA ;
  111. #define mk     ALPHA.Mk
  112. #define ik     ALPHA.Ik
  113. #define ml     ALPHA.Ml
  114. #define il     ALPHA.Il
  115. #define mr     ALPHA.Mr
  116. #define jr     ALPHA.Jr
  117. #define npfs   ALPHA.Npfs
  118.  
  119. COMMON_BLOCK struct {
  120.     double Tic;
  121.     double Times[47][3][8];
  122.     double See[3][8][3][5];
  123.     double Terrs[47][3][8];
  124.     double Csums[47][3][8];
  125.     double Fopn[47][3][8];
  126.     double Dos[47][3][8];
  127. } BETA ;
  128. #define tic     BETA.Tic
  129. #define times   BETA.Times
  130. #define see     BETA.See
  131. #define terrs   BETA.Terrs
  132. #define csums   BETA.Csums
  133. #define fopn    BETA.Fopn
  134. #define dos     BETA.Dos
  135.  
  136. COMMON_BLOCK struct {
  137.     long Ion;
  138.     long J5;
  139.     long K2;
  140.     long K3;
  141.     long Loop;
  142.     long M;
  143.     long Kr;
  144.     long It;
  145.     long N13h;
  146.     long Ibuf;
  147.     long N;
  148.     long N1;
  149.     long N2;
  150.     long N13;
  151.     long N213;
  152.     long N813;
  153.     long N14;
  154.     long N16;
  155.     long N416;
  156.     long N21;
  157.     long Nt1;
  158.     long Nt2;
  159. } SPACES ;
  160. #define  ion    SPACES.Ion
  161. #define  j5     SPACES.J5
  162. #define  k2     SPACES.K2
  163. #define  k3     SPACES.K3
  164. #define  loop   SPACES.Loop
  165. #define  m      SPACES.M
  166. #define  kr     SPACES.Kr
  167. #define  it     SPACES.It
  168. #define  n13h   SPACES.N13h
  169. #define  ibuf   SPACES.Ibuf
  170. #define  n      SPACES.N
  171. #define  n1     SPACES.N1
  172. #define  n2     SPACES.N2
  173. #define  n13    SPACES.N13
  174. #define  n213   SPACES.N213
  175. #define  n813   SPACES.N813
  176. #define  n14    SPACES.N14
  177. #define  n16    SPACES.N16
  178. #define  n416   SPACES.N416
  179. #define  n21    SPACES.N21
  180. #define  nt1    SPACES.Nt1
  181. #define  nt2    SPACES.Nt2
  182.  
  183. COMMON_BLOCK struct {
  184.     double A11;
  185.     double A12;
  186.     double A13;
  187.     double A21;
  188.     double A22;
  189.     double A23;
  190.     double A31;
  191.     double A32;
  192.     double A33;
  193.     double Ar;
  194.     double Br;
  195.     double C0;
  196.     double Cr;
  197.     double Di;
  198.     double Dk;
  199.     double Dm22;
  200.     double Dm23;
  201.     double Dm24;
  202.     double Dm25;
  203.     double Dm26;
  204.     double Dm27;
  205.     double Dm28;
  206.     double Dn;
  207.     double E3;
  208.     double E6;
  209.     double Expmax;
  210.     double Flx;
  211.     double Q;
  212.     double Qa;
  213.     double R;
  214.     double Ri;
  215.     double S;
  216.     double Scale;
  217.     double Sig;
  218.     double Stb5;
  219.     double T;
  220.     double Xnc;
  221.     double Xnei;
  222.     double Xnm;
  223. } SPACER ;
  224. #define  a11     SPACER.A11
  225. #define  a12     SPACER.A12
  226. #define  a13     SPACER.A13
  227. #define  a21     SPACER.A21
  228. #define  a22     SPACER.A22
  229. #define  a23     SPACER.A23
  230. #define  a31     SPACER.A31
  231. #define  a32     SPACER.A32
  232. #define  a33     SPACER.A33
  233. #define  ar      SPACER.Ar
  234. #define  br      SPACER.Br
  235. #define  c0      SPACER.C0
  236. #define  cr      SPACER.Cr
  237. #define  di      SPACER.Di
  238. #define  dk      SPACER.Dk
  239. #define  dm22    SPACER.Dm22
  240. #define  dm23    SPACER.Dm23
  241. #define  dm24    SPACER.Dm24
  242. #define  dm25    SPACER.Dm25
  243. #define  dm26    SPACER.Dm26
  244. #define  dm27    SPACER.Dm27
  245. #define  dm28    SPACER.Dm28
  246. #define  dn      SPACER.Dn
  247. #define  e3      SPACER.E3
  248. #define  e6      SPACER.E6
  249. #define  expmax  SPACER.Expmax
  250. #define  flx     SPACER.Flx
  251. #define  q       SPACER.Q
  252. #define  qa      SPACER.Qa
  253. #define  r       SPACER.R
  254. #define  ri      SPACER.Ri
  255. #define  s       SPACER.S
  256. #define  scale   SPACER.Scale
  257. #define  sig     SPACER.Sig
  258. #define  stb5    SPACER.Stb5
  259. #define  t       SPACER.T
  260. #define  xnc     SPACER.Xnc
  261. #define  xnei    SPACER.Xnei
  262. #define  xnm     SPACER.Xnm
  263.  
  264. COMMON_BLOCK struct {
  265.     double Time[47];
  266.     double Csum[47];
  267.     double Ww[47];
  268.     double Wt[47];
  269.     double Ticks;
  270.     double Fr[9];
  271.     double Terr1[47];
  272.     double Sumw[7];
  273.     double Start;
  274.     double Skale[47];
  275.     double Bias[47];
  276.     double Ws[95];
  277.     double Runs[47];
  278.     double Flopn[47];
  279.     long Iq[7];
  280.     long Npf;
  281.     long Npfs1[47];
  282. } SPACE0 ;
  283. #define  time    SPACE0.Time
  284. #define  csum    SPACE0.Csum
  285. #define  ww      SPACE0.Ww
  286. #define  wt      SPACE0.Wt
  287. #define  ticks   SPACE0.Ticks
  288. #define  fr      SPACE0.Fr
  289. #define  terr1   SPACE0.Terr1
  290. #define  sumw    SPACE0.Sumw
  291. #define  start   SPACE0.Start
  292. #define  skale   SPACE0.Skale
  293. #define  bias    SPACE0.Bias
  294. #define  ws      SPACE0.Ws
  295. #define  runs    SPACE0.Runs
  296. #define  flopn   SPACE0.Flopn
  297. #define  iq      SPACE0.Iq
  298. #define  npf     SPACE0.Npf
  299. #define  npfs1   SPACE0.Npfs1
  300.  
  301. COMMON_BLOCK struct {
  302.     double Wtp[3];
  303.     long Mult[3];
  304.     long Ispan[3][47];
  305.     long Ipass[3][47];
  306. } SPACEI ;
  307. #define wtp    SPACEI.Wtp
  308. #define mult   SPACEI.Mult
  309. #define ispan  SPACEI.Ispan
  310. #define ipass  SPACEI.Ipass
  311.  
  312. COMMON_BLOCK struct {
  313.     long E[96];
  314.     long F[96];
  315.     long Ix[1001];
  316.     long Ir[1001];
  317.     long Zone[300];
  318. } ISPACE ;
  319. #define e    ISPACE.E
  320. #define f    ISPACE.F
  321. #define ix   ISPACE.Ix
  322. #define ir   ISPACE.Ir
  323. #define zone ISPACE.Zone
  324.  
  325. COMMON_BLOCK struct {
  326.     double U[1001];
  327.     double V[1001];
  328.     double W[1001];
  329.     double X[1001];
  330.     double Y[1001];
  331.     double Z[1001];
  332.     double G[1001];
  333.     double Du1[101];
  334.     double Du2[101];
  335.     double Du3[101];
  336.     double Grd[1001];
  337.     double Dex[1001];
  338.     double Xi[1001];
  339.     double Ex[1001];
  340.     double Ex1[1001];
  341.     double Dex1[1001];
  342.     double Vx[1001];
  343.     double Xx[1001];
  344.     double Rx[1001];
  345.     double Rh[2048];
  346.     double Vsp[101];
  347.     double Vstp[101];
  348.     double Vxne[101];
  349.     double Vxnd[101];
  350.     double Ve3[101];
  351.     double Vlr[101];
  352.     double Vlin[101];
  353.     double B5[101];
  354.     double Plan[300];
  355.     double D[300];
  356.     double Sa[101];
  357.     double Sb[101];
  358. } SPACE1 ;
  359. #define  u    SPACE1.U
  360. #define  v    SPACE1.V
  361. #define  w    SPACE1.W
  362. #define  x    SPACE1.X
  363. #define  y    SPACE1.Y
  364. #define  z    SPACE1.Z
  365. #define  g    SPACE1.G
  366. #define  du1  SPACE1.Du1
  367. #define  du2  SPACE1.Du2
  368. #define  du3  SPACE1.Du3
  369. #define  grd  SPACE1.Grd
  370. #define  dex  SPACE1.Dex
  371. #define  xi   SPACE1.Xi
  372. #define  ex   SPACE1.Ex
  373. #define  ex1  SPACE1.Ex1
  374. #define  dex1 SPACE1.Dex1
  375. #define  vx   SPACE1.Vx
  376. #define  xx   SPACE1.Xx
  377. #define  rx   SPACE1.Rx
  378. #define  rh   SPACE1.Rh
  379. #define  vsp  SPACE1.Vsp
  380. #define  vstp SPACE1.Vstp
  381. #define  vxne SPACE1.Vxne
  382. #define  vxnd SPACE1.Vxnd
  383. #define  ve3  SPACE1.Ve3
  384. #define  vlr  SPACE1.Vlr
  385. #define  vlin SPACE1.Vlin
  386. #define  b5   SPACE1.B5
  387. #define  plan SPACE1.Plan
  388. #define  d    SPACE1.D
  389. #define  sa   SPACE1.Sa
  390. #define  sb   SPACE1.Sb
  391.  
  392. COMMON_BLOCK struct {
  393.     double P[512][4];
  394.     double Px[101][25];
  395.     double Cx[101][25];
  396.     double Vy[25][101];
  397.     double Vh[7][101];
  398.     double Vf[7][101];
  399.     double Vg[7][101];
  400.     double Vs[7][101];
  401.     double Za[7][101];
  402.     double Zp[7][101];
  403.     double Zq[7][101];
  404.     double Zr[7][101];
  405.     double Zm[7][101];
  406.     double Zb[7][101];
  407.     double Zu[7][101];
  408.     double Zv[7][101];
  409.     double Zz[7][101];
  410.     double B[64][64];
  411.     double C[64][64];
  412.     double H[64][64];
  413.     double U1[2][101][5];
  414.     double U2[2][101][5];
  415.     double U3[2][101][5];
  416. } SPACE2 ;
  417. #define  p       SPACE2.P
  418. #define  px      SPACE2.Px
  419. #define  cx      SPACE2.Cx
  420. #define  vy      SPACE2.Vy
  421. #define  vh      SPACE2.Vh
  422. #define  vf      SPACE2.Vf
  423. #define  vg      SPACE2.Vg
  424. #define  vs      SPACE2.Vs
  425. #define  za      SPACE2.Za
  426. #define  zp      SPACE2.Zp
  427. #define  zq      SPACE2.Zq
  428. #define  zr      SPACE2.Zr
  429. #define  zm      SPACE2.Zm
  430. #define  zb      SPACE2.Zb
  431. #define  zu      SPACE2.Zu
  432. #define  zv      SPACE2.Zv
  433. #define  zz      SPACE2.Zz
  434. #define  b       SPACE2.B
  435. #define  c       SPACE2.C
  436. #define  h       SPACE2.H
  437. #define  u1      SPACE2.U1
  438. #define  u2      SPACE2.U2
  439. #define  u3      SPACE2.U3
  440.  
  441. /* KERNEL routine */
  442.  
  443. CALLED_FROM_FORTRAN KERNEL()
  444. {
  445.  
  446. #pragma nodyneqv
  447. #pragma o=i
  448.  
  449.     long argument , k , l , ipnt , ipntp , i;
  450.     long lw , j , nl1 , nl2 , kx , ky , ip , kn;
  451.     long i1 , j1 , i2 , j2 , nz , ink , jn , kb5i;
  452.     long ii , lb , j4 , ng;
  453.     double tmp , temp;
  454.  
  455.     /*
  456.      *   Prologue
  457.      */
  458.  
  459.     SPACE();
  460.     argument = 0;
  461.     TEST( &argument );
  462.  
  463.     /*
  464.      *******************************************************************
  465.      *   Kernel 1 -- hydro fragment
  466.      *******************************************************************
  467.      *       DO 1 L = 1,Loop
  468.      *       DO 1 k = 1,n
  469.      *  1       X(k)= Q + Y(k)*(R*ZX(k+10) + T*ZX(k+11))
  470.      */
  471.  
  472.     for ( l=1 ; l<=loop ; l++ ) {
  473.         for ( k=0 ; k<n ; k++ ) {
  474.             x[k] = q + y[k]*( r*z[k+10] + t*z[k+11] );
  475.         }
  476.     }
  477.     argument = 1;
  478.     TEST( &argument );
  479.  
  480.     /*
  481.      *******************************************************************
  482.      *   Kernel 2 -- ICCG excerpt (Incomplete Cholesky Conjugate Gradient)
  483.      *******************************************************************
  484.      *    DO 200  L= 1,Loop
  485.      *        II= n
  486.      *     IPNTP= 0
  487.      *222   IPNT= IPNTP
  488.      *     IPNTP= IPNTP+II
  489.      *        II= II/2
  490.      *         i= IPNTP
  491.      CDIR$ IVDEP
  492.      *    DO 2 k= IPNT+2,IPNTP,2
  493.      *         i= i+1
  494.      *  2   X(i)= X(k) - V(k)*X(k-1) - V(k+1)*X(k+1)
  495.      *        IF( II.GT.1) GO TO 222
  496.      *200 CONTINUE
  497.      */
  498.  
  499.     for ( l=1 ; l<=loop ; l++ ) {
  500.         ii = n;
  501.         ipntp = 0;
  502.         do {
  503.             ipnt = ipntp;
  504.             ipntp += ii;
  505.             ii /= 2;
  506.             i = ipntp - 1;
  507. #pragma nohazard
  508.             for ( k=ipnt+1 ; k<ipntp ; k=k+2 ) {
  509.                 i++;
  510.                 x[i] = x[k] - v[k  ]*x[k-1] - v[k+1]*x[k+1];
  511.             }
  512.         } while ( ii>0 );
  513.     }
  514.     argument = 2;
  515.     TEST( &argument );
  516.  
  517.     /*
  518.      *******************************************************************
  519.      *   Kernel 3 -- inner product
  520.      *******************************************************************
  521.      *    DO 3 L= 1,Loop
  522.      *         Q= 0.0
  523.      *    DO 3 k= 1,n
  524.      *  3      Q= Q + Z(k)*X(k)
  525.      */
  526.  
  527.     for ( l=1 ; l<=loop ; l++ ) {
  528.         q = 0.0;
  529.         for ( k=0 ; k<n ; k++ ) {
  530.             q += z[k]*x[k];
  531.         }
  532.     }
  533.     argument = 3;
  534.     TEST( &argument );
  535.  
  536.     /*
  537.      *******************************************************************
  538.      *   Kernel 4 -- banded linear equations
  539.      *******************************************************************
  540.      *            m= (1001-7)/2
  541.      *    DO 444  L= 1,Loop
  542.      *    DO 444  k= 7,1001,m
  543.      *           lw= k-6
  544.      *         temp= X(k-1)
  545.      CDIR$ IVDEP
  546.      *    DO   4  j= 5,n,5
  547.      *       temp  = temp   - XZ(lw)*Y(j)
  548.      *  4        lw= lw+1
  549.      *       X(k-1)= Y(5)*temp
  550.      *444 CONTINUE
  551.      */
  552.  
  553.     m = ( 1001-7 )/2;
  554.     for ( l=1 ; l<=loop ; l++ ) {
  555.         for ( k=6 ; k<1001 ; k=k+m ) {
  556.             lw = k - 6;
  557.             temp = x[k-1];
  558. #pragma nohazard
  559.             for ( j=4 ; j<n ; j=j+5 ) {
  560.                 temp -= x[lw]*y[j];
  561.                 lw++;
  562.             }
  563.             x[k-1] = y[4]*temp;
  564.         }
  565.     }
  566.     argument = 4;
  567.     TEST( &argument );
  568.  
  569.     /*
  570.      *******************************************************************
  571.      *   Kernel 5 -- tri-diagonal elimination, below diagonal
  572.      *******************************************************************
  573.      *    DO 5 L = 1,Loop
  574.      *    DO 5 i = 2,n
  575.      *  5    X(i)= Z(i)*(Y(i) - X(i-1))
  576.      */
  577.  
  578.     for ( l=1 ; l<=loop ; l++ ) {
  579.         for ( i=1 ; i<n ; i++ ) {
  580.             x[i] = z[i]*( y[i] - x[i-1] );
  581.         }
  582.     }
  583.     argument = 5;
  584.     TEST( &argument );
  585.  
  586.     /*
  587.      *******************************************************************
  588.      *   Kernel 6 -- general linear recurrence equations
  589.      *******************************************************************
  590.      *    DO  6  L= 1,Loop
  591.      *    DO  6  i= 2,n
  592.      *    DO  6  k= 1,i-1
  593.      *        W(i)= W(i)  + B(i,k) * W(i-k)
  594.      *  6 CONTINUE
  595.      */
  596.  
  597.     for ( l=1 ; l<=loop ; l++ ) {
  598.         for ( i=1 ; i<n ; i++ ) {
  599.             for ( k=0 ; k<i ; k++ ) {
  600.                 w[i] += b[k][i] * w[(i-k)-1];
  601.             }
  602.         }
  603.     }
  604.     argument = 6;
  605.     TEST( &argument );
  606.  
  607.     /*
  608.      *******************************************************************
  609.      *   Kernel 7 -- equation of state fragment
  610.      *******************************************************************
  611.      *    DO 7 L= 1,Loop
  612.      *    DO 7 k= 1,n
  613.      *      X(k)=     U(k  ) + R*( Z(k  ) + R*Y(k  )) +
  614.      *   .        T*( U(k+3) + R*( U(k+2) + R*U(k+1)) +
  615.      *   .        T*( U(k+6) + R*( U(k+5) + R*U(k+4))))
  616.      *  7 CONTINUE
  617.      */
  618.  
  619.     for ( l=1 ; l<=loop ; l++ ) {
  620. #pragma nohazard
  621.         for ( k=0 ; k<n ; k++ ) {
  622.             x[k] = u[k] + r*( z[k] + r*y[k] ) +
  623.                    t*( u[k+3] + r*( u[k+2] + r*u[k+1] ) +
  624.                       t*( u[k+6] + r*( u[k+5] + r*u[k+4] ) ) );
  625.         }
  626.     }
  627.     argument = 7;
  628.     TEST( &argument );
  629.  
  630.     /*
  631.      *******************************************************************
  632.      *   Kernel 8 -- ADI integration
  633.      *******************************************************************
  634.      *    DO  8      L = 1,Loop
  635.      *             nl1 = 1
  636.      *             nl2 = 2
  637.      *    DO  8     kx = 2,3
  638.      CDIR$ IVDEP
  639.      *    DO  8     ky = 2,n
  640.      *          DU1(ky)=U1(kx,ky+1,nl1)  -  U1(kx,ky-1,nl1)
  641.      *          DU2(ky)=U2(kx,ky+1,nl1)  -  U2(kx,ky-1,nl1)
  642.      *          DU3(ky)=U3(kx,ky+1,nl1)  -  U3(kx,ky-1,nl1)
  643.      *    U1(kx,ky,nl2)=U1(kx,ky,nl1) +A11*DU1(ky) +A12*DU2(ky) +A13*DU3(ky)
  644.      *   .       + SIG*(U1(kx+1,ky,nl1) -2.*U1(kx,ky,nl1) +U1(kx-1,ky,nl1))
  645.      *    U2(kx,ky,nl2)=U2(kx,ky,nl1) +A21*DU1(ky) +A22*DU2(ky) +A23*DU3(ky)
  646.      *   .       + SIG*(U2(kx+1,ky,nl1) -2.*U2(kx,ky,nl1) +U2(kx-1,ky,nl1))
  647.      *    U3(kx,ky,nl2)=U3(kx,ky,nl1) +A31*DU1(ky) +A32*DU2(ky) +A33*DU3(ky)
  648.      *   .       + SIG*(U3(kx+1,ky,nl1) -2.*U3(kx,ky,nl1) +U3(kx-1,ky,nl1))
  649.      *  8 CONTINUE
  650.      */
  651.  
  652.     for ( l=1 ; l<=loop ; l++ ) {
  653.         nl1 = 0;
  654.         nl2 = 1;
  655.         for ( kx=1 ; kx<3 ; kx++ ){
  656. #pragma nohazard
  657.             for ( ky=1 ; ky<n ; ky++ ) {
  658.                du1[ky] = u1[nl1][ky+1][kx] - u1[nl1][ky-1][kx];
  659.                du2[ky] = u2[nl1][ky+1][kx] - u2[nl1][ky-1][kx];
  660.                du3[ky] = u3[nl1][ky+1][kx] - u3[nl1][ky-1][kx];
  661.                u1[nl2][ky][kx]=
  662.                   u1[nl1][ky][kx]+a11*du1[ky]+a12*du2[ky]+a13*du3[ky] + sig*
  663.                      (u1[nl1][ky][kx+1]-2.0*u1[nl1][ky][kx]+u1[nl1][ky][kx-1]);
  664.                u2[nl2][ky][kx]=
  665.                   u2[nl1][ky][kx]+a21*du1[ky]+a22*du2[ky]+a23*du3[ky] + sig*
  666.                      (u2[nl1][ky][kx+1]-2.0*u2[nl1][ky][kx]+u2[nl1][ky][kx-1]);
  667.                u3[nl2][ky][kx]=
  668.                   u3[nl1][ky][kx]+a31*du1[ky]+a32*du2[ky]+a33*du3[ky] + sig*
  669.                      (u3[nl1][ky][kx+1]-2.0*u3[nl1][ky][kx]+u3[nl1][ky][kx-1]);
  670.             }
  671.         }
  672.     }
  673.     argument = 8;
  674.     TEST( &argument );
  675.  
  676.     /*
  677.      *******************************************************************
  678.      *   Kernel 9 -- integrate predictors
  679.      *******************************************************************
  680.      *    DO 9  L = 1,Loop
  681.      *    DO 9  i = 1,n
  682.      *    PX( 1,i)= DM28*PX(13,i) + DM27*PX(12,i) + DM26*PX(11,i) +
  683.      *   .          DM25*PX(10,i) + DM24*PX( 9,i) + DM23*PX( 8,i) +
  684.      *   .          DM22*PX( 7,i) +  C0*(PX( 5,i) +      PX( 6,i))+ PX( 3,i)
  685.      *  9 CONTINUE
  686.      */
  687.  
  688.     for ( l=1 ; l<=loop ; l++ ) {
  689.         for ( i=0 ; i<n ; i++ ) {
  690.             px[i][0] = dm28*px[i][12] + dm27*px[i][11] + dm26*px[i][10] +
  691.                        dm25*px[i][ 9] + dm24*px[i][ 8] + dm23*px[i][ 7] +
  692.                        dm22*px[i][ 6] + c0*( px[i][ 4] + px[i][ 5]) + px[i][ 2];
  693.         }
  694.     }
  695.     argument = 9;
  696.     TEST( &argument );
  697.  
  698.     /*
  699.      *******************************************************************
  700.      *   Kernel 10 -- difference predictors
  701.      *******************************************************************
  702.      *    DO 10  L= 1,Loop
  703.      *    DO 10  i= 1,n
  704.      *    AR      =      CX(5,i)
  705.      *    BR      = AR - PX(5,i)
  706.      *    PX(5,i) = AR
  707.      *    CR      = BR - PX(6,i)
  708.      *    PX(6,i) = BR
  709.      *    AR      = CR - PX(7,i)
  710.      *    PX(7,i) = CR
  711.      *    BR      = AR - PX(8,i)
  712.      *    PX(8,i) = AR
  713.      *    CR      = BR - PX(9,i)
  714.      *    PX(9,i) = BR
  715.      *    AR      = CR - PX(10,i)
  716.      *    PX(10,i)= CR
  717.      *    BR      = AR - PX(11,i)
  718.      *    PX(11,i)= AR
  719.      *    CR      = BR - PX(12,i)
  720.      *    PX(12,i)= BR
  721.      *    PX(14,i)= CR - PX(13,i)
  722.      *    PX(13,i)= CR
  723.      * 10 CONTINUE
  724.      */
  725.  
  726.     for ( l=1 ; l<=loop ; l++ ) {
  727.         for ( i=0 ; i<n ; i++ ) {
  728.             ar        =      cx[i][ 4];
  729.             br        = ar - px[i][ 4];
  730.             px[i][ 4] = ar;
  731.             cr        = br - px[i][ 5];
  732.             px[i][ 5] = br;
  733.             ar        = cr - px[i][ 6];
  734.             px[i][ 6] = cr;
  735.             br        = ar - px[i][ 7];
  736.             px[i][ 7] = ar;
  737.             cr        = br - px[i][ 8];
  738.             px[i][ 8] = br;
  739.             ar        = cr - px[i][ 9];
  740.             px[i][ 9] = cr;
  741.             br        = ar - px[i][10];
  742.             px[i][10] = ar;
  743.             cr        = br - px[i][11];
  744.             px[i][11] = br;
  745.             px[i][13] = cr - px[i][12];
  746.             px[i][12] = cr;
  747.         }
  748.     }
  749.     argument = 10;
  750.     TEST( &argument );
  751.  
  752.     /*
  753.      *******************************************************************
  754.      *   Kernel 11 -- first sum
  755.      *******************************************************************
  756.      *    DO 11 L = 1,Loop
  757.      *        X(1)= Y(1)
  758.      *    DO 11 k = 2,n
  759.      * 11     X(k)= X(k-1) + Y(k)
  760.      */
  761.  
  762.     for ( l=1 ; l<=loop ; l++ ) {
  763.         x[0] = y[0];
  764.         for ( k=1 ; k<n ; k++ ) {
  765.             x[k] = x[k-1] + y[k];
  766.         }
  767.     }
  768.     argument = 11;
  769.     TEST( &argument );
  770.  
  771.     /*
  772.      *******************************************************************
  773.      *   Kernel 12 -- first difference
  774.      *******************************************************************
  775.      *    DO 12 L = 1,Loop
  776.      *    DO 12 k = 1,n
  777.      * 12     X(k)= Y(k+1) - Y(k)
  778.      */
  779.  
  780.     for ( l=1 ; l<=loop ; l++ ) {
  781.         for ( k=0 ; k<n ; k++ ) {
  782.             x[k] = y[k+1] - y[k];
  783.         }
  784.     }
  785.     argument = 12;
  786.     TEST( &argument );
  787.  
  788.     /*
  789.      *******************************************************************
  790.      *   Kernel 13 -- 2-D PIC (Particle In Cell)
  791.      *******************************************************************
  792.      *    DO  13     L= 1,Loop
  793.      *    DO  13    ip= 1,n
  794.      *              i1= P(1,ip)
  795.      *              j1= P(2,ip)
  796.      *              i1=        1 + MOD2N(i1,64)
  797.      *              j1=        1 + MOD2N(j1,64)
  798.      *         P(3,ip)= P(3,ip)  + B(i1,j1)
  799.      *         P(4,ip)= P(4,ip)  + C(i1,j1)
  800.      *         P(1,ip)= P(1,ip)  + P(3,ip)
  801.      *         P(2,ip)= P(2,ip)  + P(4,ip)
  802.      *              i2= P(1,ip)
  803.      *              j2= P(2,ip)
  804.      *              i2=            MOD2N(i2,64)
  805.      *              j2=            MOD2N(j2,64)
  806.      *         P(1,ip)= P(1,ip)  + Y(i2+32)
  807.      *         P(2,ip)= P(2,ip)  + Z(j2+32)
  808.      *              i2= i2       + E(i2+32)
  809.      *              j2= j2       + F(j2+32)
  810.      *        H(i2,j2)= H(i2,j2) + 1.0
  811.      * 13 CONTINUE
  812.      */
  813.  
  814.     for ( l=1 ; l<=loop ; l++ ) {
  815.         for ( ip=0 ; ip<n ; ip++ ) {
  816.             i1 = p[ip][0];
  817.             j1 = p[ip][1];
  818.             i1 &= 64-1;
  819.             j1 &= 64-1;
  820.             p[ip][2] += b[j1][i1];
  821.             p[ip][3] += c[j1][i1];
  822.             p[ip][0] += p[ip][2];
  823.             p[ip][1] += p[ip][3];
  824.             i2 = p[ip][0];
  825.             j2 = p[ip][1];
  826.             i2 = ( i2 & 64-1 ) - 1 ;
  827.             j2 = ( j2 & 64-1 ) - 1 ;
  828.             p[ip][0] += y[i2+32];
  829.             p[ip][1] += z[j2+32];
  830.             i2 += e[i2+32];
  831.             j2 += f[j2+32];
  832.             h[j2][i2] += 1.0;
  833.         }
  834.     }
  835.     argument = 13;
  836.     TEST( &argument );
  837.  
  838.     /*
  839.      *******************************************************************
  840.      *   Kernel 14 -- 1-D PIC (Particle In Cell)
  841.      *******************************************************************
  842.      *    DO   14   L= 1,Loop
  843.      *    DO   141  k= 1,n
  844.      *          VX(k)= 0.0
  845.      *          XX(k)= 0.0
  846.      *          IX(k)= INT(  GRD(k))
  847.      *          XI(k)= REAL( IX(k))
  848.      *         EX1(k)= EX   ( IX(k))
  849.      *        DEX1(k)= DEX  ( IX(k))
  850.      *41  CONTINUE
  851.      *    DO   142  k= 1,n
  852.      *          VX(k)= VX(k) + EX1(k) + (XX(k) - XI(k))*DEX1(k)
  853.      *          XX(k)= XX(k) + VX(k)  + FLX
  854.      *          IR(k)= XX(k)
  855.      *          RX(k)= XX(k) - IR(k)
  856.      *          IR(k)= MOD2N(  IR(k),2048) + 1
  857.      *          XX(k)= RX(k) + IR(k)
  858.      *42  CONTINUE
  859.      *    DO  14    k= 1,n
  860.      *    RH(IR(k)  )= RH(IR(k)  ) + 1.0 - RX(k)
  861.      *    RH(IR(k)+1)= RH(IR(k)+1) + RX(k)
  862.      *14  CONTINUE
  863.      */
  864.  
  865.     for ( l=1 ; l<=loop ; l++ ) {
  866.         for ( k=0 ; k<n ; k++ ) {
  867.             vx[k] = 0.0;
  868.             xx[k] = 0.0;
  869.             ix[k] = (long) grd[k];
  870.             xi[k] = (double) ix[k];
  871.             ex1[k] = ex[ ix[k] - 1 ];
  872.             dex1[k] = dex[ ix[k] - 1 ];
  873.         }
  874.         for ( k=0 ; k<n ; k++ ) {
  875.             vx[k] = vx[k] + ex1[k] + ( xx[k] - xi[k] )*dex1[k];
  876.             xx[k] = xx[k] + vx[k]  + flx;
  877.             ir[k] = xx[k];
  878.             rx[k] = xx[k] - ir[k];
  879.             ir[k] = ( ir[k] & 2048-1 ) + 1;
  880.             xx[k] = rx[k] + ir[k];
  881.         }
  882.         for ( k=0 ; k<n ; k++ ) {
  883.             rh[ ir[k]-1 ] += 1.0 - rx[k];
  884.             rh[ ir[k]   ] += rx[k];
  885.         }
  886.     }
  887.     argument = 14;
  888.     TEST( &argument );
  889.  
  890.     /*
  891.      *******************************************************************
  892.      *   Kernel 15 -- Casual Fortran.  Development version
  893.      *******************************************************************
  894.      *      DO 45  L = 1,Loop
  895.      *             NG= 7
  896.      *             NZ= n
  897.      *             AR= 0.053
  898.      *             BR= 0.073
  899.      * 15   DO 45  j = 2,NG
  900.      *      DO 45  k = 2,NZ
  901.      *             IF( j-NG) 31,30,30
  902.      * 30     VY(k,j)= 0.0
  903.      *                 GO TO 45
  904.      * 31          IF( VH(k,j+1) -VH(k,j)) 33,33,32
  905.      * 32           T= AR
  906.      *                 GO TO 34
  907.      * 33           T= BR
  908.      * 34          IF( VF(k,j) -VF(k-1,j)) 35,36,36
  909.      * 35           R= MAX( VH(k-1,j), VH(k-1,j+1))
  910.      *              S= VF(k-1,j)
  911.      *                 GO TO 37
  912.      * 36           R= MAX( VH(k,j),   VH(k,j+1))
  913.      *              S= VF(k,j)
  914.      * 37     VY(k,j)= SQRT( VG(k,j)**2 +R*R)*T/S
  915.      * 38          IF( k-NZ) 40,39,39
  916.      * 39     VS(k,j)= 0.
  917.      *                 GO TO 45
  918.      * 40          IF( VF(k,j) -VF(k,j-1)) 41,42,42
  919.      * 41           R= MAX( VG(k,j-1), VG(k+1,j-1))
  920.      *              S= VF(k,j-1)
  921.      *              T= BR
  922.      *                 GO TO 43
  923.      * 42           R= MAX( VG(k,j),   VG(k+1,j))
  924.      *              S= VF(k,j)
  925.      *              T= AR
  926.      * 43     VS(k,j)= SQRT( VH(k,j)**2 +R*R)*T/S
  927.      * 45    CONTINUE
  928.      */
  929.  
  930. #pragma intrinsic sqrt
  931.  
  932.     for ( l=1 ; l<=loop ; l++ ) {
  933.         ng = 7;
  934.         nz = n;
  935.         ar = 0.053;
  936.         br = 0.073;
  937.         for ( j=1 ; j<ng ; j++ ) {
  938.             for ( k=1 ; k<nz ; k++ ) {
  939.                 if ( (j+1) >= ng ) {
  940.                     vy[j][k] = 0.0;
  941.                     continue;
  942.                 }
  943.                 if ( vh[j+1][k] > vh[j][k] ) {
  944.                     t = ar;
  945.                 }
  946.                 else {
  947.                     t = br;
  948.                 }
  949.                 if ( vf[j][k] < vf[j][k-1] ) {
  950.                     if ( vh[j][k-1] > vh[j+1][k-1] )
  951.                         r = vh[j][k-1];
  952.                     else
  953.                         r = vh[j+1][k-1];
  954.                     s = vf[j][k-1];
  955.                 }
  956.                 else {
  957.                     if ( vh[j][k] > vh[j+1][k] )
  958.                         r = vh[j][k];
  959.                     else
  960.                         r = vh[j+1][k];
  961.                     s = vf[j][k];
  962.                 }
  963.                 vy[j][k] = sqrt( vg[j][k]*vg[j][k] + r*r )* t/s;
  964.                 if ( (k+1) >= nz ) {
  965.                     vs[j][k] = 0.0;
  966.                     continue;
  967.                 }
  968.                 if ( vf[j][k] < vf[j-1][k] ) {
  969.                     if ( vg[j-1][k] > vg[j-1][k+1] )
  970.                         r = vg[j-1][k];
  971.                     else
  972.                         r = vg[j-1][k+1];
  973.                     s = vf[j-1][k];
  974.                     t = br;
  975.                 }
  976.                 else {
  977.                     if ( vg[j][k] > vg[j][k+1] )
  978.                         r = vg[j][k];
  979.                     else
  980.                         r = vg[j][k+1];
  981.                     s = vf[j][k];
  982.                     t = ar;
  983.                 }
  984.                 vs[j][k] = sqrt( vh[j][k]*vh[j][k] + r*r )* t / s;
  985.             }
  986.         }
  987.     }
  988.     argument = 15;
  989.     TEST( &argument );
  990.  
  991.     /*
  992.      *******************************************************************
  993.      *   Kernel 16 -- Monte Carlo search loop
  994.      *******************************************************************
  995.      *          II= n/3
  996.      *          LB= II+II
  997.      *          k2= 0
  998.      *          k3= 0
  999.      *    DO 485 L= 1,Loop
  1000.      *           m= 1
  1001.      *405       i1= m
  1002.      *410       j2= (n+n)*(m-1)+1
  1003.      *    DO 470 k= 1,n
  1004.      *          k2= k2+1
  1005.      *          j4= j2+k+k
  1006.      *          j5= ZONE(j4)
  1007.      *          IF( j5-n      ) 420,475,450
  1008.      *415       IF( j5-n+II   ) 430,425,425
  1009.      *420       IF( j5-n+LB   ) 435,415,415
  1010.      *425       IF( PLAN(j5)-R) 445,480,440
  1011.      *430       IF( PLAN(j5)-S) 445,480,440
  1012.      *435       IF( PLAN(j5)-T) 445,480,440
  1013.      *440       IF( ZONE(j4-1)) 455,485,470
  1014.      *445       IF( ZONE(j4-1)) 470,485,455
  1015.      *450       k3= k3+1
  1016.      *          IF( D(j5)-(D(j5-1)*(T-D(j5-2))**2+(S-D(j5-3))**2
  1017.      *   .                        +(R-D(j5-4))**2)) 445,480,440
  1018.      *455        m= m+1
  1019.      *          IF( m-ZONE(1) ) 465,465,460
  1020.      *460        m= 1
  1021.      *465       IF( i1-m) 410,480,410
  1022.      *470 CONTINUE
  1023.      *475 CONTINUE
  1024.      *480 CONTINUE
  1025.      *485 CONTINUE
  1026.      */
  1027.  
  1028.     ii = n / 3;
  1029.     lb = ii + ii;
  1030.     k3 = k2 = 0;
  1031.     for ( l=1 ; l<=loop ; l++ ) {
  1032.         i1 = m = 1;
  1033.         label410:
  1034.         j2 = ( n + n )*( m - 1 ) + 1;
  1035.         for ( k=1 ; k<=n ; k++ ) {
  1036.             k2++;
  1037.             j4 = j2 + k + k;
  1038.             j5 = zone[j4-1];
  1039.             if ( j5 < n ) {
  1040.                 if ( j5+lb < n ) {              /* 420 */
  1041.                     tmp = plan[j5-1] - t;       /* 435 */
  1042.                 } else {
  1043.                     if ( j5+ii < n ) {          /* 415 */
  1044.                         tmp = plan[j5-1] - s;   /* 430 */
  1045.                     } else {
  1046.                         tmp = plan[j5-1] - r;   /* 425 */
  1047.                     }
  1048.                 }
  1049.             } else if( j5 == n ) {
  1050.                 break;                          /* 475 */
  1051.             } else {
  1052.                 k3++;                           /* 450 */
  1053.                 tmp=(d[j5-1]-(d[j5-2]*(t-d[j5-3])*(t-d[j5-3])+(s-d[j5-4])*
  1054.                               (s-d[j5-4])+(r-d[j5-5])*(r-d[j5-5])));
  1055.             }
  1056.             if ( tmp < 0.0 ) {
  1057.                 if ( zone[j4-2] < 0 )           /* 445 */
  1058.                     continue;                   /* 470 */
  1059.                 else if ( !zone[j4-2] )
  1060.                     break;                      /* 480 */
  1061.             } else if ( tmp ) {
  1062.                 if ( zone[j4-2] > 0 )           /* 440 */
  1063.                     continue;                   /* 470 */
  1064.                 else if ( !zone[j4-2] )
  1065.                     break;                      /* 480 */
  1066.             } else break;                       /* 485 */
  1067.             m++;                                /* 455 */
  1068.             if ( m > zone[0] )
  1069.                 m = 1;                          /* 460 */
  1070.             if ( i1-m )                         /* 465 */
  1071.                 goto label410;
  1072.             else
  1073.                 break;
  1074.         }
  1075.     }
  1076.     argument = 16;
  1077.     TEST( &argument );
  1078.  
  1079.     /*
  1080.      *******************************************************************
  1081.      *   Kernel 17 -- implicit, conditional computation
  1082.      *******************************************************************
  1083.      *          DO 62 L= 1,Loop
  1084.      *                i= n
  1085.      *                j= 1
  1086.      *              INK= -1
  1087.      *            SCALE= 5./3.
  1088.      *              XNM= 1./3.
  1089.      *               E6= 1.03/3.07
  1090.      *                   GO TO 61
  1091.      *60             E6= XNM*VSP(i)+VSTP(i)
  1092.      *          VXNE(i)= E6
  1093.      *              XNM= E6
  1094.      *           VE3(i)= E6
  1095.      *                i= i+INK
  1096.      *               IF( i.EQ.j) GO TO  62
  1097.      *61             E3= XNM*VLR(i) +VLIN(i)
  1098.      *             XNEI= VXNE(i)
  1099.      *          VXND(i)= E6
  1100.      *              XNC= SCALE*E3
  1101.      *               IF( XNM .GT.XNC) GO TO  60
  1102.      *               IF( XNEI.GT.XNC) GO TO  60
  1103.      *           VE3(i)= E3
  1104.      *               E6= E3+E3-XNM
  1105.      *          VXNE(i)= E3+E3-XNEI
  1106.      *              XNM= E6
  1107.      *                i= i+INK
  1108.      *               IF( i.NE.j) GO TO 61
  1109.      * 62 CONTINUE
  1110.      */
  1111.  
  1112.     for ( l=1 ; l<=loop ; l++ ) {
  1113.         i = n-1;
  1114.         j = 0;
  1115.         ink = -1;
  1116.         scale = 5.0 / 3.0;
  1117.         xnm = 1.0 / 3.0;
  1118.         e6 = 1.03 / 3.07;
  1119.         goto l61;
  1120. l60:    e6 = xnm*vsp[i] + vstp[i];
  1121.         vxne[i] = e6;
  1122.         xnm = e6;
  1123.         ve3[i] = e6;
  1124.         i += ink;
  1125.         if ( i==j ) goto l62;
  1126. l61:    e3 = xnm*vlr[i] + vlin[i];
  1127.         xnei = vxne[i];
  1128.         vxnd[i] = e6;
  1129.         xnc = scale*e3;
  1130.         if ( xnm > xnc ) goto l60;
  1131.         if ( xnei > xnc ) goto l60;
  1132.         ve3[i] = e3;
  1133.         e6 = e3 + e3 - xnm;
  1134.         vxne[i] = e3 + e3 - xnei;
  1135.         xnm = e6;
  1136.         i += ink;
  1137.         if ( i != j ) goto l61;
  1138. l62:;
  1139.     }
  1140.  
  1141.     argument = 17;
  1142.     TEST( &argument );
  1143.  
  1144.     /*
  1145.      *******************************************************************
  1146.      *   Kernel 18 - 2-D explicit hydrodynamics fragment
  1147.      *******************************************************************
  1148.      *       DO 75  L= 1,Loop
  1149.      *              T= 0.0037
  1150.      *              S= 0.0041
  1151.      *             KN= 6
  1152.      *             JN= n
  1153.      *       DO 70  k= 2,KN
  1154.      *       DO 70  j= 2,JN
  1155.      *        ZA(j,k)= (ZP(j-1,k+1)+ZQ(j-1,k+1)-ZP(j-1,k)-ZQ(j-1,k))
  1156.      *   .            *(ZR(j,k)+ZR(j-1,k))/(ZM(j-1,k)+ZM(j-1,k+1))
  1157.      *        ZB(j,k)= (ZP(j-1,k)+ZQ(j-1,k)-ZP(j,k)-ZQ(j,k))
  1158.      *   .            *(ZR(j,k)+ZR(j,k-1))/(ZM(j,k)+ZM(j-1,k))
  1159.      * 70    CONTINUE
  1160.      *       DO 72  k= 2,KN
  1161.      *       DO 72  j= 2,JN
  1162.      *        ZU(j,k)= ZU(j,k)+S*(ZA(j,k)*(ZZ(j,k)-ZZ(j+1,k))
  1163.      *   .                    -ZA(j-1,k) *(ZZ(j,k)-ZZ(j-1,k))
  1164.      *   .                    -ZB(j,k)   *(ZZ(j,k)-ZZ(j,k-1))
  1165.      *   .                    +ZB(j,k+1) *(ZZ(j,k)-ZZ(j,k+1)))
  1166.      *        ZV(j,k)= ZV(j,k)+S*(ZA(j,k)*(ZR(j,k)-ZR(j+1,k))
  1167.      *   .                    -ZA(j-1,k) *(ZR(j,k)-ZR(j-1,k))
  1168.      *   .                    -ZB(j,k)   *(ZR(j,k)-ZR(j,k-1))
  1169.      *   .                    +ZB(j,k+1) *(ZR(j,k)-ZR(j,k+1)))
  1170.      * 72    CONTINUE
  1171.      *       DO 75  k= 2,KN
  1172.      *       DO 75  j= 2,JN
  1173.      *        ZR(j,k)= ZR(j,k)+T*ZU(j,k)
  1174.      *        ZZ(j,k)= ZZ(j,k)+T*ZV(j,k)
  1175.      * 75    CONTINUE
  1176.      */
  1177.  
  1178.     for ( l=1 ; l<=loop ; l++ ) {
  1179.         t = 0.0037;
  1180.         s = 0.0041;
  1181.         kn = 6;
  1182.         jn = n;
  1183.         for ( k=1 ; k<kn ; k++ ) {
  1184. #pragma nohazard
  1185.           for ( j=1 ; j<jn ; j++ ) {
  1186.               za[k][j] = ( zp[k+1][j-1] +zq[k+1][j-1] -zp[k][j-1] -zq[k][j-1] )*
  1187.                          ( zr[k][j] +zr[k][j-1] ) / ( zm[k][j-1] +zm[k+1][j-1]);
  1188.               zb[k][j] = ( zp[k][j-1] +zq[k][j-1] -zp[k][j] -zq[k][j] ) *
  1189.                          ( zr[k][j] +zr[k-1][j] ) / ( zm[k][j] +zm[k][j-1]);
  1190.           }
  1191.         }
  1192.         for ( k=1 ; k<kn ; k++ ) {
  1193. #pragma nohazard
  1194.             for ( j=1 ; j<jn ; j++ ) {
  1195.                 zu[k][j] += s*( za[k][j]   *( zz[k][j] - zz[k][j+1] ) -
  1196.                                 za[k][j-1] *( zz[k][j] - zz[k][j-1] ) -
  1197.                                 zb[k][j]   *( zz[k][j] - zz[k-1][j] ) +
  1198.                                 zb[k+1][j] *( zz[k][j] - zz[k+1][j] ) );
  1199.                 zv[k][j] += s*( za[k][j]   *( zr[k][j] - zr[k][j+1] ) -
  1200.                                 za[k][j-1] *( zr[k][j] - zr[k][j-1] ) -
  1201.                                 zb[k][j]   *( zr[k][j] - zr[k-1][j] ) +
  1202.                                 zb[k+1][j] *( zr[k][j] - zr[k+1][j] ) );
  1203.             }
  1204.         }
  1205.         for ( k=1 ; k<kn ; k++ ) {
  1206. #pragma nohazard
  1207.             for ( j=1 ; j<jn ; j++ ) {
  1208.                 zr[k][j] = zr[k][j] + t*zu[k][j];
  1209.                 zz[k][j] = zz[k][j] + t*zv[k][j];
  1210.             }
  1211.         }
  1212.     }
  1213.     argument = 18;
  1214.     TEST( &argument );
  1215.  
  1216.     /*
  1217.      *******************************************************************
  1218.      *   Kernel 19 -- general linear recurrence equations
  1219.      *******************************************************************
  1220.      *               KB5I= 0
  1221.      *           DO 194 L= 1,Loop
  1222.      *           DO 191 k= 1,n
  1223.      *         B5(k+KB5I)= SA(k) +STB5*SB(k)
  1224.      *               STB5= B5(k+KB5I) -STB5
  1225.      *191        CONTINUE
  1226.      *192        DO 193 i= 1,n
  1227.      *                  k= n-i+1
  1228.      *         B5(k+KB5I)= SA(k) +STB5*SB(k)
  1229.      *               STB5= B5(k+KB5I) -STB5
  1230.      *193        CONTINUE
  1231.      *194 CONTINUE
  1232.      */
  1233.  
  1234.     kb5i = 0;
  1235.     for ( l=1 ; l<=loop ; l++ ) {
  1236.         for ( k=0 ; k<n ; k++ ) {
  1237.             b5[k+kb5i] = sa[k] + stb5*sb[k];
  1238.             stb5 = b5[k+kb5i] - stb5;
  1239.         }
  1240.         for ( i=1 ; i<=n ; i++ ) {
  1241.             k = n - i ;
  1242.             b5[k+kb5i] = sa[k] + stb5*sb[k];
  1243.             stb5 = b5[k+kb5i] - stb5;
  1244.         }
  1245.     }
  1246.     argument = 19;
  1247.     TEST( &argument );
  1248.  
  1249.     /*
  1250.      *******************************************************************
  1251.      *   Kernel 20 -- Discrete ordinates transport, conditional recurrence on xx
  1252.      *******************************************************************
  1253.      *    DO 20 L= 1,Loop
  1254.      *    DO 20 k= 1,n
  1255.      *         DI= Y(k)-G(k)/( XX(k)+DK)
  1256.      *         DN= 0.2
  1257.      *         IF( DI.NE.0.0) DN= MAX( S,MIN( Z(k)/DI, T))
  1258.      *       X(k)= ((W(k)+V(k)*DN)* XX(k)+U(k))/(VX(k)+V(k)*DN)
  1259.      *    XX(k+1)= (X(k)- XX(k))*DN+ XX(k)
  1260.      * 20 CONTINUE
  1261.      */
  1262.  
  1263.     for ( l=1 ; l<=loop ; l++ ) {
  1264.         for ( k=0 ; k<n ; k++ ) {
  1265.             di = y[k] - g[k] / ( xx[k] + dk );
  1266.             dn = 0.2;
  1267.             if ( di ) {
  1268.                 dn = z[k]/di ;
  1269.                 if ( t < dn ) dn = t;
  1270.                 if ( s > dn ) dn = s;
  1271.             }
  1272.             x[k] = ( ( w[k] + v[k]*dn )* xx[k] + u[k] ) / ( vx[k] + v[k]*dn );
  1273.             xx[k+1] = ( x[k] - xx[k] )* dn + xx[k];
  1274.         }
  1275.     }
  1276.     argument = 20;
  1277.     TEST( &argument );
  1278.  
  1279.     /*
  1280.      *******************************************************************
  1281.      *   Kernel 21 -- matrix*matrix product
  1282.      *******************************************************************
  1283.      *    DO 21 L= 1,Loop
  1284.      *    DO 21 k= 1,25
  1285.      *    DO 21 i= 1,25
  1286.      *    DO 21 j= 1,n
  1287.      *    PX(i,j)= PX(i,j) +VY(i,k) * CX(k,j)
  1288.      * 21 CONTINUE
  1289.      */
  1290.  
  1291.     for ( l=1 ; l<=loop ; l++ ) {
  1292.         for ( k=0 ; k<25 ; k++ ) {
  1293.             for ( i=0 ; i<25 ; i++ ) {
  1294. #pragma nohazard
  1295.                 for ( j=0 ; j<n ; j++ ) {
  1296.                     px[j][i] += vy[k][i] * cx[j][k];
  1297.                 }
  1298.             }
  1299.         }
  1300.     }
  1301.     argument = 21;
  1302.     TEST( &argument );
  1303.  
  1304.     /*
  1305.      *******************************************************************
  1306.      *   Kernel 22 -- Planckian distribution
  1307.      *******************************************************************
  1308.      *     EXPMAX= 20.0
  1309.      *       U(n)= 0.99*EXPMAX*V(n)
  1310.      *    DO 22 L= 1,Loop
  1311.      *    DO 22 k= 1,n
  1312.      *                                          Y(k)= U(k)/V(k)
  1313.      *       W(k)= X(k)/( EXP( Y(k)) -1.0)
  1314.      * 22 CONTINUE
  1315.      */
  1316.  
  1317. #pragma intrinsic exp
  1318.     expmax = 20.0;
  1319.     u[n-1] = 0.99*expmax*v[n-1];
  1320.     for ( l=1 ; l<=loop ; l++ ) {
  1321.         for ( k=0 ; k<n ; k++ ) {
  1322.             y[k] = u[k] / v[k];
  1323.             w[k] = x[k] / ( exp( y[k] ) -1.0 );
  1324.         }
  1325.     }
  1326.     argument = 22;
  1327.     TEST( &argument );
  1328.  
  1329.     /*
  1330.      *******************************************************************
  1331.      *   Kernel 23 -- 2-D implicit hydrodynamics fragment
  1332.      *******************************************************************
  1333.      *    DO 23  L= 1,Loop
  1334.      *    DO 23  j= 2,6
  1335.      *    DO 23  k= 2,n
  1336.      *          QA= ZA(k,j+1)*ZR(k,j) +ZA(k,j-1)*ZB(k,j) +
  1337.      *   .          ZA(k+1,j)*ZU(k,j) +ZA(k-1,j)*ZV(k,j) +ZZ(k,j)
  1338.      * 23  ZA(k,j)= ZA(k,j) +.175*(QA -ZA(k,j))
  1339.      */
  1340.  
  1341.     for ( l=1 ; l<=loop ; l++ ) {
  1342.         for ( j=1 ; j<6 ; j++ ) {
  1343.             for ( k=1 ; k<n ; k++ ) {
  1344.                 qa = za[j+1][k]*zr[j][k] + za[j-1][k]*zb[j][k] +
  1345.                      za[j][k+1]*zu[j][k] + za[j][k-1]*zv[j][k] + zz[j][k];
  1346.                 za[j][k] += 0.175*( qa - za[j][k] );
  1347.             }
  1348.         }
  1349.     }
  1350.     argument = 23;
  1351.     TEST( &argument );
  1352.  
  1353.     /*
  1354.      *******************************************************************
  1355.      *   Kernel 24 -- find location of first minimum in array
  1356.      *******************************************************************
  1357.      *     X( n/2)= -1.0E+10
  1358.      *    DO 24  L= 1,Loop
  1359.      *           m= 1
  1360.      *    DO 24  k= 2,n
  1361.      *          IF( X(k).LT.X(m))  m= k
  1362.      * 24 CONTINUE
  1363.      */
  1364.  
  1365.     x[n/2] = -1.0e+10;
  1366.     for ( l=1 ; l<=loop ; l++ ) {
  1367.         m = 0;
  1368.         for ( k=1 ; k<n ; k++ ) {
  1369.             if ( x[k] < x[m] ) m = k;
  1370.         }
  1371.     }
  1372.     argument = 24;
  1373.     TEST( &argument );
  1374.  
  1375.     /*
  1376.      *   Epilogue
  1377.      */
  1378.  
  1379.     if ( jr < 1 )
  1380.         jr = 1;
  1381.     if ( jr > mr )
  1382.         jr = 8;
  1383.     if ( il < 1 )
  1384.         il = 1;
  1385.     if ( il > ml )
  1386.         il = 3;
  1387.  
  1388.     for ( k=0 ; k<mk ; k++ ) {
  1389.         times[k][il-1][jr-1] = time[k];
  1390.         terrs[k][il-1][jr-1] = terr1[k];
  1391.         npfs[k][il-1][jr-1] = npfs1[k];
  1392.         csums[k][il-1][jr-1] = csum[k];
  1393.         dos[k][il-1][jr-1] =  runs[k];
  1394.         fopn[k][il-1][jr-1] =  flopn[k];
  1395.     }
  1396.  
  1397.     return;
  1398. }
  1399.  
  1400. /* End of File */
  1401.  
  1402.  
  1403. +++++++++++++++++++++++++++++++++++++++++++++++++
  1404. +++ ED (Edward Quillen)   (517) 336-1293      +++
  1405. +++ Vet Teaching Hosp.    System Admin.       +++
  1406. +++ Email: quillen@cps.msu.edu                +++
  1407. +++++++++++++++++++++++++++++++++++++++++++++++++
  1408.