Demo entry 6036664

n1int

   

Submitted by anonymous on Sep 20, 2016 at 16:22
Language: Fortran. Code size: 25.9 kB.

2       CONTINUE
C         --------------------------------------------------------
C         2 Startup of iteration step
          IF (.NOT.QJCRFR) THEN
C           ------------------------------------------------------
C           2.1 Scaling of variables X(N)
            CALL N1SCAL(N,X,XA,XSCAL,XW,ISCAL,QINISC,IOPT,LRWK,RWK)
            QINISC = .FALSE.
            IF(NITER.NE.0)THEN
C             ----------------------------------------------------
C             2.2 Aposteriori estimate of damping factor
              DO 2200 L1=1,N
                DXQA(L1)=DXQ(L1)
2200          CONTINUE
              IF (.NOT.QORDI) THEN
                FCNUMP = ZERO
                DO 2201 L1=1,N
                  FCNUMP=FCNUMP+(DX(L1)/XW(L1))**2
2201            CONTINUE
                TH = FC-ONE
                FCDNM = ZERO
                DO 2202 L1=1,N
                  FCDNM=FCDNM+((DXQA(L1)+TH*DX(L1))/XW(L1))**2
2202            CONTINUE
C               --------------------------------------------------
C               2.2.2 Decision criterion for Jacobian updating
C                     technique:
C                     QGENJ.EQ..TRUE. numerical differentation,
C                     QGENJ.EQ..FALSE. rank1 updating
                QGENJ = .TRUE.
                IF (FC.EQ.FCPRI) THEN
                  QGENJ = FC.LT.ONE.OR.FCA.LT.ONE.OR.DMYCOR.LE.FC*SIGMA
     $                    .OR. .NOT.QRANK1 .OR. NEW+2.GT.NBROY 
                  FCA = FC
                ELSE
                  DMYCOR = FCA*FCA*HALF*DSQRT(FCNUMP/FCDNM)
                  IF (NONLIN.LE.3) THEN
                    FCCOR = DMIN1(ONE,DMYCOR)
                  ELSE
                    FCCOR = DMIN1(ONE,HALF*DMYCOR)
                  ENDIF
                  FCA = DMAX1(DMIN1(FC,FCCOR),FCMIN)
C$Test-begin
                  IF (MPRMON.GE.5) THEN
                    WRITE(LUMON,22201) FCCOR, FC, DMYCOR, FCNUMP,
     $                                 FCDNM
22201               FORMAT (/, ' +++ aposteriori estimate +++', /,
     $                    ' FCCOR  = ', D18.10, '  FC     = ', D18.10, /,
     $                    ' DMYCOR = ', D18.10, '  FCNUMP = ', D18.10, /,
     $                    ' FCDNM  = ', D18.10, /,
     $                       ' ++++++++++++++++++++++++++++', /)
                  ENDIF
C$Test-end 
                ENDIF
                FCK2 = FCA
C               ------------------------------------------------------
C               2.2.1 Computation of the numerator of damping
C                     factor predictor
                FCNMP2 = ZERO
                DO 221 L1=1,N
                  FCNMP2=FCNMP2+(DXQA(L1)/XW(L1))**2
221             CONTINUE
                FCNUMP = FCNUMP*FCNMP2
              ENDIF
            ENDIF
          ENDIF
          QJCRFR =.FALSE.
C         --------------------------------------------------------
C         2.3 Jacobian matrix (stored to array A(M1,N))
C         --------------------------------------------------------
C         2.3.1 Jacobian generation by routine JAC or
C               difference approximation (If QGENJ.EQ..TRUE.)
C               - or -
C               Rank-1 update of Jacobian (If QGENJ.EQ..FALSE.)
          IF (QGENJ .AND. (.NOT.QSIMPL .OR. NITER.EQ.0)) THEN
            NEW = 0
            IF (JACGEN.EQ.1) THEN
               IF (MPRTIM.NE.0) CALL MONON(2)
               CALL JAC(N,M1,X,A,IFAIL)
               IF (MPRTIM.NE.0) CALL MONOFF(2)
            ELSE
              IF (MSTOR.EQ.0) THEN
                IF (MPRTIM.NE.0) CALL MONON(2)
                IF (JACGEN.EQ.3) 
     $            CALL N1JCF(FCN,N,N,X,F,A,XW,ETA,ETAMIN,ETAMAX,
     $                       ETADIF,CONV,NFCNJ,T1,IFAIL)
                IF (JACGEN.EQ.2) 
     $            CALL N1JAC(FCN, N, N, X, F, A, XW, AJDEL, AJMIN,
     $                       NFCNJ, T1, IFAIL)
                IF (MPRTIM.NE.0) CALL MONOFF(2)
              ELSE IF (MSTOR.EQ.1) THEN
                IF (MPRTIM.NE.0) CALL MONON(2)
                IF (JACGEN.EQ.3) 
     $            CALL N1JCFB(FCN,N,M1,ML,X,F,A,XW,ETA,ETAMIN,
     $                        ETAMAX,ETADIF,CONV,NFCNJ,T1,T2,T3,IFAIL)
                IF (JACGEN.EQ.2) 
     $            CALL N1JACB (FCN, N, M1, ML, X, F, A, XW,  AJDEL,
     $                         AJMIN, NFCNJ, T1, T2, T3, IFAIL)
                IF (MPRTIM.NE.0) CALL MONOFF(2)
              ENDIF
            ENDIF
            NJAC = NJAC + 1
C     Exit, If ...
            IF (JACGEN.EQ.1 .AND. IFAIL.LT.0) THEN
              IERR = 83
              GOTO 4299
            ENDIF
            IF (JACGEN.NE.1 .AND. IFAIL.NE.0) THEN
              IERR = 82
              GOTO 4299
            ENDIF
          ELSE IF (.NOT.QSIMPL) THEN
            NEW = NEW+1
          ENDIF
          IF ( NEW.EQ.0 .AND. (QLU.OR.NITER.EQ.0) ) THEN
C             ------------------------------------------------------
C             2.3.2.1 Save scaling values
              DO 2321 L1=1,N
                XWA(L1) = XW(L1)
2321          CONTINUE
C             ------------------------------------------------------
C             2.3.2.2 Prepare Jac. for use by band-solver DGBFA/DGBSL
              IF (MSTOR.EQ.1) THEN
                DO 2322 L1=1,N
                  DO 2323 L2=M2,1,-1
                    A(L2+ML,L1)=A(L2,L1)
2323              CONTINUE
2322            CONTINUE
              ENDIF
C             ------------------------------------------------------
C             2.4 Prepare solution of the linear system
C             ------------------------------------------------------
C             2.4.1 internal column scaling of matrix A
              IF (MSTOR.EQ.0) THEN
                DO 2410 K=1,N
                  S1 =-XW(K)
                  DO 2412 L1=1,N
                    A(L1,K)=A(L1,K)*S1
2412              CONTINUE
2410            CONTINUE
              ELSE IF (MSTOR.EQ.1) THEN
                DO 2413 K=1,N
                  L2=MAX0(1+M2-K,ML+1)
                  L3=MIN0(N+M2-K,M1)
                  S1 =-XW(K)
                  DO 2415 L1=L2,L3
                    A(L1,K)=A(L1,K)*S1
2415              CONTINUE
2413            CONTINUE
              ENDIF
C             ----------------------------------------------------
C             2.4.2 Row scaling of matrix A
              IF (QSCALE) THEN
                IF (MSTOR.EQ.0) THEN
                  CALL N1SCRF(N,N,A,FW)
                ELSE IF (MSTOR.EQ.1) THEN
                  CALL N1SCRB(N,M1,ML,MU,A,FW)
                ENDIF
              ELSE
                DO 242 K=1,N
                  FW(K)=ONE
242             CONTINUE
              ENDIF
          ENDIF
C         --------------------------------------------------------
C         2.4.3 Save and scale values of F(N)
          DO 243 L1=1,N
            FA(L1)=F(L1)
            T1(L1)=F(L1)*FW(L1)
243       CONTINUE
C         --------------------------------------------------------
C         3 Central part of iteration step
C         --------------------------------------------------------
C         3.1 Solution of the linear system
C         --------------------------------------------------------
C         3.1.1 Decomposition of (N,N)-matrix A
          IF (.NOT.QLINIT) THEN
            NIWLA = IWK(18)+1
            NRWLA = IWK(19)+1
            LIWL = LIWK-NIWLA+1
            LRWL = LRWK-NRWLA+1
          ENDIF
          IF ( NEW.EQ.0 .AND. (QLU.OR.NITER.EQ.0)) THEN
            IF (MPRTIM.NE.0) CALL MONON(3)
            CALL N1FACT(N,M1,ML,MU,A,IOPT,IFAIL,LIWL,IWK(NIWLA),
     $                  NILUSE,LRWL,RWK(NRWLA),NRLUSE)
            IF (MPRTIM.NE.0) CALL MONOFF(3)
            IF (.NOT.QLINIT) THEN
              NIWKFR = NIWKFR+NILUSE
              NRWKFR = NRWKFR+NRLUSE
