@@ -528,7 +528,7 @@ SUBROUTINE ros_ErrorMsg(Code,T,H,IERR)
528528 CASE (- 6 )
529529 PRINT * , ' --> No of steps exceeds maximum bound'
530530 CASE (- 7 )
531- PRINT * , ' --> Step size too small: T + 10 *H = T' , &
531+ PRINT * , ' --> Step size too small: T + 0.1 *H = T' , &
532532 ' or H < Roundoff'
533533 CASE (- 8 )
534534 PRINT * , ' --> Matrix is repeatedly singular'
@@ -653,39 +653,42 @@ SUBROUTINE ros_Integrator (Y, Tstart, Tend, T, &
653653 RETURN
654654 END IF
655655
656- ! ~~~> Compute the stages
657- Stage: DO istage = 1 , ros_S
658-
656+ ! ~~~> Compute the stages
657+ Stage: DO istage = 1 , ros_S
658+
659659 ! Current istage offset. Current istage vector is K(ioffset+1:ioffset+N)
660- ioffset = N* (istage-1 )
660+ ioffset = N* (istage-1 )
661661
662662 ! For the 1st istage the function has been computed previously
663- IF ( istage == 1 ) THEN
663+ IF ( istage == 1 ) THEN
664664 Fcn(1 :N) = Fcn0(1 :N)
665- ! istage>1 and a new function evaluation is needed at the current istage
666- ELSEIF ( ros_NewF(istage) ) THEN
665+ ! istage>1 and a new function evaluation is needed at the current istage
666+ ELSEIF ( ros_NewF(istage) ) THEN
667667 Ynew(1 :N) = Y(1 :N)
668668 DO j = 1 , istage-1
669- Ynew(1 :N) = Ynew(1 :N) &
670- + ros_A((istage-1 )* (istage-2 )/ 2 + j) * K(N* (j-1 )+ 1 :N* j)
669+ Ynew(1 :N) = Ynew(1 :N) + &
670+ ros_A((istage-1 )* (istage-2 )/ 2 + j) * K(N* (j-1 )+ 1 :N* j)
671671 END DO
672672 Tau = T + ros_Alpha(istage)* Direction* H
673673 CALL FunTemplate( Tau, Ynew, Fcn )
674674 ISTATUS(Nfun) = ISTATUS(Nfun) + 1
675- END IF ! if istage == 1 elseif ros_NewF(istage)
676- K(ioffset+1 :ioffset+ N) = Fcn(1 :N)
677- DO j = 1 , istage-1
675+ END IF ! if istage == 1 elseif ros_NewF(istage)
676+ K(ioffset+1 :ioffset+ N) = Fcn(1 :N)
677+ DO j = 1 , istage-1
678678 HC = ros_C((istage-1 )* (istage-2 )/ 2 + j)/ (Direction* H)
679- K(ioffset+1 :ioffset+ N) = K(ioffset+1 :ioffset+ N) + HC * K(N* (j-1 )+ 1 :N* j)
680- END DO
681- IF ((.NOT. Autonomous).AND. (ros_Gamma(istage).NE. ZERO)) THEN
679+ K(ioffset+1 :ioffset+ N) = K(ioffset+1 :ioffset+ N) &
680+ + HC * K(N* (j-1 )+ 1 :N* j)
681+ END DO
682+ IF ((.NOT. Autonomous).AND. (ros_Gamma(istage).NE. ZERO)) THEN
682683 HG = Direction* H* ros_Gamma(istage)
683- K(ioffset+1 :ioffset+ N) = K(ioffset+1 :ioffset+ N) + HG * dFdT(1 :N)
684- END IF
685- CALL ros_Solve(Ghimj, Pivot, K(ioffset+1 ))
684+ K(ioffset+1 :ioffset+ N) = K(ioffset+1 :ioffset+ N) &
685+ + HG * dFdT(1 :N)
686+ END IF
687+ CALL ros_Solve(Ghimj, Pivot, K(ioffset+1 ))
686688
687689 END DO Stage
688690
691+
689692! ~~~> Compute the new solution
690693 Ynew(1 :N) = Y(1 :N)
691694 DO j= 1 ,ros_S
@@ -1015,7 +1018,7 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, &
10151018 ! istage>1 and a new function evaluation is needed at the current istage
10161019 ! K = 0.0_dp ! is this fix needed? hplin 14:04 -- not. 3 hours wiser later
10171020 ELSEIF ( ros_NewF(istage) ) THEN
1018- Fcn (1 :N) = Fcn0 (1 :N)
1021+ Ynew (1 :N) = Y (1 :N)
10191022 DO j = 1 , istage-1
10201023 ! In full vector space. Just use WAXPY as normal
10211024 ! other entries in K are set to 1 previously.
@@ -1081,7 +1084,8 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, &
10811084
10821085 ! faster version:
10831086 DO i = 1 ,rNVAR
1084- K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) + HC * K(N* (j-1 )+ SPC_MAP(i))
1087+ K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) &
1088+ + HC * K(N* (j-1 )+ SPC_MAP(i))
10851089 ENDDO
10861090 ! CALL zWAXPY(N,HC,K(N*(j-1)+1),K(ioffset+1),SPC_MAP)
10871091 ! loop unrolling is consistently slower here. 18:58
@@ -1093,7 +1097,8 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, &
10931097 ! full version: CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1)
10941098 ! faster version:
10951099 DO i = 1 ,rNVAR
1096- K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) + HG * dFdT(SPC_MAP(i))
1100+ K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) &
1101+ + HG * dFdT(SPC_MAP(i))
10971102 ENDDO
10981103 ENDIF
10991104
@@ -1144,7 +1149,7 @@ SUBROUTINE ros_yIntegrator (Y, Tstart, Tend, T, &
11441149 ISTATUS(Nstp) = ISTATUS(Nstp) + 1
11451150 IF ( (Err <= ONE).OR. (H <= Hmin) ) THEN ! ~~~> Accept step
11461151 ISTATUS(Nacc) = ISTATUS(Nacc) + 1
1147- Ynew (1 :N) = Y (1 :N)
1152+ Y (1 :N) = Ynew (1 :N)
11481153 T = T + Direction* H
11491154 Hnew = MAX (Hmin,MIN (Hnew,Hmax))
11501155 IF (RejectLastH) THEN ! No step size increase after a rejected step
@@ -1542,7 +1547,8 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, &
15421547
15431548 ! faster version:
15441549 DO i = 1 ,rNVAR
1545- K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) + HC * K(N* (j-1 )+ SPC_MAP(i))
1550+ K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) &
1551+ + HC * K(N* (j-1 )+ SPC_MAP(i))
15461552 ENDDO
15471553 ! CALL zWAXPY(N,HC,K(N*(j-1)+1),K(ioffset+1),SPC_MAP)
15481554 ! loop unrolling is consistently slower here. 18:58
@@ -1554,7 +1560,8 @@ SUBROUTINE ros_yIntegratorA (Y, Tstart, Tend, T, &
15541560 ! full version: CALL WAXPY(N,HG,dFdT,1,K(ioffset+1),1)
15551561 ! faster version:
15561562 DO i = 1 ,rNVAR
1557- K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) + HG * dFdT(SPC_MAP(i))
1563+ K(ioffset+ SPC_MAP(i)) = K(ioffset+ SPC_MAP(i)) &
1564+ + HG * dFdT(SPC_MAP(i))
15581565 ENDDO
15591566 ENDIF
15601567
0 commit comments