gcc/testsuite/gfortran.dg/pr88148.f90

   1 ! { dg-do compile }
   2 ! { dg-options "-O -fno-tree-fre -fno-tree-sra -ftree-loop-vectorize" }
   3 ! { dg-additional-options "-mavx2" { target x86_64-*-* i?86-*-* } }
   4
   5 module lfk_prec
   6  integer, parameter :: dp=kind(1.d0)
   7 end module lfk_prec
   8
   9 !***********************************************
  10
  11 SUBROUTINE kernel(tk)
  12 !***********************************************************************
  13 !                                                                      *
  14 !            KERNEL     executes 24 samples of Fortran computation     *
  15 !               TK(1) - total cpu time to execute only the 24 kernels. *
  16 !               TK(2) - total Flops executed by the 24 Kernels         *
  17 !***********************************************************************
  18 !                                                                      *
  19 !     L. L. N. L.   F O R T R A N   K E R N E L S:   M F L O P S       *
  20 !                                                                      *
  21 !   These kernels measure  Fortran  numerical  computation rates for a *
  22 !   spectrum of  CPU-limited  computational  structures.  Mathematical *
  23 !   through-put is measured  in  units  of  millions of floating-point *
  24 !   operations executed per Second, called Mega-Flops/Sec.             *
  25 !                                                                      *
  26 !   This program  measures  a realistic  CPU performance range for the *
  27 !   Fortran programming system  on  a  given day.  The CPU performance *
  28 !   rates depend  strongly  on  the maturity of the Fortran compiler's *
  29 !   ability to translate Fortran code into efficient machine code.     *
  30 !   [ The CPU hardware  capability  apart  from  compiler maturity (or *
  31 !   availability), could be measured (or simulated) by programming the *
  32 !   kernels in assembly  or machine code directly.  These measurements *
  33 !   can also  serve  as a framework for tracking the maturation of the *
  34 !   Fortran compiler during system development.]                       *
  35 !                                                                      *
  36 !     Fonzi's Law: There is not now and there never will be a language *
  37 !                  in which it is the least bit difficult to write     *
  38 !                  bad programs.                                       *
  39 !                                                    F.H.MCMAHON  1972 *
  40 !***********************************************************************
  41
  42 !     l1 :=  param-dimension governs the size of most 1-d arrays
  43 !     l2 :=  param-dimension governs the size of most 2-d arrays
  44
  45 !     Loop :=  multiple pass control to execute kernel long enough to ti
  46 !    me.
  47 !     n  :=  DO loop control for each kernel.  Controls are set in subr.
  48 !     SIZES
  49
  50 !     ******************************************************************
  51 use lfk_prec
  52 implicit double precision  (a-h,o-z)
  53 !IBM  IMPLICIT  REAL*8           (A-H,O-Z)
  54
  55 REAL(kind=dp), INTENT(inout)                        :: tk
  56 INTEGER :: test !!,AND
  57
  58 COMMON/alpha/mk,ik,im,ml,il,mruns,nruns,jr,iovec,npfs(8,3,47)
  59 COMMON/beta/tic,times(8,3,47),see(5,3,8,3),terrs(8,3,47),csums(8,3  &
  60     ,47),fopn(8,3,47),dos(8,3,47)
  61
  62 COMMON/spaces/ion,j5,k2,k3,loop1,laps,loop,m,kr,lp,n13h,ibuf,nx,l,  &
  63     npass,nfail,n,n1,n2,n13,n213,n813,n14,n16,n416,n21,nt1,nt2,last,idebug  &
  64     ,mpy,loop2,mucho,mpylim,intbuf(16)
  65
  66 COMMON/spacer/a11,a12,a13,a21,a22,a23,a31,a32,a33,ar,br,c0,cr,di,dk  &
  67     ,dm22,dm23,dm24,dm25,dm26,dm27,dm28,dn,e3,e6,expmax,flx,q,qa,r,ri  &
  68     ,s,scale,sig,stb5,t,xnc,xnei,xnm
  69
  70 COMMON/space0/time(47),csum(47),ww(47),wt(47),ticks,fr(9),terr1(47  &
  71     ),sumw(7),start,skale(47),bias(47),ws(95),total(47),flopn(47),iq(7  &
  72     ),npf,npfs1(47)
  73
  74 COMMON/spacei/wtp(3),mul(3),ispan(47,3),ipass(47,3)
  75
  76 !     ******************************************************************
  77
  78
  79 INTEGER :: e,f,zone
  80 COMMON/ispace/e(96),f(96),ix(1001),ir(1001),zone(300)
  81
  82 COMMON/space1/u(1001),v(1001),w(1001),x(1001),y(1001),z(1001),g(1001)  &
  83     ,du1(101),du2(101),du3(101),grd(1001),dex(1001),xi(1001),ex(1001)  &
  84     ,ex1(1001),dex1(1001),vx(1001),xx(1001),rx(1001),rh(2048),vsp(101)  &
  85     ,vstp(101),vxne(101),vxnd(101),ve3(101),vlr(101),vlin(101),b5(101)  &
  86     ,plan(300),d(300),sa(101),sb(101)
  87
  88 COMMON/space2/p(4,512),px(25,101),cx(25,101),vy(101,25),vh(101,7),  &
  89     vf(101,7),vg(101,7),vs(101,7),za(101,7),zp(101,7),zq(101,7),zr(101  &
  90     ,7),zm(101,7),zb(101,7),zu(101,7),zv(101,7),zz(101,7),b(64,64),c(64,64)  &
  91     ,h(64,64),u1(5,101,2),u2(5,101,2),u3(5,101,2)
  92
  93 !     ******************************************************************
  94
  95 dimension zx(1023),xz(447,3),tk(6),mtmp(1)
  96 EQUIVALENCE(zx(1),z(1)),(xz(1,1),x(1))
  97 double precision temp
  98 logical ltmp
  99
 100
 101 !     ******************************************************************
 102
 103 !     STANDARD PRODUCT COMPILER DIRECTIVES MAY BE USED FOR OPTIMIZATION
 104
 105
 106
 107
 108
 109 CALL trace('KERNEL  ')
 110
 111 CALL SPACE
 112
 113 mpy= 1
 114 mpysav= mpylim
 115 loop2= 1
 116 mpylim= loop2
 117 l= 1
 118 loop= 1
 119 lp= loop
 120 it0= test(0)
 121 loop2= mpysav
 122 mpylim= loop2
 123 do
 124
 125 !***********************************************************************
 126 !***  KERNEL 1      HYDRO FRAGMENT
 127 !***********************************************************************
 128
 129   x(:n)= q+y(:n)*(r*zx(11:n+10)+t*zx(12:n+11))
 130 IF(test(1) <= 0)THEN
 131   EXIT
 132 END IF
 133 END DO
 134
 135 do
 136 !                   we must execute    DO k= 1,n  repeatedly for accurat
 137 !    e timing
 138
 139 !***********************************************************************
 140 !***  KERNEL 2      ICCG EXCERPT (INCOMPLETE CHOLESKY - CONJUGATE GRADIE
 141 !    NT)
 142 !***********************************************************************
 143
 144
 145 ii= n
 146 ipntp= 0
 147
 148 do while(ii >  1)
 149 ipnt= ipntp
 150 ipntp= ipntp+ii
 151 ii= ishft(ii,-1)
 152 i= ipntp+1
 153 !dir$ vector always
 154        x(ipntp+2:ipntp+ii+1)=x(ipnt+2:ipntp:2)-v(ipnt+2:ipntp:2) &
 155      &*x(ipnt+1:ipntp-1:2)-v(ipnt+3:ipntp+1:2)*x(ipnt+3:ipntp+1:2)
 156 END DO
 157 IF(test(2) <= 0)THEN
 158   EXIT
 159 END IF
 160 END DO
 161
 162 do
 163
 164 !***********************************************************************
 165 !***  KERNEL 3      INNER PRODUCT
 166 !***********************************************************************
 167
 168
 169 q= dot_product(z(:n),x(:n))
 170 IF(test(3) <= 0)THEN
 171   EXIT
 172 END IF
 173 END DO
 174 m= (1001-7)/2
 175
 176 !***********************************************************************
 177 !***  KERNEL 4      BANDED LINEAR EQUATIONS
 178 !***********************************************************************
 179
 180 fw= 1.000D-25
 181
 182 do
 183 !dir$ vector always
 184  xz(6,:3)= y(5)*(xz(6,:3)+matmul(y(5:n:5), xz(:n/5,:3)))
 185
 186 IF(test(4) <= 0)THEN
 187   EXIT
 188 END IF
 189 END DO
 190
 191 do
 192
 193 !***********************************************************************
 194 !***  KERNEL 5      TRI-DIAGONAL ELIMINATION, BELOW DIAGONAL (NO VECTORS
 195 !    )
 196 !***********************************************************************
 197
 198
 199 tmp= x(1)
 200 DO i= 2,n
 201   tmp= z(i)*(y(i)-tmp)
 202   x(i)= tmp
 203 END DO
 204 IF(test(5) <= 0)THEN
 205   EXIT
 206 END IF
 207 END DO
 208
 209 do
 210
 211 !***********************************************************************
 212 !***  KERNEL 6      GENERAL LINEAR RECURRENCE EQUATIONS
 213 !***********************************************************************
 214
 215
 216 DO i= 2,n
 217   w(i)= 0.0100D0+dot_product(b(i,:i-1),w(i-1:1:-1))
 218 END DO
 219 IF(test(6) <= 0)THEN
 220   EXIT
 221 END IF
 222 END DO
 223
 224 do
 225
 226 !***********************************************************************
 227 !***  KERNEL 7      EQUATION OF STATE FRAGMENT
 228 !***********************************************************************
 229
 230
 231   x(:n)= u(:n)+r*(z(:n)+r*y(:n))+t*(u(4:n+3)+r*(u(3:n+2)+r*u(2:n+1))+t*(  &
 232       u(7:n+6)+q*(u(6:n+5)+q*u(5:n+4))))
 233 IF(test(7) <= 0)THEN
 234   EXIT
 235 END IF
 236 END DO
 237
 238 do
 239
 240
 241 !***********************************************************************
 242 !***  KERNEL 8      A.D.I. INTEGRATION
 243 !***********************************************************************
 244
 245
 246 nl1= 1
 247 nl2= 2
 248 fw= 2.000D0
 249   DO ky= 2,n
 250 DO kx= 2,3
 251     du1ky= u1(kx,ky+1,nl1)-u1(kx,ky-1,nl1)
 252     du2ky= u2(kx,ky+1,nl1)-u2(kx,ky-1,nl1)
 253     du3ky= u3(kx,ky+1,nl1)-u3(kx,ky-1,nl1)
 254     u1(kx,ky,nl2)= u1(kx,ky,nl1)+a11*du1ky+a12*du2ky+a13  &
 255         *du3ky+sig*(u1(kx+1,ky,nl1)-fw*u1(kx,ky,nl1)+u1(kx-1,ky,nl1))
 256     u2(kx,ky,nl2)= u2(kx,ky,nl1)+a21*du1ky+a22*du2ky+a23  &
 257         *du3ky+sig*(u2(kx+1,ky,nl1)-fw*u2(kx,ky,nl1)+u2(kx-1,ky,nl1))
 258     u3(kx,ky,nl2)= u3(kx,ky,nl1)+a31*du1ky+a32*du2ky+a33  &
 259         *du3ky+sig*(u3(kx+1,ky,nl1)-fw*u3(kx,ky,nl1)+u3(kx-1,ky,nl1))
 260   END DO
 261 END DO
 262 IF(test(8) <= 0)THEN
 263   EXIT
 264 END IF
 265 END DO
 266
 267 do
 268
 269 !***********************************************************************
 270 !***  KERNEL 9      INTEGRATE PREDICTORS
 271 !***********************************************************************
 272
 273
 274   px(1,:n)= dm28*px(13,:n)+px(3,:n)+dm27*px(12,:n)+dm26*px(11,:n)+dm25*px(10  &
 275       ,:n)+dm24*px(9,:n)+dm23*px(8,:n)+dm22*px(7,:n)+c0*(px(5,:n)+px(6,:n))
 276 IF(test(9) <= 0)THEN
 277   EXIT
 278 END IF
 279 END DO
 280
 281 do
 282
 283 !***********************************************************************
 284 !***  KERNEL 10     DIFFERENCE PREDICTORS
 285 !***********************************************************************
 286
 287 !dir$ unroll(2)
 288           do k= 1,n
 289               br= cx(5,k)-px(5,k)
 290               px(5,k)= cx(5,k)
 291               cr= br-px(6,k)
 292               px(6,k)= br
 293               ar= cr-px(7,k)
 294               px(7,k)= cr
 295               br= ar-px(8,k)
 296               px(8,k)= ar
 297               cr= br-px(9,k)
 298               px(9,k)= br
 299               ar= cr-px(10,k)
 300               px(10,k)= cr
 301               br= ar-px(11,k)
 302               px(11,k)= ar
 303               cr= br-px(12,k)
 304               px(12,k)= br
 305               px(14,k)= cr-px(13,k)
 306               px(13,k)= cr
 307             enddo
 308 IF(test(10) <= 0)THEN
 309   EXIT
 310 END IF
 311 END DO
 312
 313 do
 314
 315 !***********************************************************************
 316 !***  KERNEL 11     FIRST SUM.   PARTIAL SUMS.              (NO VECTORS)
 317 !***********************************************************************
 318
 319
 320 temp= 0
 321 DO k= 1,n
 322   temp= temp+y(k)
 323   x(k)= temp
 324 END DO
 325 IF(test(11) <= 0)THEN
 326   EXIT
 327 END IF
 328 END DO
 329
 330 do
 331
 332 !***********************************************************************
 333 !***  KERNEL 12     FIRST DIFF.
 334 !***********************************************************************
 335
 336   x(:n)= y(2:n+1)-y(:n)
 337 IF(test(12) <= 0)THEN
 338   EXIT
 339 END IF
 340 END DO
 341 fw= 1.000D0
 342
 343 !***********************************************************************
 344 !***  KERNEL 13      2-D PIC   Particle In Cell
 345 !***********************************************************************
 346
 347
 348 do
 349
 350 ! rounding modes for integerizing make no difference here
 351           do k= 1,n
 352               i1= 1+iand(int(p(1,k)),63)
 353               j1= 1+iand(int(p(2,k)),63)
 354               p(3,k)= p(3,k)+b(i1,j1)
 355               p(1,k)= p(1,k)+p(3,k)
 356               i2= iand(int(p(1,k)),63)
 357               p(1,k)= p(1,k)+y(i2+32)
 358               p(4,k)= p(4,k)+c(i1,j1)
 359               p(2,k)= p(2,k)+p(4,k)
 360               j2= iand(int(p(2,k)),63)
 361               p(2,k)= p(2,k)+z(j2+32)
 362               i2= i2+e(i2+32)
 363               j2= j2+f(j2+32)
 364               h(i2,j2)= h(i2,j2)+fw
 365             enddo
 366 IF(test(13) <= 0)THEN
 367   EXIT
 368 END IF
 369 END DO
 370 fw= 1.000D0
 371
 372 !***********************************************************************
 373 !***  KERNEL 14      1-D PIC   Particle In Cell
 374 !***********************************************************************
 375
 376
 377
 378 do
 379
 380   ix(:n)= grd(:n)
 381 !dir$ ivdep
 382   vx(:n)= ex(ix(:n))-ix(:n)*dex(ix(:n))
 383   ir(:n)= vx(:n)+flx
 384   rx(:n)= vx(:n)+flx-ir(:n)
 385   ir(:n)= iand(ir(:n),2047)+1
 386   xx(:n)= rx(:n)+ir(:n)
 387 DO k= 1,n
 388       rh(ir(k))= rh(ir(k))+fw-rx(k)
 389       rh(ir(k)+1)= rh(ir(k)+1)+rx(k)
 390 END DO
 391 IF(test(14) <= 0)THEN
 392   EXIT
 393 END IF
 394 END DO
 395
 396 do
 397
 398 !***********************************************************************
 399 !***  KERNEL 15     CASUAL FORTRAN.  DEVELOPMENT VERSION.
 400 !***********************************************************************
 401
 402
 403 !       CASUAL ORDERING OF SCALAR OPERATIONS IS TYPICAL PRACTICE.
 404 !       THIS EXAMPLE DEMONSTRATES THE NON-TRIVIAL TRANSFORMATION
 405 !       REQUIRED TO MAP INTO AN EFFICIENT MACHINE IMPLEMENTATION.
 406
 407
 408 ng= 7
 409 nz= n
 410 ar= 0.05300D0
 411 br= 0.07300D0
 412 !$omp parallel do private(t,j,k,r,s,i,ltmp) if(nz>98)
 413 do j= 2,ng-1
 414   do k= 2,nz
 415     i= merge(k-1,k,vf(k,j) <  vf((k-1),j))
 416     t= merge(br,ar,vh(k,(j+1)) <= vh(k,j))
 417     r= MAX(vh(i,j),vh(i,j+1))
 418     s= vf(i,j)
 419     vy(k,j)= t/s*SQRT(vg(k,j)**2+r*r)
 420     if(k < nz)then
 421         ltmp=vf(k,j) >= vf(k,(j-1))
 422         i= merge(j,j-1,ltmp)
 423         t= merge(ar,br,ltmp)
 424         r= MAX(vg(k,i),vg(k+1,i))
 425         s= vf(k,i)
 426         vs(k,j)= t/s*SQRT(vh(k,j)**2+r*r)
 427     endif
 428   END do
 429   vs(nz,j)= 0.0D0
 430 END do
 431   vy(2:nz,ng)= 0.0D0
 432 IF(test(15) <= 0)THEN
 433   EXIT
 434 END IF
 435 END DO
 436 ii= n/3
 437
 438 !***********************************************************************
 439 !***  KERNEL 16     MONTE CARLO SEARCH LOOP
 440 !***********************************************************************
 441
 442 lb= ii+ii
 443 k2= 0
 444 k3= 0
 445
 446 do
 447 DO m= 1,zone(1)
 448   j2= (n+n)*(m-1)+1
 449   DO k= 1,n
 450     k2= k2+1
 451     j4= j2+k+k
 452     j5= zone(j4)
 453     IF(j5 >= n)THEN
 454       IF(j5 == n)THEN
 455         EXIT
 456       END IF
 457       k3= k3+1
 458       IF(d(j5) <  d(j5-1)*(t-d(j5-2))**2+(s-d(j5-3))**2+ (r-d(j5-4))**2)THEN
 459         go to 200
 460       END IF
 461       IF(d(j5) == d(j5-1)*(t-d(j5-2))**2+(s-d(j5-3))**2+ (r-d(j5-4))**2)THEN
 462         EXIT
 463       END IF
 464     ELSE
 465       IF(j5-n+lb <  0)THEN
 466         IF(plan(j5) <  t)THEN
 467           go to 200
 468         END IF
 469         IF(plan(j5) == t)THEN
 470           EXIT
 471         END IF
 472       ELSE
 473         IF(j5-n+ii <  0)THEN
 474           IF(plan(j5) <  s)THEN
 475             go to 200
 476           END IF
 477           IF(plan(j5) == s)THEN
 478             EXIT
 479           END IF
 480         ELSE
 481           IF(plan(j5) <  r)THEN
 482             go to 200
 483           END IF
 484           IF(plan(j5) == r)THEN
 485             EXIT
 486           END IF
 487         END IF
 488       END IF
 489     END IF
 490     IF(zone(j4-1) <= 0)THEN
 491       go to 200
 492     END IF
 493   END DO
 494   EXIT
 495   200             IF(zone(j4-1) == 0)THEN
 496     EXIT
 497   END IF
 498 END DO
 499 IF(test(16) <= 0)THEN
 500   EXIT
 501 END IF
 502 END DO
 503 dw= 5.0000D0/3.0000D0
 504
 505 !***********************************************************************
 506 !***  KERNEL 17     IMPLICIT, CONDITIONAL COMPUTATION       (NO VECTORS)
 507 !***********************************************************************
 508
 509 !         RECURSIVE-DOUBLING VECTOR TECHNIQUES CAN NOT BE USED
 510 !         BECAUSE CONDITIONAL OPERATIONS APPLY TO EACH ELEMENT.
 511
 512 fw= 1.0000D0/3.0000D0
 513 tw= 1.0300D0/3.0700D0
 514
 515 do
 516 scale= dw
 517 rtmp= fw
 518 e6= tw
 519 DO k= n,2,-1
 520   e3= rtmp*vlr(k)+vlin(k)
 521   xnei= vxne(k)
 522   vxnd(k)= e6
 523   xnc= scale*e3
 524 !                                      SELECT MODEL
 525   IF(max(rtmp,xnei) <= xnc)THEN
 526 !                                      LINEAR MODEL
 527     ve3(k)= e3
 528     rtmp= e3+e3-rtmp
 529     vxne(k)= e3+e3-xnei
 530   ELSE
 531     rtmp= rtmp*vsp(k)+vstp(k)
 532 !                                      STEP MODEL
 533     vxne(k)= rtmp
 534     ve3(k)= rtmp
 535   END IF
 536     e6= rtmp
 537 END DO
 538 xnm= rtmp
 539 IF(test(17) <= 0)THEN
 540   EXIT
 541 END IF
 542 END DO
 543
 544 do
 545
 546 !***********************************************************************
 547 !***  KERNEL 18     2-D EXPLICIT HYDRODYNAMICS FRAGMENT
 548 !***********************************************************************
 549
 550
 551 t= 0.003700D0
 552 s= 0.004100D0
 553 kn= 6
 554 jn= n
 555   zb(2:jn,2:kn)=(zr(2:jn,2:kn)+zr(2:jn,:kn-1))/(zm(2:jn,2:kn)+zm(:jn-1,2:kn)) &
 556         *(zp(:jn-1,2:kn)-zp(2:jn,2:kn)+(zq(:jn-1,2:kn)-zq(2:jn,2:kn)))
 557   za(2:jn,2:kn)=(zr(2:jn,2:kn)+zr(:jn-1,2:kn))/(zm(:jn-1,2:kn)+zm(:jn-1,3:kn+1))  &
 558         *(zp(:jn-1,3:kn+1)-zp(:jn-1,2:kn)+(zq(:jn-1,3:kn+1)-zq(:jn-1,2:kn)))
 559   zu(2:jn,2:kn)= zu(2:jn,2:kn)+ &
 560         s*(za(2:jn,2:kn)*(zz(2:jn,2:kn)-zz(3:jn+1,2:kn)) &
 561         -za(:jn-1,2:kn)*(zz(2:jn,2:kn)-zz(:jn-1,2:kn)) &
 562         -zb(2:jn,2:kn)*(zz(2:jn,2:kn)-zz(2:jn,:kn-1))+ &
 563         zb(2:jn,3:kn+1)*(zz(2:jn, 2:kn)-zz(2:jn,3:kn+1)))
 564   zv(2:jn,2:kn)= zv(2:jn,2:kn)+ &
 565         s*(za(2:jn,2:kn)*(zr(2:jn,2:kn)-zr(3:jn+1,2:kn)) &
 566         -za(:jn-1,2:kn)*(zr(2:jn,2:kn)-zr(:jn-1,2:kn)) &
 567         -zb(2:jn,2:kn)*(zr(2:jn,2:kn)-zr(2:jn,:kn-1))+ &
 568         zb(2:jn,3:kn+1)*(zr(2:jn, 2:kn)-zr(2:jn,3:kn+1)))
 569   zr(2:jn,2:kn)= zr(2:jn,2:kn)+t*zu(2:jn,2:kn)
 570   zz(2:jn,2:kn)= zz(2:jn,2:kn)+t*zv(2:jn,2:kn)
 571 IF(test(18) <= 0)THEN
 572   EXIT
 573 END IF
 574 END DO
 575
 576 do
 577
 578 !***********************************************************************
 579 !***  KERNEL 19      GENERAL LINEAR RECURRENCE EQUATIONS    (NO VECTORS)
 580 !***********************************************************************
 581
 582 kb5i= 0
 583
 584 DO k= 1,n
 585   b5(k+kb5i)= sa(k)+stb5*sb(k)
 586   stb5= b5(k+kb5i)-stb5
 587 END DO
 588 DO k= n,1,-1
 589   b5(k+kb5i)= sa(k)+stb5*sb(k)
 590   stb5= b5(k+kb5i)-stb5
 591 END DO
 592 IF(test(19) <= 0)THEN
 593   EXIT
 594 END IF
 595 END DO
 596 dw= 0.200D0
 597
 598 !***********************************************************************
 599 !***  KERNEL 20     DISCRETE ORDINATES TRANSPORT: RECURRENCE (NO VECTORS
 600 !***********************************************************************
 601
 602
 603 do
 604
 605 rtmp= xx(1)
 606 DO k= 1,n
 607   di= y(k)*(rtmp+dk)-g(k)
 608   dn=merge( max(s,min(z(k)*(rtmp+dk)/di,t)),dw,di /= 0.0)
 609   x(k)= ((w(k)+v(k)*dn)*rtmp+u(k))/(vx(k)+v(k)*dn)
 610   rtmp= ((w(k)-vx(k))*rtmp+u(k))*DN/(vx(k)+v(k)*dn)+ rtmp
 611  xx(k+1)= rtmp
 612 END DO
 613 IF(test(20) <= 0)THEN
 614   EXIT
 615 END IF
 616 END DO
 617
 618 do
 619
 620 !***********************************************************************
 621 !***  KERNEL 21     MATRIX*MATRIX PRODUCT
 622 !***********************************************************************
 623
 624     px(:25,:n)= px(:25,:n)+matmul(vy(:25,:25),cx(:25,:n))
 625 IF(test(21) <= 0)THEN
 626   EXIT
 627 END IF
 628 END DO
 629 expmax= 20.0000D0
 630
 631
 632 !***********************************************************************
 633 !***  KERNEL 22     PLANCKIAN DISTRIBUTION
 634 !***********************************************************************
 635
 636 !      EXPMAX= 234.500d0
 637 fw= 1.00000D0
 638 u(n)= 0.99000D0*expmax*v(n)
 639
 640 do
 641
 642   y(:n)= u(:n)/v(:n)
 643   w(:n)= x(:n)/(EXP(y(:n))-fw)
 644 IF(test(22) <= 0)THEN
 645   EXIT
 646 END IF
 647 END DO
 648 fw= 0.17500D0
 649
 650 !***********************************************************************
 651 !***  KERNEL 23     2-D IMPLICIT HYDRODYNAMICS FRAGMENT
 652 !***********************************************************************
 653
 654
 655 do
 656
 657       DO k= 2,n
 658          do j=2,6
 659              za(k,j)= za(k,j)+fw*(za(k,j+1)*zr(k,j)-za(k,j)+            &
 660      &          zv(k,j)*za(k-1,j)+(zz(k,j)+za(k+1,j)*                   &
 661      &          zu(k,j)+za(k,j-1)*zb(k,j)))
 662       END DO
 663     END DO
 664 IF(test(23) <= 0)THEN
 665   EXIT
 666 END IF
 667 END DO
 668 x(n/2)= -1.000D+10
 669
 670 !***********************************************************************
 671 !***  KERNEL 24     FIND LOCATION OF FIRST MINIMUM IN ARRAY
 672 !***********************************************************************
 673
 674 !      X( n/2)= -1.000d+50
 675
 676 do
 677  m= minloc(x(:n),DIM=1)
 678
 679 IF(test(24) == 0)THEN
 680   EXIT
 681 END IF
 682 END DO
 683 sum= 0.00D0
 684 som= 0.00D0
 685 DO k= 1,mk
 686   sum= sum+time(k)
 687   times(jr,il,k)= time(k)
 688   terrs(jr,il,k)= terr1(k)
 689   npfs(jr,il,k)= npfs1(k)
 690   csums(jr,il,k)= csum(k)
 691   dos(jr,il,k)= total(k)
 692   fopn(jr,il,k)= flopn(k)
 693   som= som+flopn(k)*total(k)
 694 END DO
 695 tk(1)= tk(1)+sum
 696 tk(2)= tk(2)+som
 697 !                        Dumpout Checksums:  file "chksum"
 698 !     WRITE ( 7,706) jr, il
 699 ! 706 FORMAT(1X,2I3)
 700 !     WRITE ( 7,707) ( CSUM(k), k= 1,mk)
 701 ! 707 FORMAT(5X,'&',1PE23.16,',',1PE23.16,',',1PE23.16,',')
 702
 703 CALL track('KERNEL  ')
 704 RETURN
 705 END SUBROUTINE kernel