C             Store lengths of currently required workspaces
              IWK(18) = NIWKFR-1
              IWK(19) = NRWKFR-1
            ENDIF
C       Exit Repeat If ...
            IF(IFAIL.NE.0) THEN
              IF (IFAIL.EQ.1) THEN
                IERR = 1
              ELSE
                IERR = 80
              ENDIF
              GOTO 4299
            ENDIF
          ENDIF
          QLINIT = .TRUE.
C         --------------------------------------------------------
C         3.1.2 Solution of linear (N,N)-system
          IF(NEW.EQ.0) THEN 
            IF (MPRTIM.NE.0) CALL MONON(4)
            CALL N1SOLV(N,M1,ML,MU,A,T1,IOPT,IFAIL,LIWL,IWK(NIWLA),
     $                 IDUMMY,LRWL,RWK(NRWLA),IDUMMY)
            IF (MPRTIM.NE.0) CALL MONOFF(4)
C     Exit Repeat If ...
            IF(IFAIL.NE.0)  THEN
              IERR = 81
              GOTO 4299
            ENDIF
          ELSE  
            ALFA1=ZERO
            ALFA2=ZERO
            DO 3121 I=1,N
              ALFA1=ALFA1+DX(I)*DXQ(I)/XW(I)**2
              ALFA2=ALFA2+DX(I)**2/XW(I)**2
3121        CONTINUE
            ALFA=ALFA1/ALFA2
            BETA=ONE-ALFA
            DO 3122 I=1,N
              T1(I)=(DXQ(I)+(FCA-ONE)*ALFA*DX(I))/BETA
3122        CONTINUE
            IF(NEW.EQ.1) THEN
              DO 3123 I=1,N
                DXSAVE(I,1)=DX(I)
3123          CONTINUE
            ENDIF
            DO 3124 I=1,N
              DXSAVE(I,NEW+1)=T1(I)
              DX(I)=T1(I)
              T1(I)=T1(I)/XW(I)
3124        CONTINUE
          ENDIF
C         --------------------------------------------------------
C         3.2 Evaluation of scaled natural level function SUMX
C             scaled maximum error norm CONV
C             evaluation of (scaled) standard level function
C             DLEVF ( DLEVF only, if MPRMON.GE.2 )
C             and computation of ordinary Newton corrections 
C             DX(N)
          IF (.NOT. QSIMPL) THEN
            CALL N1LVLS(N,T1,XW,F,DX,CONV,SUMX,DLEVF,MPRMON,NEW.EQ.0)
          ELSE
            CALL N1LVLS(N,T1,XWA,F,DX,CONV,SUMX,DLEVF,MPRMON,NEW.EQ.0)
          ENDIF
          DO 32 L1=1,N
            XA(L1)=X(L1)
32        CONTINUE
          SUMXA = SUMX
          DLEVXA = DSQRT(SUMXA/DBLE(FLOAT(N)))
          CONVA = CONV
          DXANRM = WNORM(N,DX,XW)
C         --------------------------------------------------------
C         3.3 A - priori estimate of damping factor FC
          IF(NITER.NE.0.AND.NONLIN.NE.1.AND.NEW.EQ.0.AND.
     $       .NOT. QORDI)THEN
C           ------------------------------------------------------
C           3.3.1 Computation of the denominator of a-priori
C                 estimate
            FCDNM = ZERO
            DO 331 L1=1,N
              FCDNM=FCDNM+((DX(L1)-DXQA(L1))/XW(L1))**2
331         CONTINUE
            FCDNM = FCDNM*SUMX
C           ------------------------------------------------------
C           3.3.2 New damping factor
            IF(FCDNM.GT.FCNUMP*FCMIN2 .OR.
     $        (NONLIN.EQ.4 .AND. FCA**2*FCNUMP .LT. 4.0D0*FCDNM)) THEN
              DMYPRI = FCA*DSQRT(FCNUMP/FCDNM)
              FCPRI = DMIN1(DMYPRI,ONE)
              IF (NONLIN.EQ.4) FCPRI = DMIN1(HALF*DMYPRI,ONE)
            ELSE
              FCPRI = ONE
