📄 livermore.c
字号:
#pragma intrinsic sqrt
for ( l=1 ; l<=loop ; l++ ) {
ng = 7;
nz = n;
ar = 0.053;
br = 0.073;
for ( j=1 ; j<ng ; j++ ) {
for ( k=1 ; k<nz ; k++ ) {
if ( (j+1) >= ng ) {
vy[j][k] = 0.0;
continue;
}
if ( vh[j+1][k] > vh[j][k] ) {
t = ar;
}
else {
t = br;
}
if ( vf[j][k] < vf[j][k-1] ) {
if ( vh[j][k-1] > vh[j+1][k-1] )
r = vh[j][k-1];
else
r = vh[j+1][k-1];
s = vf[j][k-1];
}
else {
if ( vh[j][k] > vh[j+1][k] )
r = vh[j][k];
else
r = vh[j+1][k];
s = vf[j][k];
}
vy[j][k] = sqrt( vg[j][k]*vg[j][k] + r*r )* t/s;
if ( (k+1) >= nz ) {
vs[j][k] = 0.0;
continue;
}
if ( vf[j][k] < vf[j-1][k] ) {
if ( vg[j-1][k] > vg[j-1][k+1] )
r = vg[j-1][k];
else
r = vg[j-1][k+1];
s = vf[j-1][k];
t = br;
}
else {
if ( vg[j][k] > vg[j][k+1] )
r = vg[j][k];
else
r = vg[j][k+1];
s = vf[j][k];
t = ar;
}
vs[j][k] = sqrt( vh[j][k]*vh[j][k] + r*r )* t / s;
}
}
}
argument = 15;
TEST( &argument );
/*
*******************************************************************
* Kernel 16 -- Monte Carlo search loop
*******************************************************************
* II= n/3
* LB= II+II
* k2= 0
* k3= 0
* DO 485 L= 1,Loop
* m= 1
*405 i1= m
*410 j2= (n+n)*(m-1)+1
* DO 470 k= 1,n
* k2= k2+1
* j4= j2+k+k
* j5= ZONE(j4)
* IF( j5-n ) 420,475,450
*415 IF( j5-n+II ) 430,425,425
*420 IF( j5-n+LB ) 435,415,415
*425 IF( PLAN(j5)-R) 445,480,440
*430 IF( PLAN(j5)-S) 445,480,440
*435 IF( PLAN(j5)-T) 445,480,440
*440 IF( ZONE(j4-1)) 455,485,470
*445 IF( ZONE(j4-1)) 470,485,455
*450 k3= k3+1
* IF( D(j5)-(D(j5-1)*(T-D(j5-2))**2+(S-D(j5-3))**2
* . +(R-D(j5-4))**2)) 445,480,440
*455 m= m+1
* IF( m-ZONE(1) ) 465,465,460
*460 m= 1
*465 IF( i1-m) 410,480,410
*470 CONTINUE
*475 CONTINUE
*480 CONTINUE
*485 CONTINUE
*/
ii = n / 3;
lb = ii + ii;
k3 = k2 = 0;
for ( l=1 ; l<=loop ; l++ ) {
i1 = m = 1;
label410:
j2 = ( n + n )*( m - 1 ) + 1;
for ( k=1 ; k<=n ; k++ ) {
k2++;
j4 = j2 + k + k;
j5 = zone[j4-1];
if ( j5 < n ) {
if ( j5+lb < n ) { /* 420 */
tmp = plan[j5-1] - t; /* 435 */
} else {
if ( j5+ii < n ) { /* 415 */
tmp = plan[j5-1] - s; /* 430 */
} else {
tmp = plan[j5-1] - r; /* 425 */
}
}
} else if( j5 == n ) {
break; /* 475 */
} else {
k3++; /* 450 */
tmp=(d[j5-1]-(d[j5-2]*(t-d[j5-3])*(t-d[j5-3])+(s-d[j5-4])*
(s-d[j5-4])+(r-d[j5-5])*(r-d[j5-5])));
}
if ( tmp < 0.0 ) {
if ( zone[j4-2] < 0 ) /* 445 */
continue; /* 470 */
else if ( !zone[j4-2] )
break; /* 480 */
} else if ( tmp ) {
if ( zone[j4-2] > 0 ) /* 440 */
continue; /* 470 */
else if ( !zone[j4-2] )
break; /* 480 */
} else break; /* 485 */
m++; /* 455 */
if ( m > zone[0] )
m = 1; /* 460 */
if ( i1-m ) /* 465 */
goto label410;
else
break;
}
}
argument = 16;
TEST( &argument );
/*
*******************************************************************
* Kernel 17 -- implicit, conditional computation
*******************************************************************
* DO 62 L= 1,Loop
* i= n
* j= 1
* INK= -1
* SCALE= 5./3.
* XNM= 1./3.
* E6= 1.03/3.07
* GO TO 61
*60 E6= XNM*VSP(i)+VSTP(i)
* VXNE(i)= E6
* XNM= E6
* VE3(i)= E6
* i= i+INK
* IF( i.EQ.j) GO TO 62
*61 E3= XNM*VLR(i) +VLIN(i)
* XNEI= VXNE(i)
* VXND(i)= E6
* XNC= SCALE*E3
* IF( XNM .GT.XNC) GO TO 60
* IF( XNEI.GT.XNC) GO TO 60
* VE3(i)= E3
* E6= E3+E3-XNM
* VXNE(i)= E3+E3-XNEI
* XNM= E6
* i= i+INK
* IF( i.NE.j) GO TO 61
* 62 CONTINUE
*/
for ( l=1 ; l<=loop ; l++ ) {
i = n-1;
j = 0;
ink = -1;
scale = 5.0 / 3.0;
xnm = 1.0 / 3.0;
e6 = 1.03 / 3.07;
goto l61;
l60: e6 = xnm*vsp[i] + vstp[i];
vxne[i] = e6;
xnm = e6;
ve3[i] = e6;
i += ink;
if ( i==j ) goto l62;
l61: e3 = xnm*vlr[i] + vlin[i];
xnei = vxne[i];
vxnd[i] = e6;
xnc = scale*e3;
if ( xnm > xnc ) goto l60;
if ( xnei > xnc ) goto l60;
ve3[i] = e3;
e6 = e3 + e3 - xnm;
vxne[i] = e3 + e3 - xnei;
xnm = e6;
i += ink;
if ( i != j ) goto l61;
l62:;
}
argument = 17;
TEST( &argument );
/*
*******************************************************************
* Kernel 18 - 2-D explicit hydrodynamics fragment
*******************************************************************
* DO 75 L= 1,Loop
* T= 0.0037
* S= 0.0041
* KN= 6
* JN= n
* DO 70 k= 2,KN
* DO 70 j= 2,JN
* ZA(j,k)= (ZP(j-1,k+1)+ZQ(j-1,k+1)-ZP(j-1,k)-ZQ(j-1,k))
* . *(ZR(j,k)+ZR(j-1,k))/(ZM(j-1,k)+ZM(j-1,k+1))
* ZB(j,k)= (ZP(j-1,k)+ZQ(j-1,k)-ZP(j,k)-ZQ(j,k))
* . *(ZR(j,k)+ZR(j,k-1))/(ZM(j,k)+ZM(j-1,k))
* 70 CONTINUE
* DO 72 k= 2,KN
* DO 72 j= 2,JN
* ZU(j,k)= ZU(j,k)+S*(ZA(j,k)*(ZZ(j,k)-ZZ(j+1,k))
* . -ZA(j-1,k) *(ZZ(j,k)-ZZ(j-1,k))
* . -ZB(j,k) *(ZZ(j,k)-ZZ(j,k-1))
* . +ZB(j,k+1) *(ZZ(j,k)-ZZ(j,k+1)))
* ZV(j,k)= ZV(j,k)+S*(ZA(j,k)*(ZR(j,k)-ZR(j+1,k))
* . -ZA(j-1,k) *(ZR(j,k)-ZR(j-1,k))
* . -ZB(j,k) *(ZR(j,k)-ZR(j,k-1))
* . +ZB(j,k+1) *(ZR(j,k)-ZR(j,k+1)))
* 72 CONTINUE
* DO 75 k= 2,KN
* DO 75 j= 2,JN
* ZR(j,k)= ZR(j,k)+T*ZU(j,k)
* ZZ(j,k)= ZZ(j,k)+T*ZV(j,k)
* 75 CONTINUE
*/
for ( l=1 ; l<=loop ; l++ ) {
t = 0.0037;
s = 0.0041;
kn = 6;
jn = n;
for ( k=1 ; k<kn ; k++ ) {
#pragma nohazard
for ( j=1 ; j<jn ; j++ ) {
za[k][j] = ( zp[k+1][j-1] +zq[k+1][j-1] -zp[k][j-1] -zq[k][j-1] )*
( zr[k][j] +zr[k][j-1] ) / ( zm[k][j-1] +zm[k+1][j-1]);
zb[k][j] = ( zp[k][j-1] +zq[k][j-1] -zp[k][j] -zq[k][j] ) *
( zr[k][j] +zr[k-1][j] ) / ( zm[k][j] +zm[k][j-1]);
}
}
for ( k=1 ; k<kn ; k++ ) {
#pragma nohazard
for ( j=1 ; j<jn ; j++ ) {
zu[k][j] += s*( za[k][j] *( zz[k][j] - zz[k][j+1] ) -
za[k][j-1] *( zz[k][j] - zz[k][j-1] ) -
zb[k][j] *( zz[k][j] - zz[k-1][j] ) +
zb[k+1][j] *( zz[k][j] - zz[k+1][j] ) );
zv[k][j] += s*( za[k][j] *( zr[k][j] - zr[k][j+1] ) -
za[k][j-1] *( zr[k][j] - zr[k][j-1] ) -
zb[k][j] *( zr[k][j] - zr[k-1][j] ) +
zb[k+1][j] *( zr[k][j] - zr[k+1][j] ) );
}
}
for ( k=1 ; k<kn ; k++ ) {
#pragma nohazard
for ( j=1 ; j<jn ; j++ ) {
zr[k][j] = zr[k][j] + t*zu[k][j];
zz[k][j] = zz[k][j] + t*zv[k][j];
}
}
}
argument = 18;
TEST( &argument );
/*
*******************************************************************
* Kernel 19 -- general linear recurrence equations
*******************************************************************
* KB5I= 0
* DO 194 L= 1,Loop
* DO 191 k= 1,n
* B5(k+KB5I)= SA(k) +STB5*SB(k)
* STB5= B5(k+KB5I) -STB5
*191 CONTINUE
*192 DO 193 i= 1,n
* k= n-i+1
* B5(k+KB5I)= SA(k) +STB5*SB(k)
* STB5= B5(k+KB5I) -STB5
*193 CONTINUE
*194 CONTINUE
*/
kb5i = 0;
for ( l=1 ; l<=loop ; l++ ) {
for ( k=0 ; k<n ; k++ ) {
b5[k+kb5i] = sa[k] + stb5*sb[k];
stb5 = b5[k+kb5i] - stb5;
}
for ( i=1 ; i<=n ; i++ ) {
k = n - i ;
b5[k+kb5i] = sa[k] + stb5*sb[k];
stb5 = b5[k+kb5i] - stb5;
}
}
argument = 19;
TEST( &argument );
/*
*******************************************************************
* Kernel 20 -- Discrete ordinates transport, conditional recurrence on xx
*******************************************************************
* DO 20 L= 1,Loop
* DO 20 k= 1,n
* DI= Y(k)-G(k)/( XX(k)+DK)
* DN= 0.2
* IF( DI.NE.0.0) DN= MAX( S,MIN( Z(k)/DI, T))
* X(k)= ((W(k)+V(k)*DN)* XX(k)+U(k))/(VX(k)+V(k)*DN)
* XX(k+1)= (X(k)- XX(k))*DN+ XX(k)
* 20 CONTINUE
*/
for ( l=1 ; l<=loop ; l++ ) {
for ( k=0 ; k<n ; k++ ) {
di = y[k] - g[k] / ( xx[k] + dk );
dn = 0.2;
if ( di ) {
dn = z[k]/di ;
if ( t < dn ) dn = t;
if ( s > dn ) dn = s;
}
x[k] = ( ( w[k] + v[k]*dn )* xx[k] + u[k] ) / ( vx[k] + v[k]*dn );
xx[k+1] = ( x[k] - xx[k] )* dn + xx[k];
}
}
argument = 20;
TEST( &argument );
/*
*******************************************************************
* Kernel 21 -- matrix*matrix product
*******************************************************************
* DO 21 L= 1,Loop
* DO 21 k= 1,25
* DO 21 i= 1,25
* DO 21 j= 1,n
* PX(i,j)= PX(i,j) +VY(i,k) * CX(k,j)
* 21 CONTINUE
*/
for ( l=1 ; l<=loop ; l++ ) {
for ( k=0 ; k<25 ; k++ ) {
for ( i=0 ; i<25 ; i++ ) {
#pragma nohazard
for ( j=0 ; j<n ; j++ ) {
px[j][i] += vy[k][i] * cx[j][k];
}
}
}
}
argument = 21;
TEST( &argument );
/*
*******************************************************************
* Kernel 22 -- Planckian distribution
*******************************************************************
* EXPMAX= 20.0
* U(n)= 0.99*EXPMAX*V(n)
* DO 22 L= 1,Loop
* DO 22 k= 1,n
* Y(k)= U(k)/V(k)
* W(k)= X(k)/( EXP( Y(k)) -1.0)
* 22 CONTINUE
*/
#pragma intrinsic exp
expmax = 20.0;
u[n-1] = 0.99*expmax*v[n-1];
for ( l=1 ; l<=loop ; l++ ) {
for ( k=0 ; k<n ; k++ ) {
y[k] = u[k] / v[k];
w[k] = x[k] / ( exp( y[k] ) -1.0 );
}
}
argument = 22;
TEST( &argument );
/*
*******************************************************************
* Kernel 23 -- 2-D implicit hydrodynamics fragment
*******************************************************************
* DO 23 L= 1,Loop
* DO 23 j= 2,6
* DO 23 k= 2,n
* QA= ZA(k,j+1)*ZR(k,j) +ZA(k,j-1)*ZB(k,j) +
* . ZA(k+1,j)*ZU(k,j) +ZA(k-1,j)*ZV(k,j) +ZZ(k,j)
* 23 ZA(k,j)= ZA(k,j) +.175*(QA -ZA(k,j))
*/
for ( l=1 ; l<=loop ; l++ ) {
for ( j=1 ; j<6 ; j++ ) {
for ( k=1 ; k<n ; k++ ) {
qa = za[j+1][k]*zr[j][k] + za[j-1][k]*zb[j][k] +
za[j][k+1]*zu[j][k] + za[j][k-1]*zv[j][k] + zz[j][k];
za[j][k] += 0.175*( qa - za[j][k] );
}
}
}
argument = 23;
TEST( &argument );
/*
*******************************************************************
* Kernel 24 -- find location of first minimum in array
*******************************************************************
* X( n/2)= -1.0E+10
* DO 24 L= 1,Loop
* m= 1
* DO 24 k= 2,n
* IF( X(k).LT.X(m)) m= k
* 24 CONTINUE
*/
x[n/2] = -1.0e+10;
for ( l=1 ; l<=loop ; l++ ) {
m = 0;
for ( k=1 ; k<n ; k++ ) {
if ( x[k] < x[m] ) m = k;
}
}
argument = 24;
TEST( &argument );
/*
* Epilogue
*/
for ( k=0 ; k<mk ; k++ ) {
times[k][il-1][jr-1] = time[k];
fopn[k][il-1][jr-1] = flopn[k];
terrs[k][il-1][jr-1] = terr1[k];
npfs[k][il-1][jr-1] = npfs1[k];
csums[k][il-1][jr-1] = csum[k];
dos[k][il-1][jr-1] = total[k];
}
sum = 0.0;
for ( k=0 ; k<mk ; k++ ) {
sum += time[k];
}
TK[0] += sum;
som = 0.0;
for ( k=0 ; k<mk ; k++ ) {
som += flopn[k]*total[k];
}
TK[1] += som;
TRACK( &name );
return;
}
/* End of File */
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -