From 8508e65d32006029b4defef9c47e0f3c7f707859 Mon Sep 17 00:00:00 2001 From: Ben Dudson Date: Wed, 19 Jul 2023 11:23:19 -0700 Subject: [PATCH 01/13] Removed shared labels in do loops Fixes warnings from grd/grdcomp like ``` Fortran 2018 deleted feature: Shared DO termination label 2208 at (1) ``` --- grd/grdcomp.m | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/grd/grdcomp.m b/grd/grdcomp.m index e4926226..52a8d6d7 100755 --- a/grd/grdcomp.m +++ b/grd/grdcomp.m @@ -642,33 +642,37 @@ call writedata ('gridue', runmirror) c loop over all regions -- do 800 region=1,noregs - do 2208 j=jmin(region),jmax(region) - do 2208 l=1,mseg - do 2208 k=1,npts + do j=jmin(region),jmax(region) + do l=1,mseg + do k=1,npts splcoef(k,l,j)=0.0 xknts(k,l,j)=0.0 - 2208 continue + end do + end do + end do c - do 2212 j=jmin(region),jmax(region) - do 2212 l=1,npts + do j=jmin(region),jmax(region) + do l=1,npts isegment(l,j)=0 - 2212 continue + end do + end do c - do 2209 j=jmin(region),jmax(region) + do j=jmin(region),jmax(region) nseg(j)=0 - do 2209 l=1,mseg + do l=1,mseg isys(l,j)=0 istartg(l,j)=0 iendg(l,j)=0 ncap7(l,j)=0 m(l,j)=0 - 2209 continue + end do + end do c - do 2210 k=1,npts + do k=1,npts wg(k)=1.0 xtrans(k)=0.0 ytrans(k)=0.0 - 2210 continue + end do 800 continue # end of loop over regions return From 8229a393696d83b6bb74d497fc91f8debe356c84 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 13:15:29 -0700 Subject: [PATCH 02/13] Remove DO termination statement warning from daux1 --- svr/daux1.f | 20 ++++++++++++-------- 1 file changed, 12 insertions(+), 8 deletions(-) diff --git a/svr/daux1.f b/svr/daux1.f index ef370202..6d9f025e 100755 --- a/svr/daux1.f +++ b/svr/daux1.f @@ -13,20 +13,24 @@ subroutine dewset_u (n, itol, rtol, atol, ycur, ewt) c go to (10, 20, 30, 40), itol 10 continue - do 15 i = 1,n - 15 ewt(i) = rtol(1)*dabs(ycur(i)) + atol(1) + do i = 1,n + ewt(i) = rtol(1)*dabs(ycur(i)) + atol(1) + end do return 20 continue - do 25 i = 1,n - 25 ewt(i) = rtol(1)*dabs(ycur(i)) + atol(i) + do i = 1,n + ewt(i) = rtol(1)*dabs(ycur(i)) + atol(i) + end do return 30 continue - do 35 i = 1,n - 35 ewt(i) = rtol(i)*dabs(ycur(i)) + atol(1) + do i = 1,n + ewt(i) = rtol(i)*dabs(ycur(i)) + atol(1) + end do return 40 continue - do 45 i = 1,n - 45 ewt(i) = rtol(i)*dabs(ycur(i)) + atol(i) + do i = 1,n + ewt(i) = rtol(i)*dabs(ycur(i)) + atol(i) + end do return c----------------------- end of subroutine dewset ---------------------- end From f7d0ad1e7bfdb9f9c8c8f9afb486d4fe4f011c09 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 13:20:36 -0700 Subject: [PATCH 03/13] Fixes Shared DO termination label warnings in fmombal --- api/fmombal.m | 128 +++++++++++++++++++++++++++----------------------- 1 file changed, 68 insertions(+), 60 deletions(-) diff --git a/api/fmombal.m b/api/fmombal.m index 2970b841..ea65a010 100755 --- a/api/fmombal.m +++ b/api/fmombal.m @@ -402,7 +402,7 @@ c AND U(1,1,alpha) MUST BE EXPRESSED IN TERMS OF AI, UBAR amat(1+j,1+j) = one !totmass*anorm * fanorm amat(2+j,2+j) = one !NOT USED FOR ANYTHING PHYSICAL amat(3+j,3+j) = one !NOT USED FOR ANYTHING PHYSICAL - do 10 misa = 1,miso + do misa = 1,miso nz = natom(misa) dnur1 = (fanorm*denmass(misa,z1))*nurec(misa,z1) c @@ -424,14 +424,14 @@ c AND U(1,1,alpha) MUST BE EXPRESSED IN TERMS OF AI, UBAR > rhomass*uresp(mflow,nzsum,misa,iacci) enddo i = KXA*(misa-1) - do 10 m = 1,KXA + do m = 1,KXA if( (misa.ne.mise) .or. (m.ne.mflow) )then sbar(m,misa) = > sdot(nz,zi(misa,1),miso,uresp(m,z1,misa,iforc),KXA) amat(m+i,1+j) = > -sdot(nz,zi(misa,1),miso,uresp(m,z1,misa,iacci),KXA) endif - do 10 mpmisb = 1,KXA*miso + do mpmisb = 1,KXA*miso do nzb = 1,nz usum(nzb) = uresp(m,nzb,misa,ilam1)*elab(ilam1,misa,mpmisb) > + uresp(m,nzb,misa,ilam2)*elab(ilam2,misa,mpmisb) @@ -450,7 +450,9 @@ c PERFORM SUMS OVER MISA (ISOTOPE MASS A) > + (denmass(misa,nzsum)/totmass)*usum(nzsum) enddo endif - 10 continue + end do + end do + end do c c ADD DIAGONAL ELEMENTS c @@ -694,43 +696,45 @@ real nuion(miso,0:nzch), nurec(miso,nzch), acci = acci0*anorm sumflow = 0. sumaflow= 0. - do 55 misa = 1,miso + do misa = 1,miso vtherm=sqrt(2.0*tempa(misa)/amu(misa)/promas) - do 55 m = 1,2 !Only look at flow, heat flow - if( m.eq.1 )then - write(*,500)m,amu(misa),sbar(m,misa), - > 1000.*tempa(misa)/xj7kv,vtherm,caplam(m,misa) - write(*,510) - if( misa.ne.1 ) - > write(*,520)0,uneut(misa),den(misa,0),nuion(misa,0) - else if(m.eq.2)then - write(*,505)m,amu(misa),sbar(m,misa) - write(*,530) - endif - nzmax = natom(misa) - do 55 nz = 1,nzmax - force = fmom(m,nz,misa) - friction = zi(misa,nz) * ( ela(m,1,misa)*usol(1,nz,misa) - > + ela(m,2,misa)*usol(2,nz,misa) - > + ela(m,3,misa)*usol(3,nz,misa) + caplam(m,misa) ) - friction = friction - denmass(misa,nz)*usol(m,nz,misa)* - > (nuion(misa,nz) + nurec(misa,nz))*al32(m) - if( nz.gt.1 )friction = friction + usol(m,nz-1,misa) - > *denmass(misa,nz-1) * nuion(misa,nz-1) * al32(m) - if( nz.lt.nzmax )friction = friction + usol(m,nz+1,misa) - > *denmass(misa,nz+1) * nurec(misa,nz+1) * al32(m) - if( m.eq.1 )then - force = force + denmass(misa,nz)*acci - sumflow = sumflow + denmass(misa,nz)*usol(1,nz,misa) - sumaflow= sumaflow+ denmass(misa,nz)*abs(usol(1,nz,misa)) - write(*,550)nz,usol(m,nz,misa),force,friction, - > den(misa,nz),zi(misa,nz),nuion(misa,nz),nurec(misa,nz) - else if( m.eq.2 )then - pres = den(misa,nz)*tempa(misa) - write(*,550)nz,-2.5*usol(m,nz,misa)*pres, - > force,friction,gradp(misa,nz),den(misa,nz)*gradt(misa,nz) - endif - 55 continue + do m = 1,2 !Only look at flow, heat flow + if( m.eq.1 )then + write(*,500)m,amu(misa),sbar(m,misa), + > 1000.*tempa(misa)/xj7kv,vtherm,caplam(m,misa) + write(*,510) + if( misa.ne.1 ) + > write(*,520)0,uneut(misa),den(misa,0),nuion(misa,0) + else if(m.eq.2)then + write(*,505)m,amu(misa),sbar(m,misa) + write(*,530) + endif + nzmax = natom(misa) + do nz = 1,nzmax + force = fmom(m,nz,misa) + friction = zi(misa,nz) * ( ela(m,1,misa)*usol(1,nz,misa) + > + ela(m,2,misa)*usol(2,nz,misa) + > + ela(m,3,misa)*usol(3,nz,misa) + caplam(m,misa) ) + friction = friction - denmass(misa,nz)*usol(m,nz,misa)* + > (nuion(misa,nz) + nurec(misa,nz))*al32(m) + if( nz.gt.1 )friction = friction + usol(m,nz-1,misa) + > *denmass(misa,nz-1) * nuion(misa,nz-1) * al32(m) + if( nz.lt.nzmax )friction = friction + usol(m,nz+1,misa) + > *denmass(misa,nz+1) * nurec(misa,nz+1) * al32(m) + if( m.eq.1 )then + force = force + denmass(misa,nz)*acci + sumflow = sumflow + denmass(misa,nz)*usol(1,nz,misa) + sumaflow= sumaflow+ denmass(misa,nz)*abs(usol(1,nz,misa)) + write(*,550)nz,usol(m,nz,misa),force,friction, + > den(misa,nz),zi(misa,nz),nuion(misa,nz),nurec(misa,nz) + else if( m.eq.2 )then + pres = den(misa,nz)*tempa(misa) + write(*,550)nz,-2.5*usol(m,nz,misa)*pres, + > force,friction,gradp(misa,nz),den(misa,nz)*gradt(misa,nz) + endif + end do + end do + end do sumaflow = max(sumaflow,totmass*abs(umass)) 500 format(/' UBAR (m=',i1,', mass=',1pe9.2,') = ',1pe10.3, > ' (m/sec) TEMP = ',1pe10.3,' (eV) VTHERM = ',1pe10.3, @@ -940,7 +944,7 @@ call scopy(k2,a(napk+1),1,a(na+1),1) kblock = km1 krow = k endif - do 90 j=1,kblock + do j=1,kblock c***** Determine element of maximum modulus in next column. c***** Only allow permutations within blocks of k rows. if( iflag.ne.0 )goto 1010 @@ -992,7 +996,7 @@ call saxpy(krow-jp1+1,-sca,a(ms+jp1),1, c***** Subtract pivot row right side from other right-hand side(s) if(sca.ne.0.0)call saxpy(nbrt,-sca,b(msb,1),nk,b(ls,1),nk) enddo - 90 continue + end do c***** Move block 2 to 1, block 3 to 2, 0 out block 3 of rows used in c***** this stage that will also be used in next. if( iblock.ne.n .and. iflag.eq.0 )then @@ -1013,7 +1017,7 @@ call scopy(k2,a(nrow(l)+k+1),1,a(nrow(l)+1),1) iflag=1 return endif - do 100 lb = 1,nbrt + do lb = 1,nbrt nkb = nrow(nk)/k3 + 1 x(nk,lb)=b(nkb,lb)/a(nrow(nk)+k) neq=nk-1 @@ -1039,7 +1043,8 @@ call scopy(k2,a(nrow(l)+k+1),1,a(nrow(l)+1),1) neq=neq-1 enddo if(neq.gt.0) goto 180 - 100 continue + 100 continue + end do iflag=0 return end @@ -1111,7 +1116,7 @@ c SAME (SEE SUBROUTINE ZSOURCE FOR DETAILS) c kx3 = 3*KXA kxsq3 = kx3 * KXA - do 100 misa = 1,miso + do misa = 1,miso nz = natom(misa) nkz = nz * KXA nresp = NBA*nkz !Number of response array elements per isotope @@ -1146,11 +1151,11 @@ call zsource(asource,zi,den,fmom(1,1,misa),denmass,misa,nz) amat(i)=zero enddo - do 10 m = 1,KXA - do 10 mp= 1,KXA + do m = 1,KXA + do mp= 1,KXA i1 = mp + kx3*(m-1) i2 = m + nforc - do 10 jz= 1,nz + do jz= 1,nz index = i1 + (jz-1)*kxsq3 !consecutive block-by-block index is = i2 + (jz-1)*KXA c @@ -1187,7 +1192,9 @@ c UPPER DIAGONAL ELEMENTS (RECOMBINATION FROM HIGHER CHARGE STATE) > asource(is) = asource(is) - amat(index)*usave(mp,jp,misa) endif endif - 10 continue + end do + end do + end do c c SOLVE TRIDIAGONAL SYSTEM WITH NBA RIGHT SIDES STORED IN SOURCE c STORE FOR FAST (INCREMENTAL) RECALCULATION (LDIR = 0,1) @@ -1222,7 +1229,7 @@ c URESP(KXA=1,2,3,NATOM,n=1,..,NBA). noff = 1 + nkz*(ntype-1) call scopy(nkz,xsol(noff),1,uresp(1,1,misa,ntype),1) enddo - 100 continue + end do return end c---- End of subroutine zrespond --------------------------------------- @@ -1232,16 +1239,17 @@ implicit real (a-h,o-z), integer (i-n) Use(Reduced_ion_constants) real zi(miso,*), den(miso,0:nz), source(KXA,nz,NBA), fmom(KXA,*) real denmass(miso,0:nz) - do 20 jz= 1,nz - do 20 m = 1,KXA - if( m.eq.1 )then - source(m,jz,ilam1) = one - source(m,jz,iacci) = denmass(misa,jz) * anorm / zi(misa,jz) - else if( (m.eq.ilam2) .or. (m.eq.ilam3) )then - source(m,jz,m) = one - endif - source(m,jz,iforc) = fmom(m,jz)/zi(misa,jz) - 20 continue + do jz= 1,nz + do m = 1,KXA + if( m.eq.1 )then + source(m,jz,ilam1) = one + source(m,jz,iacci) = denmass(misa,jz) * anorm / zi(misa,jz) + else if( (m.eq.ilam2) .or. (m.eq.ilam3) )then + source(m,jz,m) = one + endif + source(m,jz,iforc) = fmom(m,jz)/zi(misa,jz) + end do + end do return end c---- End of subroutine zsource ---------------------------------------- From eef5ee660c0af637c07e45df176e83523ce7721e Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 13:46:10 -0700 Subject: [PATCH 04/13] Remove DO termination statement warning from vodpk --- svr/vodpk.m | 385 +++++++++++++++++++++++++++++++--------------------- 1 file changed, 228 insertions(+), 157 deletions(-) diff --git a/svr/vodpk.m b/svr/vodpk.m index 7caed32b..a8c3b1da 100755 --- a/svr/vodpk.m +++ b/svr/vodpk.m @@ -1708,9 +1708,10 @@ CALL SCOPY (N, Y, 1, RWORK(LYH), 1) NQ = 1 H = ONE CALL EWSET (N, ITOL, RTOL, ATOL, RWORK(LYH), RWORK(LEWT)) - DO 120 I = 1, N + DO I = 1, N IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 621 - 120 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) + RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) + END DO IF (H0 .NE. ZERO) GO TO 180 C Call VHIN to set initial step size H0 to be attempted. --------------- CALL VHIN (N, T, RWORK(LYH), RWORK(LF0), F, RPAR, IPAR, TOUT, @@ -1817,9 +1818,10 @@ CALL XERRWV (MSG, 56, 113, 0, 0, 0, 0, 0, ZERO, ZERO) CALL XERRWV (MSG, 56, 113, 0, 0, 0, 0, 2, TN, RCFL) ENDIF 255 CONTINUE - DO 260 I = 1, N + DO I = 1, N IF (RWORK(I+LEWT-1) .LE. ZERO) GO TO 510 - 260 RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) + RWORK(I+LEWT-1) = ONE/RWORK(I+LEWT-1) + END DO 270 TOLSF = UROUND*VNORML (N, RWORK(LYH), RWORK(LEWT)) IF (TOLSF .LE. ONE) GO TO 280 TOLSF = TOLSF*TWO @@ -1978,12 +1980,13 @@ CALL XERRWV(MSG, 50, 208, 0, 0, 0, 0, 2, TN, H) C Compute IMXER if relevant. ------------------------------------------- 560 BIG = ZERO IMXER = 1 - DO 570 I = 1, N + DO I = 1, N SIZE = ABS(RWORK(I+LACOR-1)*RWORK(I+LEWT-1)) IF (BIG .GE. SIZE) GO TO 570 BIG = SIZE IMXER = I 570 CONTINUE + END DO IWORK(16) = IMXER C Set Y vector, T, and optional outputs. ------------------------------- 580 CONTINUE @@ -2279,19 +2282,22 @@ CALL JAC (F, N, TN, Y, YH, EWT, SAVF, ACOR, HRL1, CRATE = ONE NSLP = NST IF (IERPJ .NE. 0) GO TO 420 - 250 DO 260 I = 1, N - 260 ACOR(I) = ZERO + 250 DO I = 1, N + ACOR(I) = ZERO + END DO 270 IF (MITER .NE. 0) GO TO 350 C----------------------------------------------------------------------- C In the case of functional iteration, update Y directly from C the result of the last function evaluation. C----------------------------------------------------------------------- - DO 290 I = 1, N + DO I = 1, N SAVF(I) = RL1*(H*SAVF(I) - YH(I,2)) - 290 Y(I) = SAVF(I) - ACOR(I) + Y(I) = SAVF(I) - ACOR(I) + END DO DEL = VNORML (N, Y, EWT) - DO 300 I = 1, N - 300 Y(I) = YH(I,1) + SAVF(I) + DO I = 1, N + Y(I) = YH(I,1) + SAVF(I) + END DO CALL SCOPY (N, SAVF, 1, ACOR, 1) GO TO 400 C----------------------------------------------------------------------- @@ -2300,8 +2306,9 @@ CALL SCOPY (N, SAVF, 1, ACOR, 1) C A as coefficient matrix. The correction is scaled by the factor C 2/(1+RC) to account for changes in H*RL1 since the last JAC call. C----------------------------------------------------------------------- - 350 DO 360 I = 1, N - 360 VSAV(I) = HRL1*SAVF(I) - (RL1*YH(I,2) + ACOR(I)) + 350 DO I = 1, N + VSAV(I) = HRL1*SAVF(I) - (RL1*YH(I,2) + ACOR(I)) + END DO CALL VSOLPK (Y, SAVF, VSAV, EWT, WM, IWM, F, PSOL, IERSL, 1 RPAR, IPAR) NNI = NNI + 1 @@ -2314,8 +2321,9 @@ CALL SSCAL (N, CSCALE, VSAV, 1) IF (IERSL .GT. 0) GO TO 410 DEL = VNORML (N, VSAV, EWT) CALL SAXPY (N, ONE, VSAV, 1, ACOR, 1) - DO 380 I = 1, N - 380 Y(I) = YH(I,1) + ACOR(I) + DO I = 1, N + Y(I) = YH(I,1) + ACOR(I) + END DO C----------------------------------------------------------------------- C Test for convergence. If M.gt.0, an estimate of the convergence C rate constant is stored in CRATE, and this is used in the test. @@ -2598,16 +2606,18 @@ C and V(*,k). C The initial residual is the vector b. Apply scaling to b, and test C for an immediate return with x = 0 or x = b. C----------------------------------------------------------------------- - DO 10 I = 1, N - 10 V(I,1) = B(I)*WGHT(I) + DO I = 1, N + V(I,1) = B(I)*WGHT(I) + END DO BNRM0 = SNRM2 (N, V, 1) BNRM = BNRM0 IF (BNRM0 .GT. DELTA) GO TO 30 IF (MNEWT .GT. 0) GO TO 20 CALL SCOPY (N, B, 1, X, 1) RETURN - 20 DO 25 I = 1, N - 25 X(I) = 0.0E0 + 20 DO I = 1, N + X(I) = 0.0E0 + END DO RETURN 30 CONTINUE C Apply inverse of left preconditioner to vector b. -------------------- @@ -2619,17 +2629,19 @@ CALL PSOL (N, TN, Y, SAVF, WK, HB0, WP, IWP, B, 1, IF (IER .NE. 0) GO TO 300 40 IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 55 C Calculate norm of scaled vector V(*, 1) and normalize it. ------------ - DO 50 I = 1, N - 50 V(I,1) = B(I)*WGHT(I) + DO I = 1, N + V(I,1) = B(I)*WGHT(I) + END DO BNRM = SNRM2 (N, V, 1) DELTA = DELTA*(BNRM/BNRM0) 55 TEM = 1.0E0/BNRM CALL SSCAL (N, TEM, V(1,1), 1) C Zero out the HES array. ---------------------------------------------- - DO 65 J = 1, MAXL - DO 60 I = 1, MAXLP1 - 60 HES (I,J) = 0.0E0 - 65 CONTINUE + DO J = 1, MAXL + DO I = 1, MAXLP1 + HES (I,J) = 0.0E0 + END DO + END DO C----------------------------------------------------------------------- C Main loop to compute the vectors V(*,2) to V(*,MAXL). C The running product PROD is needed for the convergence test. @@ -2662,20 +2674,22 @@ CALL SHEQR (HES , MAXLP1, LL, Q, INFO, LL) IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN IF (LL .EQ. KMP+1) THEN CALL SCOPY (N, V(1,1), 1, DL, 1) - DO 75 I = 1, KMP + DO I = 1, KMP IP1 = I + 1 I2 = I*2 S = Q(I2) C = Q(I2-1) - DO 70 K = 1, N - 70 DL(K) = S*DL(K) + C*V(K,IP1) - 75 CONTINUE + DO K = 1, N + DL(K) = S*DL(K) + C*V(K,IP1) + END DO + END DO ENDIF S = Q(2*LL) C = Q(2*LL-1)/SNORMW LLP1 = LL + 1 - DO 80 K = 1, N - 80 DL(K) = S*DL(K) + C*V(K,LLP1) + DO K = 1, N + DL(K) = S*DL(K) + C*V(K,LLP1) + END DO DLNRM = SNRM2 (N, DL, 1) RHO = RHO*DLNRM ENDIF @@ -2706,17 +2720,20 @@ CALL SSCAL (N, TEM, V(1,LL+1), 1) 200 CONTINUE LL = LGMR LLP1 = LL + 1 - DO 210 K = 1, LLP1 - 210 B(K) = 0.0E0 + DO K = 1, LLP1 + B(K) = 0.0E0 + END DO B(1) = BNRM CALL SHELS (HES , MAXLP1, LL, Q, B) - DO 220 K = 1, N - 220 X(K) = 0.0E0 - DO 230 I = 1, LL + DO K = 1, N + X(K) = 0.0E0 + END DO + DO I = 1, LL CALL SAXPY (N, B(I), V(1,I), 1, X, 1) - 230 CONTINUE - DO 240 I = 1, N - 240 X(I) = X(I)/WGHT(I) + END DO + DO I = 1, N + X(I) = X(I)/WGHT(I) + END DO IF (JPRE .LE. 1) RETURN CALL PSOL (N, TN, Y, SAVF, WK, HB0, WP, IWP, X, 2, 1 IER, RPAR, IPAR) @@ -2844,14 +2861,16 @@ DIMENSION Y(*), SAVF(*), V(*), WGHT(*), FTEM(*), Z(*), C----------------------------------------------------------------------- NFAIL = 0 C Set vtem = D * v. ---------------------------------------------------- - DO 10 I = 1, N - 10 VTEM(I) = V(I)/WGHT(I) + DO I = 1, N + VTEM(I) = V(I)/WGHT(I) + END DO IF (JPRE .GE. 2) GO TO 30 C C JPRE = 0 or 1. Save y in Z and increment Y by VTEM. ----------------- CALL SCOPY (N, Y, 1, Z, 1) - DO 20 I = 1, N - 20 Y(I) = Z(I) + VTEM(I)*scaleps + DO I = 1, N + Y(I) = Z(I) + VTEM(I)*scaleps + END DO SIG = scaleps RSIG = 1.0E0/scaleps FAC = HB0*RSIG @@ -2864,14 +2883,16 @@ CALL PSOL (N, TN, Y, SAVF, FTEM, HB0, WP, IWP, VTEM, 2, NPSL = NPSL + 1 IF (IERP .NE. 0) GO TO 100 C Calculate l-2 norm of (D-inverse) * VTEM. ---------------------------- - DO 40 I = 1, N - 40 Z(I) = VTEM(I)*WGHT(I) + DO I = 1, N + Z(I) = VTEM(I)*WGHT(I) + END DO RSIG = SNRM2 (N, Z, 1)/scaleps SIG = 1.0E0/RSIG C Save y in Z and increment Y by VTEM/norm. ---------------------------- CALL SCOPY (N, Y, 1, Z, 1) - 45 DO 50 I = 1, N - 50 Y(I) = Z(I) + VTEM(I)*SIG + 45 DO I = 1, N + Y(I) = Z(I) + VTEM(I)*SIG + END DO FAC = HB0*RSIG C C For all JPRE, call F with incremented Y argument, and restore Y. ----- @@ -2888,10 +2909,12 @@ CALL F (N, TN, Y, FTEM, RPAR, IPAR, IFAIL) ENDIF CALL SCOPY (N, Z, 1, Y, 1) C Set Z = (I - HB0*Jacobian) * VTEM, using difference quotient. -------- - DO 70 I = 1, N - 70 Z(I) = FTEM(I) - SAVF(I) - DO 80 I = 1, N - 80 Z(I) = VTEM(I) - FAC*Z(I) + DO I = 1, N + Z(I) = FTEM(I) - SAVF(I) + END DO + DO I = 1, N + Z(I) = VTEM(I) - FAC*Z(I) + END DO C Apply inverse of left preconditioner to Z, if nontrivial. ------------ IF (JPRE .EQ. 0 .OR. JPRE .EQ. 2) GO TO 85 CALL PSOL (N, TN, Y, SAVF, FTEM, HB0, WP, IWP, Z, 1, @@ -2900,8 +2923,9 @@ CALL PSOL (N, TN, Y, SAVF, FTEM, HB0, WP, IWP, Z, 1, IF (IERP .NE. 0) GO TO 100 85 CONTINUE C Apply D-inverse to Z and return. ------------------------------------- - DO 90 I = 1, N - 90 Z(I) = Z(I)*WGHT(I) + DO I = 1, N + Z(I) = Z(I)*WGHT(I) + END DO IER = 0 RETURN C @@ -2994,8 +3018,9 @@ DIMENSION Y(*), SAVF(*), B(*), WGHT(*), X(*), IF (MNEWT .GT. 0) GO TO 10 CALL SCOPY (N, B, 1, X, 1) RETURN - 10 DO 20 I = 1, N - 20 X(I) = 0.0E0 + 10 DO I = 1, N + X(I) = 0.0E0 + END DO RETURN C Apply inverse of left preconditioner to vector b. -------------------- 30 IER = 0 @@ -3075,40 +3100,52 @@ C listed (local) variables to be saved between calls to this integrator. 1 LRVK1 /3/, LIVK1 /11/ C IF (JOB .EQ. 2) GO TO 100 - DO 10 I = 1, LENRV1 - 10 RSAV(I) = RVOD1(I) - DO 12 I = 1, LENRV2 - 12 RSAV(LENRV1+I) = RVOD2(I) + DO I = 1, LENRV1 + RSAV(I) = RVOD1(I) + END DO + DO I = 1, LENRV2 + RSAV(LENRV1+I) = RVOD2(I) + END DO IOFF = LENRV1 + LENRV2 - DO 14 I = 1, LRVK1 - 14 RSAV(IOFF+I) = RVPK1(I) -C - DO 20 I = 1, LENIV1 - 20 ISAV(I) = IVOD1(I) - DO 22 I = 1, LENIV2 - 22 ISAV(LENIV1+I) = IVOD2(I) + DO I = 1, LRVK1 + RSAV(IOFF+I) = RVPK1(I) + END DO +C + DO I = 1, LENIV1 + ISAV(I) = IVOD1(I) + END DO + DO I = 1, LENIV2 + ISAV(LENIV1+I) = IVOD2(I) + END DO IOFF = LENIV1 + LENIV2 - DO 24 I = 1, LIVK1 - 24 ISAV(IOFF+I) = IVPK1(I) + DO I = 1, LIVK1 + ISAV(IOFF+I) = IVPK1(I) + END DO C RETURN C 100 CONTINUE - DO 110 I = 1, LENRV1 - 110 RVOD1(I) = RSAV(I) - DO 112 I = 1, LENRV2 - 112 RVOD2(I) = RSAV(LENRV1+I) + DO I = 1, LENRV1 + RVOD1(I) = RSAV(I) + END DO + DO I = 1, LENRV2 + RVOD2(I) = RSAV(LENRV1+I) + END DO IOFF = LENRV1 + LENRV2 - DO 114 I = 1, LRVK1 - 114 RVPK1(I) = RSAV(IOFF+I) -C - DO 120 I = 1, LENIV1 - 120 IVOD1(I) = ISAV(I) - DO 122 I = 1, LENIV2 - 122 IVOD2(I) = ISAV(LENIV1+I) + DO I = 1, LRVK1 + RVPK1(I) = RSAV(IOFF+I) + END DO +C + DO I = 1, LENIV1 + IVOD1(I) = ISAV(I) + END DO + DO I = 1, LENIV2 + IVOD2(I) = ISAV(LENIV1+I) + END DO IOFF = LENIV1 + LENIV2 - DO 124 I = 1, LIVK1 - 124 IVPK1(I) = ISAV(IOFF+I) + DO I = 1, LIVK1 + IVPK1(I) = ISAV(IOFF+I) + END DO C RETURN C----------------------- End of Subroutine VKSRC ----------------------- @@ -3209,16 +3246,18 @@ C listed (local) variables to be saved between calls to this integrator. 50 CONTINUE C Estimate the second derivative as a difference quotient in f. -------- T1 = T0 + HG - DO 60 I = 1, N - 60 Y(I) = Y0(I) + HG*YDOT(I) + DO I = 1, N + Y(I) = Y0(I) + HG*YDOT(I) + END DO IFAIL = 0 CALL F (N, T1, Y, TEMP, RPAR, IPAR, IFAIL) IF (IFAIL .NE. 0) THEN HNEW = PT2*HG GO TO 75 ENDIF - DO 70 I = 1, N - 70 TEMP(I) = (TEMP(I) - YDOT(I))/HG + DO I = 1, N + TEMP(I) = (TEMP(I) - YDOT(I))/HG + END DO YDDNRM = VNORML (N, TEMP, EWT) C Get the corresponding new value of h. -------------------------------- IF (YDDNRM*HUB*HUB .GT. TWO) THEN @@ -3347,11 +3386,13 @@ C listed (local) variables to be saved between calls to this integrator. IC = 1 IF (K .EQ. 0) GO TO 15 JJ1 = L - K - DO 10 JJ = JJ1, NQ - 10 IC = IC*JJ + DO JJ = JJ1, NQ + IC = IC*JJ + END DO 15 C = float(IC) - DO 20 I = 1, N - 20 DKY(I) = C*YH(I,L) + DO I = 1, N + DKY(I) = C*YH(I,L) + END DO IF (K .EQ. NQ) GO TO 55 JB2 = NQ - K DO 50 JB = 1, JB2 @@ -3360,11 +3401,13 @@ C listed (local) variables to be saved between calls to this integrator. IC = 1 IF (K .EQ. 0) GO TO 35 JJ1 = JP1 - K - DO 30 JJ = JJ1, J - 30 IC = IC*JJ + DO JJ = JJ1, J + IC = IC*JJ + END DO 35 C = float(IC) - DO 40 I = 1, N - 40 DKY(I) = C*YH(I,JP1) + S*DKY(I) + DO I = 1, N + DKY(I) = C*YH(I,JP1) + S*DKY(I) + END DO 50 CONTINUE IF (K .EQ. 0) RETURN 55 R = H**(-K) @@ -3583,8 +3626,9 @@ CALL VJUST (YH, LDYH, 1) I1 = 1 + (NEWQ + 1)*LDYH I2 = (MAXORD + 1)*LDYH IF (I1 .GT. I2) GO TO 120 - DO 110 I = I1, I2 - 110 YH1(I) = ZERO + DO I = I1, I2 + YH1(I) = ZERO + END DO 120 IF (NEWQ .LE. MAXORD) GO TO 140 FLOTL = float(LMAX) IF (MAXORD .LT. NQ-1) THEN @@ -3612,10 +3656,10 @@ CALL VJUST (YH, LDYH, -1) IF (NEWQ .LE. MAXORD) GO TO 50 C Rescale the history array for a change in H by a factor of ETA. ------ 150 R = ONE - DO 180 J = 2, L + DO J = 2, L R = R*ETA CALL SSCAL (N, R, YH(1,J), 1 ) - 180 CONTINUE + END DO H = HSCAL*ETA HSCAL = H RC = RC*ETA @@ -3628,11 +3672,12 @@ CALL SSCAL (N, R, YH(1,J), 1 ) C----------------------------------------------------------------------- 200 TN = TN + H I1 = NQNYH + 1 - DO 220 JB = 1, NQ + DO JB = 1, NQ I1 = I1 - LDYH - DO 210 I = I1, NQNYH - 210 YH1(I) = YH1(I) + YH1(I+LDYH) - 220 CONTINUE + DO I = I1, NQNYH + YH1(I) = YH1(I) + YH1(I+LDYH) + END DO + END DO CALL VSET RL1 = ONE/EL(2) RC = RC*(RL1/PRL1) @@ -3655,11 +3700,12 @@ C The VNLS routine failed to achieve convergence (NFLAG .NE. 0). ETAMAX = ONE TN = TOLD I1 = NQNYH + 1 - DO 430 JB = 1, NQ + DO JB = 1, NQ I1 = I1 - LDYH - DO 420 I = I1, NQNYH - 420 YH1(I) = YH1(I) - YH1(I+LDYH) - 430 CONTINUE + DO I = I1, NQNYH + YH1(I) = YH1(I) - YH1(I+LDYH) + END DO + END DO IF (NFLAG .LT. -1) GO TO 680 IF (ABS(H) .LE. HMIN*ONEPSM) GO TO 670 IF (NCF .EQ. MXNCF) GO TO 670 @@ -3721,9 +3767,10 @@ ccTDR call xerrab("") NST = NST + 1 HU = H NQU = NQ - DO 470 IBACK = 1, NQ + DO IBACK = 1, NQ I = L - IBACK - 470 TAU(I+1) = TAU(I) + TAU(I+1) = TAU(I) + END DO TAU(1) = H DO 480 J = 1, L CALL SAXPY (N, EL(J), ACOR, 1, YH(1,J), 1 ) @@ -3751,11 +3798,12 @@ CALL SCOPY (N, ACOR, 1, YH(1,LMAX), 1 ) NFLAG = -2 TN = TOLD I1 = NQNYH + 1 - DO 520 JB = 1, NQ + DO JB = 1, NQ I1 = I1 - LDYH - DO 510 I = I1, NQNYH - 510 YH1(I) = YH1(I) - YH1(I+LDYH) - 520 CONTINUE + DO I = I1, NQNYH + YH1(I) = YH1(I) - YH1(I+LDYH) + END DO + END DO IF (ABS(H) .LE. HMIN*ONEPSM) GO TO 660 ETAMAX = ONE IF (KFLAG .LE. KFC) GO TO 530 @@ -3789,8 +3837,9 @@ CALL VJUST (YH, LDYH, -1) CALL F (N, TN, Y, SAVF, RPAR, IPAR, IFAIL) NFE = NFE + 1 IF (IFAIL .NE. 0) GO TO 670 - DO 550 I = 1, N - 550 YH(I,2) = H*SAVF(I) + DO I = 1, N + YH(I,2) = H*SAVF(I) + END DO NQWAIT = 10 GO TO 200 C----------------------------------------------------------------------- @@ -3817,8 +3866,9 @@ CALL F (N, TN, Y, SAVF, RPAR, IPAR, IFAIL) IF (L .EQ. LMAX) GO TO 580 C Compute ratio of new H to current H at current order plus one. ------- CNQUOT = (TQ(5)/CONP)*(H/TAU(2))**L - DO 575 I = 1, N - 575 SAVF(I) = ACOR(I) - CNQUOT*YH(I,LMAX) + DO I = 1, N + SAVF(I) = ACOR(I) - CNQUOT*YH(I,LMAX) + END DO DUP = VNORML (N, SAVF, EWT)/TQ(3) ETAQP1 = ONE/((BIAS3*DUP)**(ONE/(FLOTL + ONE)) + ADDON) 580 IF (ETAQ .GE. ETAQP1) GO TO 590 @@ -3968,57 +4018,65 @@ GO TO (100, 200), METH 110 HSUM = H EM(1) = ONE FLOTNQ = FLOTL - ONE - DO 115 I = 2, L - 115 EM(I) = ZERO + DO I = 2, L + EM(I) = ZERO + END DO DO 150 J = 1, NQM1 IF ((J .NE. NQM1) .OR. (NQWAIT .NE. 1)) GO TO 130 S = ONE CSUM = ZERO - DO 120 I = 1, NQM1 + DO I = 1, NQM1 CSUM = CSUM + S*EM(I)/float(I+1) - 120 S = -S + S = -S + END DO TQ(1) = EM(NQM1)/(FLOTNQ*CSUM) 130 RXI = H/HSUM - DO 140 IBACK = 1, J + DO IBACK = 1, J I = (J + 2) - IBACK - 140 EM(I) = EM(I) + EM(I-1)*RXI + EM(I) = EM(I) + EM(I-1)*RXI + END DO HSUM = HSUM + TAU(J) 150 CONTINUE C Compute integral from -1 to 0 of polynomial and of x times it. ------- S = ONE EM0 = ZERO CSUM = ZERO - DO 160 I = 1, NQ + DO I = 1, NQ FLOTI = float(I) EM0 = EM0 + S*EM(I)/FLOTI CSUM = CSUM + S*EM(I)/(FLOTI+ONE) - 160 S = -S + S = -S + END DO C In EL, form coefficients of normalized integrated polynomial. -------- S = ONE/EM0 EL(1) = ONE - DO 170 I = 1, NQ - 170 EL(I+1) = S*EM(I)/float(I) + DO I = 1, NQ + EL(I+1) = S*EM(I)/float(I) + END DO XI = HSUM/H TQ(2) = XI*EM0/CSUM TQ(5) = XI/EL(L) IF (NQWAIT .NE. 1) GO TO 300 C For higher order control constant, multiply polynomial by 1+x/xi(q). - RXI = ONE/XI - DO 180 IBACK = 1, NQ + DO IBACK = 1, NQ I = (L + 1) - IBACK - 180 EM(I) = EM(I) + EM(I-1)*RXI + EM(I) = EM(I) + EM(I-1)*RXI + END DO C Compute integral of polynomial. -------------------------------------- S = ONE CSUM = ZERO - DO 190 I = 1, L + DO I = 1, L CSUM = CSUM + S*EM(I)/float(I+1) - 190 S = -S + S = -S + END DO TQ(3) = FLOTL*EM0/CSUM GO TO 300 C C Set coefficients for BDF methods. ------------------------------------ - 200 DO 210 I = 3, L - 210 EL(I) = ZERO + 200 DO I = 3, L + EL(I) = ZERO + END DO EL(1) = ONE EL(2) = ONE ALPH0 = -ONE @@ -4033,18 +4091,20 @@ GO TO (100, 200), METH RXI = H/HSUM JP1 = J + 1 ALPH0 = ALPH0 - ONE/float(JP1) - DO 220 IBACK = 1, JP1 + DO IBACK = 1, JP1 I = (J + 3) - IBACK - 220 EL(I) = EL(I) + EL(I-1)*RXI + EL(I) = EL(I) + EL(I-1)*RXI + END DO 230 CONTINUE ALPH0 = ALPH0 - ONE/float(NQ) RXIS = -EL(2) - ALPH0 HSUM = HSUM + TAU(NQM1) RXI = H/HSUM AHATN0 = -EL(2) - RXI - DO 235 IBACK = 1, NQ + DO IBACK = 1, NQ I = (NQ + 2) - IBACK - 235 EL(I) = EL(I) + EL(I-1)*RXIS + EL(I) = EL(I) + EL(I-1)*RXIS + END DO 240 T1 = ONE - AHATN0 + ALPH0 T2 = ONE + float(NQ)*T1 TQ(2) = ABS(ALPH0*T2/T1) @@ -4132,8 +4192,9 @@ GO TO (100, 200), METH 100 CONTINUE IF (IORD .EQ. 1) GO TO 180 C Order decrease. ------------------------------------------------------ - DO 110 J = 1, LMAX - 110 EL(J) = ZERO + DO J = 1, LMAX + EL(J) = ZERO + END DO EL(2) = ONE HSUM = ZERO DO 130 J = 1, NQM2 @@ -4141,25 +4202,29 @@ GO TO (100, 200), METH HSUM = HSUM + TAU(J) XI = HSUM/HSCAL JP1 = J + 1 - DO 120 IBACK = 1, JP1 + DO IBACK = 1, JP1 I = (J + 3) - IBACK - 120 EL(I) = EL(I)*XI + EL(I-1) + EL(I) = EL(I)*XI + EL(I-1) + END DO 130 CONTINUE C Construct coefficients of integrated polynomial. --------------------- - DO 140 J = 2, NQM1 - 140 EL(J+1) = float(NQ)*EL(J)/float(J) + DO J = 2, NQM1 + EL(J+1) = float(NQ)*EL(J)/float(J) + END DO C Subtract correction terms from YH array. ----------------------------- DO 170 J = 3, NQ - DO 160 I = 1, N - 160 YH(I,J) = YH(I,J) - YH(I,L)*EL(J) + DO I = 1, N + YH(I,J) = YH(I,J) - YH(I,L)*EL(J) + END DO 170 CONTINUE RETURN C Order increase. ------------------------------------------------------ C Zero out next column in YH array. ------------------------------------ 180 CONTINUE LP1 = L + 1 - DO 190 I = 1, N - 190 YH(I,LP1) = ZERO + DO I = 1, N + YH(I,LP1) = ZERO + END DO RETURN C----------------------------------------------------------------------- C Stiff option... @@ -4168,8 +4233,9 @@ GO TO (100, 200), METH 200 CONTINUE IF (IORD .EQ. 1) GO TO 300 C Order decrease. ------------------------------------------------------ - DO 210 J = 1, LMAX - 210 EL(J) = ZERO + DO J = 1, LMAX + EL(J) = ZERO + END DO EL(3) = ONE HSUM = ZERO DO 230 J = 1, NQM2 @@ -4177,19 +4243,22 @@ GO TO (100, 200), METH HSUM = HSUM + TAU(J) XI = HSUM/HSCAL JP1 = J + 1 - DO 220 IBACK = 1, JP1 + DO IBACK = 1, JP1 I = (J + 4) - IBACK - 220 EL(I) = EL(I)*XI + EL(I-1) + EL(I) = EL(I)*XI + EL(I-1) + END DO 230 CONTINUE C Subtract correction terms from YH array. ----------------------------- DO 250 J = 3, NQ - DO 240 I = 1, N - 240 YH(I,J) = YH(I,J) - YH(I,L)*EL(J) + DO I = 1, N + YH(I,J) = YH(I,J) - YH(I,L)*EL(J) + END DO 250 CONTINUE RETURN C Order increase. ------------------------------------------------------ - 300 DO 310 J = 1, LMAX - 310 EL(J) = ZERO + 300 DO J = 1, LMAX + EL(J) = ZERO + END DO EL(3) = ONE ALPH0 = -ONE ALPH1 = ONE @@ -4205,17 +4274,19 @@ GO TO (100, 200), METH PROD = PROD*XI ALPH0 = ALPH0 - ONE/float(JP1) ALPH1 = ALPH1 + ONE/XI - DO 320 IBACK = 1, JP1 + DO IBACK = 1, JP1 I = (J + 4) - IBACK - 320 EL(I) = EL(I)*XIOLD + EL(I-1) + EL(I) = EL(I)*XIOLD + EL(I-1) + END DO XIOLD = XI 330 CONTINUE 340 CONTINUE T1 = (-ALPH0 - ALPH1)/PROD C Load column L + 1 in YH array. --------------------------------------- LP1 = L + 1 - DO 350 I = 1, N - 350 YH(I,LP1) = T1*YH(I,LMAX) + DO I = 1, N + YH(I,LP1) = T1*YH(I,LMAX) + END DO C Add correction terms to YH array. ------------------------------------ NQP1 = NQ + 1 DO 370 J = 3, NQP1 From e0997b99586f212628909ec9fd9924771c4c0720 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 13:46:54 -0700 Subject: [PATCH 05/13] Fixes Shared DO termination label warnings in inelrates --- api/inelrates.m | 22 ++++++++++++---------- 1 file changed, 12 insertions(+), 10 deletions(-) diff --git a/api/inelrates.m b/api/inelrates.m index fa538b3e..ccc6ef0a 100755 --- a/api/inelrates.m +++ b/api/inelrates.m @@ -51,23 +51,25 @@ call xerrab("") read(us,inrad) close (us) c - do 102 i = 1,ncaset + do i = 1,ncaset terad(i)=terad(i)*1.602e-19/ctemp - 102 continue + end do c - do 103 i = 1,ncasent - rntau(i)=rntau(i)*crni*rl/cs - 103 continue + do i = 1,ncasent + rntau(i)=rntau(i)*crni*rl/cs + end do c c - do 111 i = 1,ncaset - do 111 j = 1,ncaseno - do 111 k = 1,ncasent + do i = 1,ncaset + do j = 1,ncaseno + do k = 1,ncasent c - radrate(i,j,k)=radrate(i,j,k)*crni*rl*1.602e-19 + radrate(i,j,k)=radrate(i,j,k)*crni*rl*1.602e-19 . /cs/ctemp/1.5 c - 111 continue + end do + end do + end do c c----------------------------------------------------------------------c elseif (impflag .eq. 2) then From 2259a1248e325c5f43fd895b469f17c71f8c2cb1 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 13:56:35 -0700 Subject: [PATCH 06/13] Removes DO termination statement warning from inelrates --- api/inelrates.m | 27 +++++++++++++++------------ 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/api/inelrates.m b/api/inelrates.m index ccc6ef0a..da9a8200 100755 --- a/api/inelrates.m +++ b/api/inelrates.m @@ -107,19 +107,21 @@ call xerrab("") c c ionization rate c - do 200 i = 1,ntev + do i = 1,ntev read(us,1000) tevb(i),(rsi(i,k),k=0,nz-1) - do 201 k = 0,nz-1 - 201 rsi(i,k) = rsi(i,k)*crni*rl/cs - 200 continue + do k = 0,nz-1 + rsi(i,k) = rsi(i,k)*crni*rl/cs + end do + end do c c recombination rate c - do 210 i = 1,ntev + do i = 1,ntev read(us,1000) tevb(i),(rre(i,k),k=1,nz) - do 211 k = 1,nz - 211 rre(i,k) = rre(i,k)*crni*rl/cs - 210 continue + do k = 1,nz + rre(i,k) = rre(i,k)*crni*rl/cs + end do + end do c c radiative power rate c @@ -134,11 +136,12 @@ c write(*,*) i,(rpwr(i,k),k=0,nz) c c CX recomb. rate c - do 230 i = 1,ntev + do i = 1,ntev read(us,1000) tevb(i),(rrcx(i,k),k=1,nz) - do 231 k = 1,nz - 231 rrcx(i,k) = rrcx(i,k)*crni*rl/cs - 230 continue + do k = 1,nz + rrcx(i,k) = rrcx(i,k)*crni*rl/cs + end do + end do c c ... Scale temperature scale for multi-charge rates to code units. do i = 1, ntev From 3de8a333804fa84a1e0e46aea3b84d4064580313 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 15:32:07 -0700 Subject: [PATCH 07/13] Remove DO termination statement warning from svrut3 --- svr/svrut3.m | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/svr/svrut3.m b/svr/svrut3.m index 5321083a..468157ac 100755 --- a/svr/svrut3.m +++ b/svr/svrut3.m @@ -356,20 +356,24 @@ real a(nn,ndiag),am(0:nn,ndiagm), c c Find the column in A where the main diagonal is stored c - do 1 k=1,ndiag - 1 if(idiag(k).eq.0) mdiag=k + do k=1,ndiag + if(idiag(k).eq.0) mdiag=k + end do c iprcnd=0 if(iprcnd.eq.1) then c c Use the identity matrix as the pre-conditioner c - do 111 i=1,nn - do 111 k=1,ndiag - 111 am(i,k) = 0.e0 + do i=1,nn + do k=1,ndiag + am(i,k) = 0.e0 + end do + end do c - do 2 i=1,nn - 2 am(i,mdiag) = 1.e0 + do i=1,nn + am(i,mdiag) = 1.e0 + end do c else c @@ -708,8 +712,9 @@ integer idiag(ndiagm) c c Find column in am(i,j) where main diagonal is located c - do 10 k=1,ndiag - 10 if(idiag(k).eq.0) mdiag=k + do k=1,ndiag + if(idiag(k).eq.0) mdiag=k + end do c c Compute the b' vector which is the solution to L(b')=y c From ac7fcb76cd152ece072b5dbe18b09e2958317ec7 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 15:33:19 -0700 Subject: [PATCH 08/13] Removes DO termination statement warning from nksol --- svr/nksol.m | 175 ++++++++++++++++++++++++++++++---------------------- 1 file changed, 101 insertions(+), 74 deletions(-) diff --git a/svr/nksol.m b/svr/nksol.m index d1c04115..d4798677 100755 --- a/svr/nksol.m +++ b/svr/nksol.m @@ -1412,8 +1412,9 @@ real function vnormnk (n, v, s) real v, s, sum dimension v(n), s(n) sum = 0.0e0 - do 10 i = 1,n - 10 sum = sum + (v(i)*s(i))**2 + do i = 1,n + sum = sum + (v(i)*s(i))**2 + end do vnormnk = sqrt(sum) return c----------------------- end of function vnormnk ------------------------- @@ -1510,8 +1511,9 @@ call pset (n, u, savf, su, sf, x, f, wm(locwmp), iwm(locimp), c----------------------------------------------------------------------- c load x with -f(u). c----------------------------------------------------------------------- - do 100 i = 1,n - 100 x(i) = -savf(i) + do i = 1,n + x(i) = -savf(i) + end do c----------------------------------------------------------------------- c call solpk to solve j*x = -f using the appropriate krylov c algorithm. @@ -1616,8 +1618,9 @@ go to (100,200,300), methk ib = iv + n*mmax ihes = ib + n iwk = ihes + mmax*mmax - do 110 i = 1,n - 110 wm(i+ib-1) = x(i) + do i = 1,n + wm(i+ib-1) = x(i) + end do call spiom (n, u, savf, wm(ib), su, sf, mmax, kmp, eps, * f, jac, psol, npsl, x, wm(iv), wm(ihes), iwm, miom, * wm(locwmp), iwm(locimp), wm(iwk), ipflg, iflag, rhom) @@ -1639,8 +1642,9 @@ call spiom (n, u, savf, wm(ib), su, sf, mmax, kmp, eps, ihsv = ihes + mmax*(mmaxp1+1) + 1 iwk = ihsv + mmax*mmaxp1 iq = iwk + n - do 210 i = 1,n - 210 wm(i+ib-1) = x(i) + do i = 1,n + wm(i+ib-1) = x(i) + end do call spigmr (n, u, savf, wm(ib), su, sf, mmax, mmaxp1, kmp, * eps, f, jac, psol, npsl, x, wm(iv), wm(ihes), wm(iq), * wm(ihsv), mgmr, wm(locwmp), iwm(locimp), wm(iwk), methn, @@ -1670,8 +1674,9 @@ call scopy (mmax, wm(ib), 1, wm(ihes), 1) ihsv = ihes + mmax*(mmaxp1+1) + 1 iwk = ihsv + mmax*mmaxp1 iq = iwk + n - do 310 i = 1,n - 310 wm(i+ib-1) = x(i) + do i = 1,n + wm(i+ib-1) = x(i) + end do ccc - this is direct solver do i=1,n @@ -1803,16 +1808,18 @@ c orthogonal vectors v(*,1) to v(*,miom). miom = 0 npsl = 0 c zero out the hes array. ---------------------------------------------- - do 15 j = 1,mmax - do 10 i = 1,mmax - 10 hes(i,j) = 0.0e0 - 15 continue + do j = 1,mmax + do i = 1,mmax + hes(i,j) = 0.0e0 + end do + end do c----------------------------------------------------------------------- c the initial residual is the vector b. apply scaling to b, and c then normalize as v(*,1). c----------------------------------------------------------------------- - do 20 i = 1,n - 20 v(i,1) = b(i)*sf(i) + do i = 1,n + v(i,1) = b(i)*sf(i) + end do bnrm = snrm2(n, v, 1) tem = 1.0e0/bnrm call sscal (n, tem, v(1,1), 1) @@ -1880,17 +1887,20 @@ c compute the approximation x(m) to the solution. c----------------------------------------------------------------------- 200 continue m = miom - do 210 k = 1,m - 210 b(k) = 0.0e0 + do k = 1,m + b(k) = 0.0e0 + end do b(1) = bnrm call shesl(hes,mmax,m,ipvt,b) - do 220 k = 1,n - 220 x(k) = 0.0e0 - do 230 i = 1,m + do k = 1,n + x(k) = 0.0e0 + end do + do i = 1,m call saxpy (n, b(i), v(1,i), 1, x, 1) - 230 continue - do 240 i = 1,n - 240 x(i) = x(i)/su(i) + end do + do i = 1,n + x(i) = x(i)/su(i) + end do if (ipflg .eq. 1) then ier = 0 call psol (n, u, savf, su, sf, f, jac, wk, wmp, iwmp, x, ier) @@ -1978,8 +1988,9 @@ c f(u) = 0. c----------------------------------------------------------------------- c mfdif = 1. call user-supplied routine for computing j*v. c----------------------------------------------------------------------- - do 110 i = 1,n - 110 vtem(i) = v(i)/su(i) + do i = 1,n + vtem(i) = v(i)/su(i) + end do if (ipflg .eq. 1) then ier = 0 call psol (n, u, savf, su, sf, f, jac, z, wmp, iwmp, vtem, ier) @@ -1989,9 +2000,9 @@ call psol (n, u, savf, su, sf, f, jac, z, wmp, iwmp, vtem, ier) call jac (n, u, savf, vtem, z, wmp, iwmp) nje = nje + 1 c ... Add a 1/(pseudo time-step) carried in u(n+2): TDR 8/21/95 - do 118 i = 1,n + do i = 1,n z(i) = z(i) + vtem(i)*sf(i)*u(n+2) # sf in numerator or dem?? - 118 continue + end do c 3-15-93. this scaling should already be performed in the jac routine. c do 120 i = 1,n c 120 z(i) = z(i)*sf(i) @@ -2001,24 +2012,28 @@ call jac (n, u, savf, vtem, z, wmp, iwmp) c----------------------------------------------------------------------- 200 continue c set vtem = (su-inverse) * v. - do 210 i = 1,n + do i = 1,n vtem(i) = v(i)/su(i) - 210 continue + end do if (ipflg .eq. 0) then c ipflg = 0. save u in z and increment u by sigma*vtem. - do 215 i = 1,n - 215 z(i) = u(i)*su(i) + do i = 1,n + z(i) = u(i)*su(i) + end do utv = sdot(n, z, 1, v, 1) suitv = zero - do 220 i = 1,n - 220 suitv = suitv + abs(v(i)) + do i = 1,n + suitv = suitv + abs(v(i)) + end do c vtv = 1 in this case. sigma = sqteta*max(abs(utv),suitv) sigma = sign(sigma,utv) - do 230 i = 1,n - 230 z(i) = u(i) - do 240 i = 1,n - 240 u(i) = z(i) + sigma*vtem(i) + do i = 1,n + z(i) = u(i) + end do + do i = 1,n + u(i) = z(i) + sigma*vtem(i) + end do else c ipflg = 1. apply inverse of right preconditioner ier = 0 @@ -2028,18 +2043,21 @@ call psol (n, u, savf, su, sf, f, jac, ftem, wmp, iwmp, if (ier .ne. 0) return vtv = zero suitv = zero - do 260 i = 1,n + do i = 1,n tmp = vtem(i)*su(i) z(i) = tmp*su(i) vtv = vtv + tmp*tmp - 260 suitv = suitv + abs(tmp) + suitv = suitv + abs(tmp) + end do utv = sdot(n, u, 1, z, 1) sigma = sqteta*max(abs(utv),suitv)/vtv sigma = sign(sigma,utv) - do 270 i = 1,n - 270 z(i) = u(i) - do 280 i = 1,n - 280 u(i) = z(i) + sigma*vtem(i) + do i = 1,n + z(i) = u(i) + end do + do i = 1,n + u(i) = z(i) + sigma*vtem(i) + end do endif call f(n, u, ftem) nfe = nfe + 1 @@ -2048,15 +2066,16 @@ ccc u(i) = z(i) - sigma*vtem(i) # change sign of perturbation ccc 281 continue ccc call f(n, u, vtem) # here vtem is second ftem ccc nfe = nfe + 1 # End 2nd order Jac - do 290 i = 1,n - 290 u(i) = z(i) + do i = 1,n + u(i) = z(i) + end do c ... u(n+2) contains the 1/(pseudo time-step) nufak: TDR 8/20/95 - do 300 i = 1,n - 300 z(i) = (ftem(i) - savf(i))/sigma - vtem(i)*u(n+2) -ccc 300 z(i) = (ftem(i) - vtem(i))/(2*sigma) - vtem(i)*u(n+2) # 2nd order Jac -cccc vtem(i)*su(i)*u(n+2)/sf(i) - do 310 i = 1,n - 310 z(i) = z(i)*sf(i) + do i = 1,n + z(i) = (ftem(i) - savf(i))/sigma - vtem(i)*u(n+2) + end do + do i = 1,n + z(i) = z(i)*sf(i) + end do return c----------------------- end of subroutine atv ------------------------- end @@ -2438,18 +2457,19 @@ dimension savf(n), x(n), u(*), wmp(*), wk(n), iwmp(*) mgmr = 0 npsl = 0 c zero out the hes and hessv arrays ------------------------------------ - do 15 j = 1,mmax - do 10 i = 1,mmaxp1 + do j = 1,mmax + do i = 1,mmaxp1 hes(i,j) = 0.0e0 hessv(i,j) = 0.0e0 - 10 continue - 15 continue + end do + end do c----------------------------------------------------------------------- c the initial residual is the vector b. apply scaling to b, and c then normalize as v(*,1). c----------------------------------------------------------------------- - do 20 i = 1,n - 20 v(i,1) = b(i)*sf(i) + do i = 1,n + v(i,1) = b(i)*sf(i) + end do bnrm = snrm2(n, v, 1) tem = 1.0e0/bnrm call sscal (n, tem, v(1,1), 1) @@ -2510,18 +2530,21 @@ call sscal (n, tem, v(1,m+1), 1) 200 continue m = mgmr mp1 = m + 1 - do 210 k = 1,mp1 - 210 b(k) = 0.0e0 + do k = 1,mp1 + b(k) = 0.0e0 + end do b(1) = bnrm call shels(hes,mmaxp1,m,q,b) if ((methn .eq. 0) .or. (methn .eq. 2)) then - do 220 k = 1,n - 220 x(k) = 0.0e0 - do 230 i = 1,m + do k = 1,n + x(k) = 0.0e0 + end do + do i = 1,m call saxpy(n,b(i),v(1,i),1,x,1) - 230 continue - do 240 i = 1,n - 240 x(i) = x(i)/su(i) + end do + do i = 1,n + x(i) = x(i)/su(i) + end do if (ipflg .eq. 1) then ier = 0 call psol (n, u, savf, su, sf, f, jac, wk, wmp, iwmp, x, ier) @@ -2889,8 +2912,9 @@ call sscal(m, t, ycp, 1) do 100 i = 1,m call saxpy (n, ynew(i), v(1,i), 1, xnew, 1) 100 continue - do 110 i = 1,n - 110 xnew(i) = xnew(i)/su(i) + do i = 1,n + xnew(i) = xnew(i)/su(i) + end do if (ipflg .ge. 1) then ier = 0 call psol (n, u, savf, su, sf, f, jac, wk, wmp, iwmp, xnew, ier) @@ -3306,8 +3330,9 @@ call sswap(n, u, 1, unew, 1) f1nprv = f1nrmp rl = min(two*rl,rlmax) call scopy(n, u, 1, unew, 1) - do 140 i = 1,n - 140 u(i) = unew(i) + rl*p(i) + do i = 1,n + u(i) = unew(i) + rl*p(i) + end do call f(n, u, savf) nfe = nfe + 1 call sswap(n, u, 1, unew, 1) @@ -3335,8 +3360,9 @@ call sswap(n, u, 1, unew, 1) rlincr = rldiff/two rl = rllo + rlincr call scopy(n, u, 1, unew, 1) - do 160 i = 1,n - 160 u(i) = unew(i) + rl*p(i) + do i = 1,n + u(i) = unew(i) + rl*p(i) + end do call f(n, u, savf) nfe = nfe + 1 call sswap(n, u, 1, unew, 1) @@ -3365,8 +3391,9 @@ call sswap(n, u, 1, unew, 1) c u value that satisfied alpha-condition and continue. c increment counter on number of steps beta-condition not satisfied. call scopy(n, u, 1, unew, 1) - do 170 i = 1,n - 170 u(i) = unew(i) + rllo*p(i) + do i = 1,n + u(i) = unew(i) + rllo*p(i) + end do call f(n, u, savf) nfe = nfe + 1 call sswap(n, u, 1, unew, 1) From 0a4b4ab42acd4adfbd6761c5a4664ee5b579d2c8 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 15:58:14 -0700 Subject: [PATCH 09/13] Removes DO termination statement warning from daspk.m --- svr/daspk.m | 360 ++++++++++++++++++++++++++++++---------------------- 1 file changed, 209 insertions(+), 151 deletions(-) diff --git a/svr/daspk.m b/svr/daspk.m index 687149f9..5967cd3a 100755 --- a/svr/daspk.m +++ b/svr/daspk.m @@ -1734,8 +1734,9 @@ CALL SDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) CALL SINVWT(NEQ,RWORK(LWT),IER) IF (IER .NE. 0) GO TO 713 IF (INFO(16) .NE. 0) THEN - DO 305 I = 1, NEQ - 305 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) + DO I = 1, NEQ + RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) + END DO ENDIF C C Compute unit roundoff and HMIN. @@ -1857,9 +1858,10 @@ CALL SDAWTS(NEQ,INFO(2),RTOL,ATOL,Y,RWORK(LWT),RPAR,IPAR) CALL SINVWT(NEQ,RWORK(LWT),IER) IF (IER .NE. 0) GO TO 713 IF (INFO(16) .NE. 0) THEN - DO 357 I = 1, NEQ - 357 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) - ENDIF + DO I = 1, NEQ + RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) + END DO + ENDIF C C Reset the initial stepsize to be used by SDSTP. C Use H0, if this was input. Otherwise, recompute H0, @@ -2053,8 +2055,9 @@ CALL SINVWT(NEQ,RWORK(LWT),IER) GO TO 527 ENDIF IF (INFO(16) .NE. 0) THEN - DO 515 I = 1, NEQ - 515 RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) + DO I = 1, NEQ + RWORK(LVT+I-1) = MAX(IWORK(LID+I-1),0)*RWORK(LWT+I-1) + END DO ENDIF C C Test for too much accuracy requested. @@ -2069,9 +2072,10 @@ CALL SINVWT(NEQ,RWORK(LWT),IER) ATOL(1)=R*ATOL(1) IDID=-2 GO TO 527 -523 DO 524 I=1,NEQ +523 DO I=1,NEQ RTOL(I)=R*RTOL(I) -524 ATOL(I)=R*ATOL(I) + ATOL(I)=R*ATOL(I) + END DO IDID=-2 GO TO 527 525 CONTINUE @@ -2806,7 +2810,7 @@ DIMENSION RPAR(*),IPAR(*), RDAT(*), IDAT(*) TEMP1=H GAMMA(1)=0.0E0 SIGMA(1)=1.0E0 - DO 210 I=2,KP1 + DO I=2,KP1 TEMP2=PSI(I-1) PSI(I-1)=TEMP1 BETA(I)=BETA(I-1)*PSI(I-1)/TEMP2 @@ -2814,7 +2818,7 @@ DIMENSION RPAR(*),IPAR(*), RDAT(*), IDAT(*) ALPHA(I)=H/TEMP1 SIGMA(I)=(I-1)*SIGMA(I-1)*ALPHA(I) GAMMA(I)=GAMMA(I-1)+ALPHA(I-1)/H -210 CONTINUE + END DO PSI(KP1)=TEMP1 230 CONTINUE C @@ -2822,10 +2826,10 @@ DIMENSION RPAR(*),IPAR(*), RDAT(*), IDAT(*) C ALPHAS = 0.0E0 ALPHA0 = 0.0E0 - DO 240 I = 1,K + DO I = 1,K ALPHAS = ALPHAS - 1.0E0/I ALPHA0 = ALPHA0 - ALPHA(I) -240 CONTINUE + END DO C C Compute leading coefficient CJ C @@ -2840,10 +2844,11 @@ DIMENSION RPAR(*),IPAR(*), RDAT(*), IDAT(*) C Change PHI to PHI STAR C IF(KP1 .LT. NSP1) GO TO 280 - DO 270 J=NSP1,KP1 - DO 260 I=1,NEQ -260 PHI(I,J)=BETA(J)*PHI(I,J) -270 CONTINUE + DO J=NSP1,KP1 + DO I=1,NEQ + PHI(I,J)=BETA(J)*PHI(I,J) + END DO + END DO 280 CONTINUE C C Update time @@ -2892,16 +2897,18 @@ CALL NLS(X,Y,YPRIME,NEQ, EST = ERK KNEW=K IF(K .EQ. 1)GO TO 430 - DO 405 I = 1,NEQ -405 DELTA(I) = PHI(I,KP1) + E(I) + DO I = 1,NEQ + DELTA(I) = PHI(I,KP1) + E(I) + END DO ERKM1=SIGMA(K)*SDWNRM(NEQ,DELTA,VT,RPAR,IPAR) TERKM1 = K*ERKM1 IF(K .GT. 2)GO TO 410 IF(TERKM1 .LE. 0.5*TERK)GO TO 420 GO TO 430 410 CONTINUE - DO 415 I = 1,NEQ -415 DELTA(I) = PHI(I,K) + DELTA(I) + DO I = 1,NEQ + DELTA(I) = PHI(I,K) + DELTA(I) + END DO ERKM2=SIGMA(K-1)*SDWNRM(NEQ,DELTA,VT,RPAR,IPAR) TERKM2 = (K-1)*ERKM2 IF(MAX(TERKM1,TERKM2).GT.TERK)GO TO 430 @@ -3027,15 +3034,19 @@ ccc call xerrab("") C 575 CONTINUE IF(KOLD.EQ.IWM(LMXORD))GO TO 585 - DO 580 I=1,NEQ -580 PHI(I,KP2)=E(I) + DO I=1,NEQ + PHI(I,KP2)=E(I) + END DO 585 CONTINUE - DO 590 I=1,NEQ -590 PHI(I,KP1)=PHI(I,KP1)+E(I) - DO 595 J1=2,KP1 + DO I=1,NEQ + PHI(I,KP1)=PHI(I,KP1)+E(I) + END DO + DO J1=2,KP1 J=KP1-J1+1 - DO 595 I=1,NEQ -595 PHI(I,J)=PHI(I,J)+PHI(I,J+1) + DO I=1,NEQ + PHI(I,J)=PHI(I,J)+PHI(I,J+1) + END DO + END DO JSTART = 1 RETURN C @@ -3057,14 +3068,16 @@ ccc call xerrab("") C X=XOLD IF(KP1.LT.NSP1)GO TO 630 - DO 620 J=NSP1,KP1 + DO J=NSP1,KP1 TEMP1=1.0E0/BETA(J) - DO 610 I=1,NEQ -610 PHI(I,J)=TEMP1*PHI(I,J) -620 CONTINUE + DO I=1,NEQ + PHI(I,J)=TEMP1*PHI(I,J) + END DO + END DO 630 CONTINUE - DO 640 I=2,KP1 -640 PSI(I-1)=PSI(I)-H + DO I=2,KP1 + PSI(I-1)=PSI(I)-H + END DO C C C Test whether failure is due to nonlinear solver @@ -3155,9 +3168,10 @@ CALL SDATRP(X,X,Y,YPRIME,NEQ,K,PHI,PSI) C 690 IF (KOLD .EQ. 0) THEN PSI(1) = H - DO 695 I = 1,NEQ -695 PHI(I,2) = R*PHI(I,2) - ENDIF + DO I = 1,NEQ + PHI(I,2) = R*PHI(I,2) + END DO + ENDIF GO TO 200 C C------END OF SUBROUTINE SDSTP------------------------------------------ @@ -3435,11 +3449,12 @@ IMPLICIT real(A-H,O-Z) INTEGER I,NEQ,IER DIMENSION WT(*) C - DO 10 I = 1,NEQ + DO I = 1,NEQ IF (WT(I) .LE. 0.0E0) GO TO 30 - 10 CONTINUE - DO 20 I = 1,NEQ - 20 WT(I) = 1.0E0/WT(I) + END DO + DO I = 1,NEQ + WT(I) = 1.0E0/WT(I) + END DO IER = 0 RETURN C @@ -3492,20 +3507,22 @@ DIMENSION YOUT(*),YPOUT(*) DIMENSION PHI(NEQ,*),PSI(*) KOLDP1=KOLD+1 TEMP1=XOUT-X - DO 10 I=1,NEQ + DO I=1,NEQ YOUT(I)=PHI(I,1) -10 YPOUT(I)=0.0E0 + YPOUT(I)=0.0E0 + END DO C=1.0E0 D=0.0E0 GAMMA=TEMP1/PSI(1) - DO 30 J=2,KOLDP1 + DO J=2,KOLDP1 D=D*GAMMA+C/PSI(J-1) C=C*GAMMA GAMMA=(TEMP1+PSI(J-1))/PSI(J) - DO 20 I=1,NEQ + DO I=1,NEQ YOUT(I)=YOUT(I)+C*PHI(I,J) -20 YPOUT(I)=YPOUT(I)+D*PHI(I,J) -30 CONTINUE + YPOUT(I)=YPOUT(I)+D*PHI(I,J) + END DO + END DO RETURN C C------END OF SUBROUTINE SDATRP----------------------------------------- @@ -3538,13 +3555,14 @@ DIMENSION V(*),RWT(*) DIMENSION RPAR(*),IPAR(*) SDWNRM = 0.0E0 VMAX = 0.0E0 - DO 10 I = 1,NEQ + DO I = 1,NEQ IF(ABS(V(I)*RWT(I)) .GT. VMAX) VMAX = ABS(V(I)*RWT(I)) -10 CONTINUE + END DO IF(VMAX .LE. 0.0E0) GO TO 30 SUM = 0.0E0 - DO 20 I = 1,NEQ -20 SUM = SUM + ((V(I)*RWT(I))/VMAX)**2 + DO I = 1,NEQ + SUM = SUM + ((V(I)*RWT(I))/VMAX)**2 + END DO SDWNRM = VMAX*SQRT(SUM/NEQ) 30 CONTINUE RETURN @@ -4002,8 +4020,9 @@ CALL SCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR) IVIO = 1 RATIO1 = TAU/PNRM RATIO = RATIO*RATIO1 - DO 20 I = 1,NEQ - 20 P(I) = P(I)*RATIO1 + DO I = 1,NEQ + P(I) = P(I)*RATIO1 + END DO PNRM = TAU IF (KPRIN .GE. 2) THEN MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)' @@ -4327,14 +4346,16 @@ DIMENSION PHI(NEQ,*),GAMMA(*) C Predict the solution and derivative and compute the tolerance C for the Newton iteration. C - DO 310 I=1,NEQ + DO I=1,NEQ Y(I)=PHI(I,1) -310 YPRIME(I)=0.0E0 - DO 330 J=2,KP1 - DO 320 I=1,NEQ + YPRIME(I)=0.0E0 + END DO + DO J=2,KP1 + DO I=1,NEQ Y(I)=Y(I)+PHI(I,J) -320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) -330 CONTINUE + YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) + END DO + END DO PNORM = SDWNRM (NEQ,Y,WT,RPAR,IPAR) TOLNEW = 100.E0*UROUND*PNORM C @@ -4383,12 +4404,14 @@ CALL SNSD(X,Y,YPRIME,NEQ,RES,PDUM,WT,RPAR,IPAR,DUMSVR, C large, then consider the corrector iteration to have failed. C 375 IF(NONNEG .EQ. 0) GO TO 390 - DO 377 I = 1,NEQ -377 DELTA(I) = MIN(Y(I),0.0E0) + DO I = 1,NEQ + DELTA(I) = MIN(Y(I),0.0E0) + END DO DELNRM = SDWNRM(NEQ,DELTA,WT,RPAR,IPAR) IF(DELNRM .GT. EPCON) GO TO 380 - DO 378 I = 1,NEQ -378 E(I) = E(I) - DELTA(I) + DO I = 1,NEQ + E(I) = E(I) - DELTA(I) + END DO GO TO 390 C C @@ -4512,8 +4535,9 @@ DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) C Initialize Newton counter M and accumulation vector E. C M = 0 - DO 100 I=1,NEQ -100 E(I)=0.0E0 + DO I=1,NEQ + E(I)=0.0E0 + END DO C C Corrector loop. C @@ -4523,8 +4547,9 @@ DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) C If necessary, multiply residual by convergence factor. C IF (MULDEL .EQ. 1) THEN - DO 320 I = 1,NEQ -320 DELTA(I) = DELTA(I) * CONFAC + DO I = 1,NEQ + DELTA(I) = DELTA(I) * CONFAC + END DO ENDIF C C Compute a new iterate (back-substitution). @@ -4534,10 +4559,11 @@ CALL SSLVD(NEQ,DELTA,WM,IWM) C C Update Y, E, and YPRIME. C - DO 340 I=1,NEQ + DO I=1,NEQ Y(I)=Y(I)-DELTA(I) E(I)=E(I)-DELTA(I) -340 YPRIME(I)=YPRIME(I)-CJ*DELTA(I) + YPRIME(I)=YPRIME(I)-CJ*DELTA(I) + END DO C C Test for convergence of the iteration. C @@ -4670,8 +4696,9 @@ GO TO (100,200,300,400,500),MTYPE C Dense user-supplied matrix. C 100 LENPD=IWM(LNPD) - DO 110 I=1,LENPD -110 WM(I)=0.0E0 + DO I=1,LENPD + WM(I)=0.0E0 + END DO CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) GO TO 230 C @@ -4681,7 +4708,7 @@ CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) 200 IRES=0 NROW=0 SQUR = SQRT(UROUND) - DO 210 I=1,NEQ + DO I=1,NEQ DEL=SQUR*MAX(ABS(Y(I)),ABS(H*YPRIME(I)), * ABS(1.E0/EWT(I))) DEL=SIGN(DEL,H*YPRIME(I)) @@ -4694,12 +4721,13 @@ CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR) IF (IRES .LT. 0) RETURN DELINV=1.0E0/DEL - DO 220 L=1,NEQ -220 WM(NROW+L)=(E(L)-DELTA(L))*DELINV + DO L=1,NEQ + WM(NROW+L)=(E(L)-DELTA(L))*DELINV + END DO NROW=NROW+NEQ Y(I)=YSAVE YPRIME(I)=YPSAVE -210 CONTINUE + END DO C C C Do dense-matrix LU decomposition on J. @@ -4716,8 +4744,9 @@ C Dummy section for IWM(MTYPE)=3. C Banded user-supplied matrix. C 400 LENPD=IWM(LNPD) - DO 410 I=1,LENPD -410 WM(I)=0.0E0 + DO I=1,LENPD + WM(I)=0.0E0 + END DO CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) MEBAND=2*IWM(LML)+IWM(LMU)+1 GO TO 550 @@ -4734,8 +4763,8 @@ CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) IPSAVE=ISAVE+MSAVE IRES=0 SQUR=SQRT(UROUND) - DO 540 J=1,MBA - DO 510 N=J,NEQ,MBAND + DO J=1,MBA + DO N=J,NEQ,MBAND K= (N-J)/MBAND + 1 WM(ISAVE+K)=Y(N) WM(IPSAVE+K)=YPRIME(N) @@ -4744,11 +4773,12 @@ CALL JACD(X,Y,YPRIME,WM,CJ,RPAR,IPAR) DEL=SIGN(DEL,H*YPRIME(N)) DEL=(Y(N)+DEL)-Y(N) Y(N)=Y(N)+DEL -510 YPRIME(N)=YPRIME(N)+CJ*DEL + YPRIME(N)=YPRIME(N)+CJ*DEL + END DO IWM(LNRE)=IWM(LNRE)+1 CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR) IF (IRES .LT. 0) RETURN - DO 530 N=J,NEQ,MBAND + DO N=J,NEQ,MBAND K= (N-J)/MBAND + 1 Y(N)=WM(ISAVE+K) YPRIME(N)=WM(IPSAVE+K) @@ -4760,10 +4790,11 @@ CALL RES(X,Y,YPRIME,CJ,E,IRES,RPAR,IPAR) I1=MAX(1,(N-IWM(LMU))) I2=MIN(NEQ,(N+IWM(LML))) II=N*MEB1-IWM(LML) - DO 520 I=I1,I2 -520 WM(II+I)=(E(I)-DELTA(I))*DELINV -530 CONTINUE -540 CONTINUE + DO I=I1,I2 + WM(II+I)=(E(I)-DELTA(I))*DELINV + END DO + END DO + END DO C C C Do LU decomposition of banded J. @@ -5338,8 +5369,9 @@ CALL SCNSTR (NEQ, Y, YNEW, ICNSTR, TAU, RLX, IRET, IVAR) IVIO = 1 RATIO1 = TAU/PNRM RATIO = RATIO*RATIO1 - DO 20 I = 1,NEQ - 20 P(I) = P(I)*RATIO1 + DO I = 1,NEQ + P(I) = P(I)*RATIO1 + END DO PNRM = TAU IF (KPRIN .GE. 2) THEN MSG = '------ CONSTRAINT VIOL., PNRM = (R1), INDEX = (I1)' @@ -5691,14 +5723,16 @@ DIMENSION GAMMA(*),RPAR(*),IPAR(*) C Predict the solution and derivative and compute the tolerance C for the Newton iteration. C - DO 310 I=1,NEQ + DO I=1,NEQ Y(I)=PHI(I,1) -310 YPRIME(I)=0.0E0 - DO 330 J=2,KP1 - DO 320 I=1,NEQ - Y(I)=Y(I)+PHI(I,J) -320 YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) -330 CONTINUE + YPRIME(I)=0.0E0 + END DO + DO J=2,KP1 + DO I=1,NEQ + Y(I)=Y(I)+PHI(I,J) + YPRIME(I)=YPRIME(I)+GAMMA(J)*PHI(I,J) + END DO + END DO EPLIN = EPLI*EPCON TOLNEW = EPLIN C @@ -5746,12 +5780,14 @@ CALL SNSK(X,Y,YPRIME,NEQ,RES,PSOL,WT,RPAR,IPAR,SAVR, C large, then consider the corrector iteration to have failed. C IF(NONNEG .EQ. 0) GO TO 390 - DO 360 I = 1,NEQ - 360 DELTA(I) = MIN(Y(I),0.0E0) + DO I = 1,NEQ + DELTA(I) = MIN(Y(I),0.0E0) + END DO DELNRM = SDWNRM(NEQ,DELTA,WT,RPAR,IPAR) IF(DELNRM .GT. EPCON) GO TO 380 - DO 370 I = 1,NEQ - 370 E(I) = E(I) - DELTA(I) + DO I = 1,NEQ + E(I) = E(I) - DELTA(I) + END DO GO TO 390 C C @@ -5885,8 +5921,9 @@ DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) C Initialize Newton counter M and accumulation vector E. C M = 0 - DO 100 I=1,NEQ -100 E(I) = 0.0E0 + DO I=1,NEQ + E(I) = 0.0E0 + END DO C C Corrector loop. C @@ -5896,14 +5933,16 @@ DIMENSION WM(*),IWM(*), RPAR(*),IPAR(*) C If necessary, multiply residual by convergence factor. C IF (MULDEL .EQ. 1) THEN - DO 320 I = 1,NEQ -320 DELTA(I) = DELTA(I) * CONFAC - ENDIF + DO I = 1,NEQ + DELTA(I) = DELTA(I) * CONFAC + END DO + ENDIF C C Save residual in SAVR. C - DO 340 I = 1,NEQ -340 SAVR(I) = DELTA(I) + DO I = 1,NEQ + SAVR(I) = DELTA(I) + END DO C C Compute a new iterate. Store the correction in DELTA. C @@ -5914,10 +5953,11 @@ CALL SSLVK (NEQ, Y, X, YPRIME, SAVR, DELTA, WT, WM, IWM, C C Update Y, E, and YPRIME. C - DO 360 I=1,NEQ + DO I=1,NEQ Y(I) = Y(I) - DELTA(I) E(I) = E(I) - DELTA(I) -360 YPRIME(I) = YPRIME(I) - CJ*DELTA(I) + YPRIME(I) = YPRIME(I) - CJ*DELTA(I) + END DO C C Test for convergence of the iteration. C @@ -6057,8 +6097,9 @@ DIMENSION Y(*), YPRIME(*), SAVR(*), X(*), EWT(*), LZ = LDL + NEQ CALL SSCAL (NEQ, RSQRTN, EWT, 1) CALL SCOPY (NEQ, X, 1, WM(LR), 1) - DO 110 I = 1,NEQ - 110 X(I) = 0.E0 + DO I = 1,NEQ + X(I) = 0.E0 + END DO C----------------------------------------------------------------------- C Top of loop for the restart algorithm. Initial pass approximates C X and sets up a transformed system to perform subsequent restarts @@ -6081,8 +6122,9 @@ CALL SSPIGM (NEQ, TN, Y, YPRIME, SAVR, WM(LR), EWT, MAXL, MAXLP1, NLI = NLI + LGMR NPS = NPS + NPSL NRE = NRE + NRES - DO 120 I = 1,NEQ - 120 X(I) = X(I) + WM(LZ+I-1) + DO I = 1,NEQ + X(I) = X(I) + WM(LZ+I-1) + END DO IF ((IFLAG .EQ. 1) .AND. (NRSTS .LT. NRMAX) .AND. (IRES .EQ. 0)) 1 GO TO 115 C----------------------------------------------------------------------- @@ -6259,8 +6301,9 @@ DIMENSION Y(*), YPRIME(*), SAVR(*), R(*), WGHT(*), Z(*), C The initial guess for Z is 0. The initial residual is therefore C the vector R. Initialize Z to 0. C----------------------------------------------------------------------- - DO 10 I = 1,NEQ - 10 Z(I) = 0.0E0 + DO I = 1,NEQ + Z(I) = 0.0E0 + END DO C----------------------------------------------------------------------- C Apply inverse of left preconditioner to vector R if NRSTS .EQ. 0. C Form V(*,1), the scaled preconditioned right hand side. @@ -6270,11 +6313,13 @@ CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, WK, CJ, WGHT, WP, IWP, 1 R, EPLIN, IER, RPAR, IPAR) NPSL = 1 IF (IER .NE. 0) GO TO 300 - DO 30 I = 1,NEQ - 30 V(I,1) = R(I)*WGHT(I) + DO I = 1,NEQ + V(I,1) = R(I)*WGHT(I) + END DO ELSE - DO 35 I = 1,NEQ - 35 V(I,1) = R(I) + DO I = 1,NEQ + V(I,1) = R(I) + END DO ENDIF C----------------------------------------------------------------------- C Calculate norm of scaled vector V(*,1) and normalize it @@ -6291,10 +6336,11 @@ CALL SSCAL (NEQ, TEM, V(1,1), 1) C----------------------------------------------------------------------- C Zero out the HES array. C----------------------------------------------------------------------- - DO 65 J = 1,MAXL - DO 60 I = 1,MAXLP1 - 60 HES(I,J) = 0.0E0 - 65 CONTINUE + DO J = 1,MAXL + DO I = 1,MAXLP1 + HES(I,J) = 0.0E0 + END DO + END DO C----------------------------------------------------------------------- C Main loop to compute the vectors V(*,2) to V(*,MAXL). C The running product PROD is needed for the convergence test. @@ -6328,20 +6374,22 @@ CALL SHEQR (HES, MAXLP1, LL, Q, INFO, LL) IF ((LL.GT.KMP) .AND. (KMP.LT.MAXL)) THEN IF (LL .EQ. KMP+1) THEN CALL SCOPY (NEQ, V(1,1), 1, DL, 1) - DO 75 I = 1,KMP + DO I = 1,KMP IP1 = I + 1 I2 = I*2 S = Q(I2) C = Q(I2-1) - DO 70 K = 1,NEQ - 70 DL(K) = S*DL(K) + C*V(K,IP1) - 75 CONTINUE + DO K = 1,NEQ + DL(K) = S*DL(K) + C*V(K,IP1) + END DO + END DO ENDIF S = Q(2*LL) C = Q(2*LL-1)/SNORMW LLP1 = LL + 1 - DO 80 K = 1,NEQ - 80 DL(K) = S*DL(K) + C*V(K,LLP1) + DO K = 1,NEQ + DL(K) = S*DL(K) + C*V(K,LLP1) + END DO DLNRM = SNRM2 (NEQ, DL, 1) RHO = RHO*DLNRM ENDIF @@ -6361,8 +6409,9 @@ CALL SSCAL (NEQ, TEM, V(1,LL+1), 1) IF (RHO .LT. RNRM) GO TO 150 120 CONTINUE IFLAG = 2 - DO 130 I = 1,NEQ - 130 Z(I) = 0.E0 + DO I = 1,NEQ + Z(I) = 0.E0 + END DO RETURN 150 IFLAG = 1 C----------------------------------------------------------------------- @@ -6378,18 +6427,20 @@ C Calculate DL from the V(I)'s. C CALL SCOPY (NEQ, V(1,1), 1, DL, 1) MAXLM1 = MAXL - 1 - DO 175 I = 1,MAXLM1 + DO I = 1,MAXLM1 IP1 = I + 1 I2 = I*2 S = Q(I2) C = Q(I2-1) - DO 170 K = 1,NEQ - 170 DL(K) = S*DL(K) + C*V(K,IP1) - 175 CONTINUE + DO K = 1,NEQ + DL(K) = S*DL(K) + C*V(K,IP1) + END DO + END DO S = Q(2*MAXL) C = Q(2*MAXL-1)/SNORMW - DO 180 K = 1,NEQ - 180 DL(K) = S*DL(K) + C*V(K,MAXLP1) + DO K = 1,NEQ + DL(K) = S*DL(K) + C*V(K,MAXLP1) + END DO ENDIF C C Scale DL by RNRM*PROD to obtain the residual RL. @@ -6405,17 +6456,20 @@ CALL SSCAL(NEQ, TEM, DL, 1) 200 CONTINUE LL = LGMR LLP1 = LL + 1 - DO 210 K = 1,LLP1 - 210 R(K) = 0.0E0 + DO K = 1,LLP1 + R(K) = 0.0E0 + END DO R(1) = RNRM CALL SHELS (HES, MAXLP1, LL, Q, R) - DO 220 K = 1,NEQ - 220 Z(K) = 0.0E0 - DO 230 I = 1,LL + DO K = 1,NEQ + Z(K) = 0.0E0 + END DO + DO I = 1,LL CALL SAXPY (NEQ, R(I), V(1,I), 1, Z, 1) - 230 CONTINUE - DO 240 I = 1,NEQ - 240 Z(I) = Z(I)/WGHT(I) + END DO + DO I = 1,NEQ + Z(I) = Z(I)/WGHT(I) + END DO C Load RHO into RHOK. RHOK = RHO RETURN @@ -6525,16 +6579,18 @@ DIMENSION Y(*), YPRIME(*), SAVR(*), V(*), WGHT(*), YPTEM(*), C----------------------------------------------------------------------- C Set VTEM = D * V. C----------------------------------------------------------------------- - DO 10 I = 1,NEQ - 10 VTEM(I) = V(I)/WGHT(I) + DO I = 1,NEQ + VTEM(I) = V(I)/WGHT(I) + END DO IER = 0 C----------------------------------------------------------------------- C Store Y in Z and increment Z by VTEM. C Store YPRIME in YPTEM and increment YPTEM by VTEM*CJ. C----------------------------------------------------------------------- - DO 20 I = 1,NEQ + DO I = 1,NEQ YPTEM(I) = YPRIME(I) + VTEM(I)*CJ*scaleps - 20 Z(I) = Y(I) + VTEM(I)*scaleps + Z(I) = Y(I) + VTEM(I)*scaleps + END DO C----------------------------------------------------------------------- C Call RES with incremented Y, YPRIME arguments C stored in Z, YPTEM. VTEM is overwritten with new residual. @@ -6547,8 +6603,9 @@ CALL RES(TN,Z,YPTEM,CJ,VTEM,IRES,RPAR,IPAR) C Set Z = (dF/dY) * VBAR using difference quotient. C (VBAR is old value of VTEM before calling RES) C----------------------------------------------------------------------- - DO 70 I = 1,NEQ - 70 Z(I) = (VTEM(I) - SAVR(I))/scaleps + DO I = 1,NEQ + Z(I) = (VTEM(I) - SAVR(I))/scaleps + END DO C----------------------------------------------------------------------- C Apply inverse of left preconditioner to Z. C----------------------------------------------------------------------- @@ -6559,8 +6616,9 @@ CALL PSOL (NEQ, TN, Y, YPRIME, SAVR, YPTEM, CJ, WGHT, WP, IWP, C----------------------------------------------------------------------- C Apply D-inverse to Z and return. C----------------------------------------------------------------------- - DO 90 I = 1,NEQ - 90 Z(I) = Z(I)*WGHT(I) + DO I = 1,NEQ + Z(I) = Z(I)*WGHT(I) + END DO RETURN C C------END OF SUBROUTINE SATV------------------------------------------- From bc8bf0da4f57cd908fad578a91c11daa4392ffb4 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 15:59:20 -0700 Subject: [PATCH 10/13] Removes missed DO termination statement warning from daspk.m --- svr/daspk.m | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/svr/daspk.m b/svr/daspk.m index 5967cd3a..00ddb266 100755 --- a/svr/daspk.m +++ b/svr/daspk.m @@ -1897,9 +1897,10 @@ C Load H and RWORK(LH) with H0. C Load Y and H*YPRIME into PHI(*,1) and PHI(*,2). C ITEMP = LPHI + NEQ - DO 380 I = 1,NEQ + DO I = 1,NEQ RWORK(LPHI + I - 1) = Y(I) -380 RWORK(ITEMP + I - 1) = H*YPRIME(I) + RWORK(ITEMP + I - 1) = H*YPRIME(I) + END DO C GO TO 500 C @@ -2984,8 +2985,9 @@ ccc call xerrab("") IF(KNEW.EQ.KM1)GO TO 540 IF(K.EQ.IWM(LMXORD)) GO TO 550 IF(KP1.GE.NS.OR.KDIFF.EQ.1)GO TO 550 - DO 510 I=1,NEQ -510 DELTA(I)=E(I)-PHI(I,KP2) + DO I=1,NEQ + DELTA(I)=E(I)-PHI(I,KP2) + END DO ERKP1 = (1.0E0/(K+2))*SDWNRM(NEQ,DELTA,VT,RPAR,IPAR) TERKP1 = (K+2)*ERKP1 IF(K.GT.1)GO TO 520 From 5db73ee850cfc83d487afc4a5b8bbc05a5b96863 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 16:20:19 -0700 Subject: [PATCH 11/13] Fixes Shared DO termination label warnings in flxcomp.m --- flx/flxcomp.m | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/flx/flxcomp.m b/flx/flxcomp.m index d2536cd9..1a529432 100755 --- a/flx/flxcomp.m +++ b/flx/flxcomp.m @@ -20,13 +20,14 @@ call refine c *** clear the working arrays - do 20 kk=1,jdim + do kk=1,jdim ilast(kk)=0 npoint(kk)=0 - do 20 ll=1,npts - xcurve(ll,kk)=0.0 - ycurve(ll,kk)=0.0 - 20 continue + do ll=1,npts + xcurve(ll,kk)=0.0 + ycurve(ll,kk)=0.0 + end do + end do c initialize the contour indexing variables ncmin=0 From 30adb4aa177188e69899da6f194ac8d1a2c54778 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 17:06:10 -0700 Subject: [PATCH 12/13] Fixes DO-loop errors for comutil.m Still need to address: - Arithmetic IF-statements - Replace ASSIGN statements - Assigned GOTO-statement --- com/comutil.m | 203 +++++++++++++++++++++++++++++--------------------- 1 file changed, 117 insertions(+), 86 deletions(-) diff --git a/com/comutil.m b/com/comutil.m index 1aed0001..adcd39a1 100755 --- a/com/comutil.m +++ b/com/comutil.m @@ -1501,12 +1501,13 @@ CALL INTRhV(TY,NY+KY,YVAL,ILOY,LEFTY,MFLAG) C ... GET COEFFICIENTS OF NONZERO BASIS FUNCTIONS IN (3) C IPOS = IXY - 1 - DO 50 K=KMIN,KMAX - DO 50 J=JMIN,JMAX + DO K=KMIN,KMAX + DO J=JMIN,JMAX IPOS = IPOS + 1 WORK(IPOS) = B1VAhL(XVAL,IDX,TX,NX,KX,FCN(1,J,K),INX, + WORK(IW),IERR) - 50 CONTINUE + END DO + END DO IHAVEX = 1 C ENDIF @@ -5646,40 +5647,43 @@ CALL XERMShG ('SLATEC', 'BNDAChC', 'MDG.LT.IR, PROBABLE ERROR.', C ALG. STEPS 6-7 IF (JT.LE.IR) GO TO 30 C ALG. STEPS 8-9 - DO 10 I=1,MT + DO I=1,MT IG1=JT+MT-I IG2=IR+MT-I - DO 10 J=1,NBP1 - G(IG1,J)=G(IG2,J) - 10 CONTINUE + DO J=1,NBP1 + G(IG1,J)=G(IG2,J) + END DO + END DO C ALG. STEP 10 IE=JT-IR - DO 20 I=1,IE + DO I=1,IE IG=IR+I-1 - DO 20 J=1,NBP1 - G(IG,J)=ZERO - 20 CONTINUE + DO J=1,NBP1 + G(IG,J)=ZERO + END DO + END DO C ALG. STEP 11 IR=JT C ALG. STEP 12 30 MU=MIN(NB-1,IR-IP-1) IF (MU.EQ.0) GO TO 60 C ALG. STEP 13 - DO 50 L=1,MU + DO L=1,MU C ALG. STEP 14 K=MIN(L,JT-IP) C ALG. STEP 15 LP1=L+1 IG=IP+L - DO 40 I=LP1,NB + DO I=LP1,NB JG=I-K G(IG,JG)=G(IG,I) - 40 CONTINUE + END DO C ALG. STEP 16 - DO 50 I=1,K - JG=NBP1-I - G(IG,JG)=ZERO - 50 CONTINUE + DO I=1,K + JG=NBP1-I + G(IG,JG)=ZERO + END DO + END DO C ALG. STEP 17 60 IP=JT C ALG. STEPS 18-19 @@ -5908,21 +5912,21 @@ real G(MDG,*),X(*),zero,rnorm,rsq,s GO TO (10,90,50), MODE C ********************* MODE = 1 C ALG. STEP 26 - 10 DO 20 J=1,N + 10 DO J=1,N X(J)=G(J,NB+1) - 20 CONTINUE + END DO RSQ=ZERO NP1=N+1 IRM1=IR-1 IF (NP1.GT.IRM1) GO TO 40 - DO 30 J=NP1,IRM1 + DO J=NP1,IRM1 RSQ=RSQ+G(J,NB+1)**2 - 30 CONTINUE + END DO RNORM=SQRT(RSQ) 40 CONTINUE C ********************* MODE = 3 C ALG. STEP 27 - 50 DO 80 II=1,N + 50 DO II=1,N I=N+1-II C ALG. STEP 28 S=ZERO @@ -5931,29 +5935,31 @@ GO TO (10,90,50), MODE IF (I.EQ.N) GO TO 70 C ALG. STEP 30 IE=MIN(N+1-I,NB) - DO 60 J=2,IE + DO J=2,IE JG=J+L IX=I-1+J S=S+G(I,JG)*X(IX) - 60 CONTINUE + END DO C ALG. STEP 31 70 IF (G(I,L+1)) 80,130,80 80 X(I)=(X(I)-S)/G(I,L+1) + END DO C ALG. STEP 32 RETURN C ********************* MODE = 2 - 90 DO 120 J=1,N + 90 DO J=1,N S=ZERO IF (J.EQ.1) GO TO 110 I1=MAX(1,J-NB+1) I2=J-1 - DO 100 I=I1,I2 + DO I=I1,I2 L=J-I+1+MAX(0,I-IP) S=S+X(I)*G(I,L) - 100 CONTINUE + END DO 110 L=MAX(0,J-IP) IF (G(J,L+1)) 120,130,120 120 X(J)=(X(J)-S)/G(J,L+1) + END DO RETURN C 130 CONTINUE @@ -6000,20 +6006,23 @@ real T(*),VNIKX(K,*),A(20,20),x,fkmd,diff,v CALL BSPLVhN(T,K+1-NDERIV,1,X,ILEFT,VNIKX(NDERIV,NDERIV)) IF (NDERIV .LE. 1) GO TO 99 IDERIV = NDERIV - DO 15 I=2,NDERIV + DO I=2,NDERIV IDERVM = IDERIV-1 - DO 11 J=IDERIV,K - 11 VNIKX(J-1,IDERVM) = VNIKX(J,IDERIV) + DO J=IDERIV,K + VNIKX(J-1,IDERVM) = VNIKX(J,IDERIV) + END DO IDERIV = IDERVM CALL BSPLVhN(T,0,2,X,ILEFT,VNIKX(IDERIV,IDERIV)) - 15 CONTINUE -C - DO 20 I=1,K - DO 19 J=1,K - 19 A(I,J) = 0. - 20 A(I,I) = 1. + END DO +C + DO I=1,K + DO J=1,K + A(I,J) = 0. + END DO + A(I,I) = 1. + END DO KMD = K - DO 40 M=2,NDERIV + DO M=2,NDERIV KMD = KMD-1 FKMD = KMD I = ILEFT @@ -6023,20 +6032,24 @@ CALL BSPLVhN(T,0,2,X,ILEFT,VNIKX(IDERIV,IDERIV)) DIFF = T(IPKMD) - T(I) IF (JM1 .EQ. 0) GO TO 26 IF (DIFF .EQ. 0.) GO TO 25 - DO 24 L=1,J - 24 A(L,J) = (A(L,J) - A(L,J-1))/DIFF*FKMD + DO L=1,J + A(L,J) = (A(L,J) - A(L,J-1))/DIFF*FKMD + END DO 25 J = JM1 I = I - 1 GO TO 21 26 IF (DIFF .EQ. 0.) GO TO 30 A(1,1) = A(1,1)/DIFF*FKMD C - 30 DO 40 I=1,K + 30 DO I=1,K V = 0. JLOW = MAX(I,M) - DO 35 J=JLOW,K - 35 V = A(I,J)*VNIKX(J,M) + V - 40 VNIKX(I,M) = V + DO J=JLOW,K + V = A(I,J)*VNIKX(J,M) + V + END DO + VNIKX(I,M) = V + END DO + END DO 99 RETURN END c----------------------------------------------------------------------- @@ -6079,11 +6092,12 @@ GO TO (10,20),INDEX DELTAM(J) = X - T(IMJP1) VMPREV = 0. JP1 = J+1 - DO 26 L=1,J + DO L=1,J JP1ML = JP1-L VM = VNIKX(L)/(DELTAP(L) + DELTAM(JP1ML)) VNIKX(L) = VM*DELTAP(L) + VMPREV - 26 VMPREV = VM*DELTAM(JP1ML) + VMPREV = VM*DELTAM(JP1ML) + END DO VNIKX(JP1) = VMPREV J = JP1 IF (J .LT. JHIGH) GO TO 20 @@ -6150,13 +6164,15 @@ real U(IUE,*), C(*),cl,sm,one,clinv,up,b,ul1m1,sdot CL=ABS(U(1,LPIVOT)) IF (MODE.EQ.2) GO TO 60 C ****** CONSTRUCT THE TRANSFORMATION. ****** - DO 10 J=L1,M - 10 CL=MAX(ABS(U(1,J)),CL) + DO J=L1,M + CL=MAX(ABS(U(1,J)),CL) + END DO IF (CL) 130,130,20 20 CLINV=ONE/CL SM=(U(1,LPIVOT)*CLINV)**2 - DO 30 J=L1,M - 30 SM=SM+(U(1,J)*CLINV)**2 + DO J=L1,M + SM=SM+(U(1,J)*CLINV)**2 + END DO CL=CL*SQRT(SM) IF (U(1,LPIVOT)) 50,50,40 40 CL=-CL @@ -6181,15 +6197,17 @@ real U(IUE,*), C(*),cl,sm,one,clinv,up,b,ul1m1,sdot I3=I2+INCR I4=I3 SM=C(I2)*up - DO 90 I=L1,M + DO I=L1,M SM=SM+C(I3)*U(1,I) - 90 I3=I3+ICE + I3=I3+ICE + END DO IF (SM) 100,120,100 100 SM=SM*B C(I2)=C(I2)+SM*up - DO 110 I=L1,M + DO I=L1,M C(I4)=C(I4)+SM*U(1,I) - 110 I4=I4+ICE + I4=I4+ICE + END DO 120 CONTINUE 130 RETURN 140 CONTINUE @@ -6390,27 +6408,28 @@ CALL XERMShG ('SLATEC', 'HFTIh', RETURN 6 CONTINUE C - DO 80 J=1,LDIAG + DO J=1,LDIAG IF (J.EQ.1) GO TO 20 C C UPDATE SQUARED COLUMN LENGTHS AND FIND LMAX C .. LMAX=J - DO 10 L=J,N + DO L=J,N H(L)=H(L)-A(J-1,L)**2 IF (H(L).GT.H(LMAX)) LMAX=L - 10 CONTINUE + END DO IF (FACTOR*H(LMAX) .GT. HMAX*RELEPS) GO TO 50 C C COMPUTE SQUARED COLUMN LENGTHS AND FIND LMAX C .. 20 LMAX=J - DO 40 L=J,N + DO L=J,N H(L)=0. - DO 30 I=J,M - 30 H(L)=H(L)+A(I,L)**2 + DO I=J,M + H(L)=H(L)+A(I,L)**2 + END DO IF (H(L).GT.H(LMAX)) LMAX=L - 40 CONTINUE + END DO HMAX=H(LMAX) C .. C LMAX HAS BEEN DETERMINED @@ -6420,22 +6439,24 @@ CALL XERMShG ('SLATEC', 'HFTIh', 50 CONTINUE IP(J)=LMAX IF (IP(J).EQ.J) GO TO 70 - DO 60 I=1,M + DO I=1,M TMP=A(I,J) A(I,J)=A(I,LMAX) - 60 A(I,LMAX)=TMP + A(I,LMAX)=TMP + END DO H(LMAX)=H(J) C C COMPUTE THE J-TH TRANSFORMATION AND APPLY IT TO A AND B. C .. 70 CALL H12h (1,J,J+1,M,A(1,J),1,H(J),A(1,J+1),1,MDA,N-J) - 80 CALL H12h (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB) + CALL H12h (2,J,J+1,M,A(1,J),1,H(J),B,1,MDB,NB) + END DO C C DETERMINE THE PSEUDORANK, K, USING THE TOLERANCE, TAU. C .. - DO 90 J=1,LDIAG - IF (ABS(A(J,J)).LE.TAU) GO TO 100 - 90 CONTINUE + DO J=1,LDIAG + IF (ABS(A(J,J)).LE.TAU) GO TO 100 + END DO K=LDIAG GO TO 110 100 K=J-1 @@ -6444,58 +6465,67 @@ CALL XERMShG ('SLATEC', 'HFTIh', C COMPUTE THE NORMS OF THE RESIDUAL VECTORS. C IF (NB.LE.0) GO TO 140 - DO 130 JB=1,NB + DO JB=1,NB TMP=SZERO IF (KP1.GT.M) GO TO 130 - DO 120 I=KP1,M - 120 TMP=TMP+B(I,JB)**2 + DO I=KP1,M + TMP=TMP+B(I,JB)**2 + END DO 130 RNORM(JB)=SQRT(TMP) + END DO 140 CONTINUE C SPECIAL FOR PSEUDORANK = 0 IF (K.GT.0) GO TO 160 IF (NB.LE.0) GO TO 270 - DO 150 JB=1,NB - DO 150 I=1,N - 150 B(I,JB)=SZERO + DO JB=1,NB + DO I=1,N + B(I,JB)=SZERO + END DO + END DO GO TO 270 C C IF THE PSEUDORANK IS LESS THAN N COMPUTE HOUSEHOLDER C DECOMPOSITION OF FIRST K ROWS. C .. 160 IF (K.EQ.N) GO TO 180 - DO 170 II=1,K + DO II=1,K I=KP1-II - 170 CALL H12h (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1) + CALL H12h (1,I,KP1,N,A(I,1),MDA,G(I),A,MDA,1,I-1) + END DO 180 CONTINUE C C IF (NB.LE.0) GO TO 270 - DO 260 JB=1,NB + DO JB=1,NB C C SOLVE THE K BY K TRIANGULAR SYSTEM. C .. - DO 210 L=1,K + DO L=1,K SM=DZERO I=KP1-L IF (I.EQ.K) GO TO 200 IP1=I+1 - DO 190 J=IP1,K - 190 SM=SM+A(I,J)*DBLE(B(J,JB)) + DO J=IP1,K + SM=SM+A(I,J)*DBLE(B(J,JB)) + END DO 200 SM1=SM - 210 B(I,JB)=(B(I,JB)-SM1)/A(I,I) + B(I,JB)=(B(I,JB)-SM1)/A(I,I) + END DO C C COMPLETE COMPUTATION OF SOLUTION VECTOR. C .. IF (K.EQ.N) GO TO 240 - DO 220 J=KP1,N - 220 B(J,JB)=SZERO - DO 230 I=1,K - 230 CALL H12h (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1) + DO J=KP1,N + B(J,JB)=SZERO + END DO + DO I=1,K + CALL H12h (2,I,KP1,N,A(I,1),MDA,G(I),B(1,JB),1,MDB,1) + END DO C C RE-ORDER THE SOLUTION VECTOR TO COMPENSATE FOR THE C COLUMN INTERCHANGES. C .. - 240 DO 250 JJ=1,LDIAG + 240 DO JJ=1,LDIAG J=LDIAG+1-JJ IF (IP(J).EQ.J) GO TO 250 L=IP(J) @@ -6503,7 +6533,8 @@ CALL XERMShG ('SLATEC', 'HFTIh', B(L,JB)=B(J,JB) B(J,JB)=TMP 250 CONTINUE - 260 CONTINUE + END DO + END DO C .. C THE SOLUTION VECTORS, X, ARE NOW C IN THE FIRST N ROWS OF THE ARRAY B(,). From 8b9e3ddf1e41535d1eacf006f356d9a4589bdd64 Mon Sep 17 00:00:00 2001 From: Andreas Holm <60451789+holm10@users.noreply.github.com> Date: Wed, 18 Oct 2023 17:33:59 -0700 Subject: [PATCH 13/13] Remove obsolete definition of condition in exmain.c --- bbb/exmain.c | 2 -- 1 file changed, 2 deletions(-) diff --git a/bbb/exmain.c b/bbb/exmain.c index 6828f483..751bf78e 100755 --- a/bbb/exmain.c +++ b/bbb/exmain.c @@ -37,8 +37,6 @@ void int_handler() { printf("or a single line to be evaluated by Python.\n"); #pragma omp master { -int condition; -condition=1; while(1){ #ifdef HAS_READLINE ret = readline("Debug>>> ");