C$Test-begin
              DMYPRI = -1.0D0
C$Test-end
            ENDIF
C$Test-begin
            IF (MPRMON.GE.5) THEN
              WRITE(LUMON,33201) FCPRI, FC, FCA, DMYPRI, FCNUMP,
     $                           FCDNM
33201         FORMAT (/, ' +++ apriori estimate +++', /,
     $                ' FCPRI  = ', D18.10, '  FC     = ', D18.10, /,
     $                ' FCA    = ', D18.10, '  DMYPRI = ', D18.10, /,
     $                ' FCNUMP = ', D18.10, '  FCDNM  = ', D18.10, /,
     $                   ' ++++++++++++++++++++++++', /)
            ENDIF
C$Test-end 
            FC = DMAX1(FCPRI,FCMIN)
            IF (QBDAMP) THEN
              FCBH = FCA*FCBND
              IF (FC.GT.FCBH) THEN
                FC = FCBH
                IF (MPRMON.GE.4)
     $            WRITE(LUMON,*) ' *** incr. rest. act. (a prio) ***'
              ENDIF
              FCBH = FCA/FCBND
              IF (FC.LT.FCBH) THEN
                FC = FCBH
                IF (MPRMON.GE.4)
     $            WRITE(LUMON,*) ' *** decr. rest. act. (a prio) ***'
              ENDIF
            ENDIF
          ENDIF
CWEI
          IF (IORMON.GE.2) THEN
            SUMXA2=SUMXA1
            SUMXA1=SUMXA0
            SUMXA0=DLEVXA
            IF (SUMXA0.EQ.ZERO) SUMXA0=SMALL
C           Check convergence rates (linear and superlinear)
C           ICONV : Convergence indicator
C                   =0: No convergence indicated yet
C                   =1: Damping factor is 1.0d0
C                   =2: Superlinear convergence detected (alpha >=1.2)
C                   =3: Quadratic convergence detected (alpha > 1.8)
            FCMON = DMIN1(FC,FCMON)
            IF (FCMON.LT.ONE) THEN
              ICONV = 0
              ALPHAE = ZERO
            ENDIF
            IF (FCMON.EQ.ONE .AND. ICONV.EQ.0) ICONV=1
            IF (NITER.GE.1) THEN
              CLIN1 = CLIN0
              CLIN0 = SUMXA0/SUMXA1
            ENDIF
            IF (ICONV.GE.1.AND.NITER.GE.2) THEN
              ALPHAK = ALPHAE
              ALPHAE = ZERO
              IF (CLIN1.LE.0.95D0) ALPHAE = DLOG(CLIN0)/DLOG(CLIN1)
              IF (ALPHAK.NE.ZERO) ALPHAK =0.5D0*(ALPHAE+ALPHAK)
              ALPHAA = DMIN1(ALPHAK,ALPHAE)
              CALPHK = CALPHA
              CALPHA = ZERO
              IF (ALPHAE.NE.ZERO) CALPHA = SUMXA1/SUMXA2**ALPHAE
              SUMXTE = DSQRT(CALPHA*CALPHK)*SUMXA1**ALPHAK-SUMXA0
              IF (ALPHAA.GE.1.2D0 .AND. ICONV.EQ.1) ICONV = 2
              IF (ALPHAA.GT.1.8D0) ICONV = 3
              IF (MPRMON.GE.4)  WRITE(LUMON,32001) ICONV, ALPHAE, 
     $                            CALPHA, CLIN0, ALPHAK, SUMXTE
32001         FORMAT(' ** ICONV: ',I1,'  ALPHA: ',D9.2,
     $               '  CONST-ALPHA: ',D9.2,'  CONST-LIN: ',D9.2,' **',
     $               /,' **',11X,'ALPHA-POST: ',D9.2,' CHECK: ',D9.2,
     $               25X,'**')
              IF ( ICONV.GE.2 .AND. ALPHAA.LT.0.9D0 ) THEN
                 IF (IORMON.EQ.3) THEN
                   IERR = 4
                   GOTO 4299
                 ELSE
                   QMSTOP = .TRUE.
                 ENDIF 
              ENDIF
            ENDIF
          ENDIF
          FCMON = FC
C
C         --------------------------------------------------------
C         3.4 Save natural level for later computations of
C             corrector and print iterate
          FCNUMK = SUMX
          IF (MPRMON.GE.2) THEN
            IF (MPRTIM.NE.0) CALL MONON(5)
            CALL N1PRV1(DLEVF,DLEVXA,FCKEEP,NITER,NEW,MPRMON,LUMON,
     $                  QMIXIO)
            IF (MPRTIM.NE.0) CALL MONOFF(5)
          ENDIF
          NRED = 0
          QNEXT = .FALSE.
          QREP  = .FALSE.   
C         QREP = ITER .GT. ITMAX   or  QREP = ITER .GT. 0
C
C         Damping-factor reduction loop
C         ================================
C         DO (Until)
34        CONTINUE
C           ------------------------------------------------------
C           3.5 Preliminary new iterate
            DO 35 L1=1,N
              X(L1)=XA(L1)+DX(L1)*FC
35          CONTINUE
C           -----------------------------------------------------
C           3.5.2 Exit, if problem is specified as being linear
C     Exit Repeat If ...
            IF( NONLIN.EQ.1 )THEN
              IERR = 0
              GOTO 4299
            ENDIF
C           ------------------------------------------------------
C           3.6.1 Computation of the residual vector
            IF (MPRTIM.NE.0) CALL MONON(1)
            CALL FCN(N,X,F,IFAIL)
            IF (MPRTIM.NE.0) CALL MONOFF(1)
            NFCN = NFCN+1
C     Exit, If ...
            IF(IFAIL.LT.0)THEN
              IERR = 82
              GOTO 4299
            ENDIF
            IF(IFAIL.EQ.1 .OR. IFAIL.EQ.2) THEN
              IF (QORDI) THEN
                IERR = 82
                GOTO 4299
              ENDIF
              IF (IFAIL.EQ.1) THEN
                FCREDU = HALF
              ELSE
                FCREDU = F(1)
C     Exit, If ...
                IF (FCREDU.LE.0 .OR. FCREDU.GE.1) THEN
                  IERR = 82
                  GOTO 4299
                ENDIF
              ENDIF
              IF (MPRMON.GE.2) THEN
36101           FORMAT(8X,I2,' FCN could not be evaluated  ',
     $                 8X,F7.5,4X,I2)
                WRITE(LUMON,36101)NITER,FC,NEW
              ENDIF
              FCH = FC
              FC = FCREDU*FC
              IF (FCH.GT.FCMIN) FC = DMAX1(FC,FCMIN)
              IF (QBDAMP) THEN
                FCBH = FCH/FCBND
                IF (FC.LT.FCBH) THEN
                  FC = FCBH
                  IF (MPRMON.GE.4) WRITE(LUMON,*)
     $               ' *** decr. rest. act. (FCN redu.) ***'
                ENDIF
              ENDIF
              IF (FC.LT.FCMIN) THEN
                IERR = 3
                GOTO 4299
              ENDIF  
C     Break DO (Until) ...
              GOTO 3109
            ENDIF
            IF (QORDI) THEN
C             -------------------------------------------------------
C             3.6.2 Convergence test for ordinary Newton iteration
C     Exit Repeat If ...
              IF( DXANRM.LE.RTOL )THEN
                IERR = 0
                GOTO 4299
              ENDIF
            ELSE
              DO 361 L1=1,N
               T1(L1)=F(L1)*FW(L1)
361           CONTINUE
C             --------------------------------------------------
C             3.6.3 Solution of linear (N,N)-system
              IF (MPRTIM.NE.0) CALL MONON(4)
              CALL N1SOLV(N,M1,ML,MU,A,T1,IOPT,IFAIL,LIWL,IWK(NIWLA),
     $                   IDUMMY,LRWL,RWK(NRWLA),IDUMMY)
              IF (MPRTIM.NE.0) CALL MONOFF(4)
C     Exit Repeat If ...
              IF(IFAIL.NE.0)  THEN
                IERR = 81
                GOTO 4299
              ENDIF
              IF(NEW.GT.0) THEN 
                DO 3630 I=1,N
                  DXQ(I) = T1(I)*XWA(I)
3630            CONTINUE                   
                DO 363 ILOOP=1,NEW 
                  SUM1=ZERO
                  SUM2=ZERO
                  DO 3631 I=1,N
                    SUM1=SUM1+(DXQ(I)*DXSAVE(I,ILOOP))/ XW(I)**2
                    SUM2=SUM2+(DXSAVE(I,ILOOP)/XW(I))**2
3631              CONTINUE
                  BETA=SUM1/SUM2
                  DO 3632 I=1,N
                    DXQ(I)=DXQ(I)+BETA*DXSAVE(I,ILOOP+1)
                    T1(I) = DXQ(I)/XW(I)
3632              CONTINUE
363             CONTINUE
              ENDIF
C             ----------------------------------------------------
C             3.6.4 Evaluation of scaled natural level function
C                   SUMX
C                   scaled maximum error norm CONV and evaluation
C                   of (scaled) standard level function DLEVFN
              IF (.NOT. QSIMPL) THEN
                CALL N1LVLS(N,T1,XW,F,DXQ,CONV,SUMX,DLEVFN,MPRMON,
     $                      NEW.EQ.0)
              ELSE
                CALL N1LVLS(N,T1,XWA,F,DXQ,CONV,SUMX,DLEVFN,MPRMON,
     $                      NEW.EQ.0)
              ENDIF  
              DXNRM = WNORM(N,DXQ,XW)
C             -----------------------------------------------------
C             3.6.5 Convergence test
C     Exit Repeat If ...
              IF ( DXNRM.LE.RTOL .AND. DXANRM.LE.RSMALL .AND. 
     $            FC.EQ.ONE ) THEN
                IERR = 0
                GOTO 4299
              ENDIF
C           
              FCA = FC
C             ----------------------------------------------------
C             3.6.6 Evaluation of reduced damping factor
              TH = FCA-ONE
              FCDNM = ZERO
              DO 39 L1=1,N
                FCDNM=FCDNM+((DXQ(L1)+TH*DX(L1))/XW(L1))**2
39            CONTINUE
              IF (FCDNM.NE.ZERO) THEN
                DMYCOR = FCA*FCA*HALF*DSQRT(FCNUMK/FCDNM)
              ELSE
                DMYCOR = 1.0D+35
              ENDIF
              IF (NONLIN.LE.3) THEN
                FCCOR = DMIN1(ONE,DMYCOR)
              ELSE
                FCCOR = DMIN1(ONE,HALF*DMYCOR)
              ENDIF
C$Test-begin
              IF (MPRMON.GE.5) THEN
                WRITE(LUMON,39001) FCCOR, FC, DMYCOR, FCNUMK,
     $                             FCDNM, FCA
39001           FORMAT (/, ' +++ corrector computation +++', /,
     $                ' FCCOR  = ', D18.10, '  FC     = ', D18.10, /,
     $                ' DMYCOR = ', D18.10, '  FCNUMK = ', D18.10, /,
     $                ' FCDNM  = ', D18.10, '  FCA    = ', D18.10, /,
     $                   ' +++++++++++++++++++++++++++++', /)
              ENDIF
C$Test-end 
            ENDIF
C           ------------------------------------------------------
C           3.7 Natural monotonicity test
            IF(SUMX.GT.SUMXA .AND. .NOT.QORDI)THEN
C             ----------------------------------------------------
C             3.8 Output of iterate
              IF(MPRMON.GE.3) THEN
                IF (MPRTIM.NE.0) CALL MONON(5)
                CALL N1PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC,
     $                      NITER,MPRMON,LUMON,QMIXIO,'*')
                IF (MPRTIM.NE.0) CALL MONOFF(5)
              ENDIF
              IF (QMSTOP) THEN
                IERR = 4
                GOTO 4299
              ENDIF
              FCH = DMIN1(FCCOR,HALF*FC)
              IF (FC.GT.FCMIN) THEN
                FC=DMAX1(FCH,FCMIN)
              ELSE
                FC=FCH
              ENDIF
              IF (QBDAMP) THEN
                FCBH = FCA/FCBND
                IF (FC.LT.FCBH) THEN
                  FC = FCBH
                  IF (MPRMON.GE.4)
     $              WRITE(LUMON,*) ' *** decr. rest. act. (a post) ***'
                ENDIF
              ENDIF
CWEI
              FCMON = FC
C
C$Test-begin
                IF (MPRMON.GE.5) THEN
                  WRITE(LUMON,39002) FC
39002             FORMAT (/, ' +++ corrector setting 1 +++', /,
     $                    ' FC     = ', D18.10, /,
     $                       ' +++++++++++++++++++++++++++', /)
                ENDIF
C$Test-end 
              QREP = .TRUE.
              NCORR = NCORR+1
              NRED = NRED+1
C             ----------------------------------------------------
C             3.10 If damping factor is too small:
C                  Refresh Jacobian,if current Jacobian was computed
C                  by a Rank1-update, else fail exit
              QJCRFR  = FC.LT.FCMIN.OR.NEW.GT.0.AND.NRED.GT.1
C     Exit Repeat If ...
              IF(QJCRFR.AND.NEW.EQ.0)THEN
                IERR = 3
                GOTO 4299
              ENDIF
            ELSE
              IF (.NOT.QORDI.AND..NOT.QREP .AND. FCCOR.GT.SIGMA2*FC)
     $        THEN
                IF(MPRMON.GE.3) THEN
                  IF (MPRTIM.NE.0) CALL MONON(5)
                  CALL N1PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC,
     $                        NITER,MPRMON,LUMON,QMIXIO,'+')
                  IF (MPRTIM.NE.0) CALL MONOFF(5)
                ENDIF
                FC = FCCOR
C$Test-begin
                IF (MPRMON.GE.5) THEN
                  WRITE(LUMON,39003) FC
39003             FORMAT (/, ' +++ corrector setting 2 +++', /,
     $                    ' FC     = ', D18.10, /,
     $                       ' +++++++++++++++++++++++++++', /)
                ENDIF
C$Test-end 
                QREP = .TRUE.
              ELSE
                QNEXT = .TRUE.
              ENDIF
            ENDIF
3109      CONTINUE
          IF(.NOT.(QNEXT.OR.QJCRFR)) GOTO  34
C         UNTIL ( expression - negated above)
C         End of damping-factor reduction loop
C         =======================================
          IF(QJCRFR)THEN
C           ------------------------------------------------------
C           3.11 Restore former values for repeting iteration
C                step
            NREJR1 = NREJR1+1
            DO 3111 L1=1,N
              X(L1)=XA(L1)
3111        CONTINUE
            DO 3112 L1=1,N
              F(L1)=FA(L1)
3112        CONTINUE
            IF(MPRMON.GE.2)THEN
31130           FORMAT(8X,I2,' Not accepted damping factor ',
     $                 8X,F7.5,4X,I2)
                WRITE(LUMON,31130)NITER,FC,NEW
            ENDIF
            FC = FCKEEP
            FCA = FCK2
            IF(NITER.EQ.0)THEN
              FC = FCMIN
            ENDIF
            QGENJ = .TRUE.
          ELSE
C           ------------------------------------------------------
C           4 Preparations to start the following iteration step
C           ------------------------------------------------------
C           4.1 Print values
            IF(MPRMON.GE.3 .AND. .NOT.QORDI) THEN
              IF (MPRTIM.NE.0) CALL MONON(5)
              CALL N1PRV2(DLEVFN,DSQRT(SUMX/DBLE(FLOAT(N))),FC,NITER+1,
     $                    MPRMON,LUMON,QMIXIO,'*')
              IF (MPRTIM.NE.0) CALL MONOFF(5)
            ENDIF
C           Print the natural level of the current iterate and return
C           it in one-step mode
            SUMXS = SUMX
            SUMX = SUMXA
            IF(MPRSOL.GE.2.AND.NITER.NE.0) THEN
              IF (MPRTIM.NE.0) CALL MONON(5)
              CALL N1SOUT(N,XA,2,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL)
              IF (MPRTIM.NE.0) CALL MONOFF(5)
            ELSE IF(MPRSOL.GE.1.AND.NITER.EQ.0)THEN
              IF (MPRTIM.NE.0) CALL MONON(5)
              CALL N1SOUT(N,XA,1,IOPT,RWK,LRWK,IWK,LIWK,MPRSOL,LUSOL)
              IF (MPRTIM.NE.0) CALL MONOFF(5)
            ENDIF
            NITER = NITER+1
            DLEVF = DLEVFN
C     Exit Repeat If ...
            IF(NITER.GE.NITMAX)THEN
              IERR = 2
              GOTO 4299
            ENDIF
            FCKEEP = FC
C           ------------------------------------------------------
C           4.2 Return, if in one-step mode
C
C Exit Subroutine If ...
            IF (MODE.EQ.1) THEN
              IWK(18)=NIWLA-1
              IWK(19)=NRWLA-1
              IOPT(1)=1
              RETURN
            ENDIF
          ENDIF
        GOTO 2
C       End Repeat
4299    CONTINUE

This snippet took 0.05 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).