C ****************************************** C VERSION OF JANUARY 26, 1996 C ****************************************** C SUBROUTINE DECOMR(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,FAC1,E1,LDE1,IP1,IER,IJOB,CALHES,IPHES) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1), & IP1(NM1),IPHES(N),HELP(5000) LOGICAL CALHES COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX DO J=1,N DO I=1,N E1(I,J)=-FJAC(I,J) END DO E1(J,J)=E1(J,J)+FAC1 END DO CALL DEC (N,LDE1,E1,IP1,IER) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX C DO J=1,N C DO I=1,MBJAC C E1(I+MLE,J)=-FJAC(I,J) C END DO C E1(MDIAG,J)=E1(MDIAG,J)+FAC1 C END DO C CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER) DO J=2,N DO I=1,3 E1(I+1,J)=-FJAC(I,J) END DO E1(3,J)=E1(3,J)+FAC1 END DO CALL DECB (N-1,LDE1,E1(1,2),1,1,IP1,IER) DO I=1,N-1 HELP(I)=-FJAC(5,I+1) END DO CALL SOLB (N-1,LDE1,E1(1,2),1,1,HELP,IP1) SUM=0.0D0 DO I=1,N-1 E1(5,I+1)=HELP(I) SUM=SUM+FJAC(4,I+1)*HELP(I) END DO E1(5,1)=FAC1-FJAC(2,1)+SUM RETURN C C C ----------------------------------------------------------- C 3 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO J=1,N DO I=1,N E1(I,J)=-FJAC(I,J) END DO DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS) E1(I,J)=E1(I,J)+FAC1*FMAS(I-J+MBDIAG,J) END DO END DO CALL DEC (N,LDE1,E1,IP1,IER) RETURN C C ----------------------------------------------------------- C 4 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX C DO J=1,N C DO I=1,MBJAC C E1(I+MLE,J)=-FJAC(I,J) C END DO C DO I=1,MBB C IB=I+MDIFF C E1(IB,J)=E1(IB,J)+FAC1*FMAS(I,J) C END DO C END DO C CALL DECB (N,LDE1,E1,MLE,MUE,IP1,IER) DO J=2,N DO I=1,3 E1(I+1,J)=-FJAC(I,J) END DO E1(3,J)=E1(3,J)+FAC1*FMAS(1,J) END DO CALL DECB (N-1,LDE1,E1(1,2),1,1,IP1,IER) DO I=1,N-1 HELP(I)=-FJAC(5,I+1) END DO CALL SOLB (N-1,LDE1,E1(1,2),1,1,HELP,IP1) SUM=0.0D0 DO I=1,N-1 E1(5,I+1)=HELP(I) SUM=SUM+FJAC(4,I+1)*HELP(I) END DO E1(5,1)=FAC1*FMAS(1,1)-FJAC(2,1)+SUM RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE DECOMR C C *********************************************************** C SUBROUTINE DECOMC(N,FJAC,LDJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,IP2,IER,IJOB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),HELP1(5000),HELP2(5000), & E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1) COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX DO J=1,N DO I=1,N E2R(I,J)=-FJAC(I,J) E2I(I,J)=0.D0 END DO E2R(J,J)=E2R(J,J)+ALPHN E2I(J,J)=BETAN END DO CALL DECC (N,LDE1,E2R,E2I,IP2,IER) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX C DO J=1,N C DO I=1,MBJAC C IMLE=I+MLE C E2R(IMLE,J)=-FJAC(I,J) C E2I(IMLE,J)=0.D0 C END DO C E2R(MDIAG,J)=E2R(MDIAG,J)+ALPHN C E2I(MDIAG,J)=BETAN C END DO C CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER) DO J=2,N DO I=1,3 IMLE=I+1 E2R(IMLE,J)=-FJAC(I,J) E2I(IMLE,J)=0.D0 END DO E2R(3,J)=E2R(3,J)+ALPHN E2I(3,J)=BETAN END DO CALL DECBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,IP2,IER) DO I=1,N-1 HELP1(I)=-FJAC(5,I+1) HELP2(I)=0.0D0 END DO CALL SOLBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,HELP1,HELP2,IP2) SUMR=0.0D0 SUMI=0.0D0 DO I=1,N-1 E2R(5,I+1)=HELP1(I) E2I(5,I+1)=HELP2(I) SUMR=SUMR+FJAC(4,I+1)*HELP1(I) SUMI=SUMI+FJAC(4,I+1)*HELP2(I) END DO E2R(5,1)=ALPHN-FJAC(2,1)+SUMR E2I(5,1)=BETAN+SUMI RETURN C C ----------------------------------------------------------- C 3 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO J=1,N DO I=1,N E2R(I,J)=-FJAC(I,J) E2I(I,J)=0.D0 END DO END DO DO J=1,N DO I=MAX(1,J-MUMAS),MIN(N,J+MLMAS) BB=FMAS(I-J+MBDIAG,J) E2R(I,J)=E2R(I,J)+ALPHN*BB E2I(I,J)=BETAN*BB END DO END DO CALL DECC(N,LDE1,E2R,E2I,IP2,IER) RETURN C C ----------------------------------------------------------- C 4 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX C DO J=1,N C DO I=1,MBJAC C IMLE=I+MLE C E2R(IMLE,J)=-FJAC(I,J) C E2I(IMLE,J)=0.D0 C END DO C DO I=MAX(1,MUMAS+2-J),MIN(MBB,MUMAS+1-J+N) C IB=I+MDIFF C BB=FMAS(I,J) C E2R(IB,J)=E2R(IB,J)+ALPHN*BB C E2I(IB,J)=BETAN*BB C END DO C END DO C CALL DECBC (N,LDE1,E2R,E2I,MLE,MUE,IP2,IER) DO J=2,N DO I=1,3 IMLE=I+1 E2R(IMLE,J)=-FJAC(I,J) E2I(IMLE,J)=0.D0 END DO BB=FMAS(1,J) E2R(3,J)=E2R(3,J)+ALPHN*BB E2I(3,J)=BETAN*BB END DO CALL DECBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,IP2,IER) DO I=1,N-1 HELP1(I)=-FJAC(5,I+1) HELP2(I)=0.0D0 END DO CALL SOLBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,HELP1,HELP2,IP2) SUMR=0.0D0 SUMI=0.0D0 DO I=1,N-1 E2R(5,I+1)=HELP1(I) E2I(5,I+1)=HELP2(I) SUMR=SUMR+FJAC(4,I+1)*HELP1(I) SUMI=SUMI+FJAC(4,I+1)*HELP2(I) END DO BB=FMAS(1,1) E2R(5,1)=ALPHN*BB-FJAC(2,1)+SUMR E2I(5,1)=BETAN*BB+SUMI RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE DECOMC C C *********************************************************** C SUBROUTINE SLVRAR(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,FAC1,E1,LDE1,Z1,F1,IP1,IPHES,IER,IJOB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1), & IP1(NM1),IPHES(N),Z1(N),F1(N) COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX DO I=1,N Z1(I)=Z1(I)-F1(I)*FAC1 END DO CALL SOL (N,LDE1,E1,Z1,IP1) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX DO I=1,N Z1(I)=Z1(I)-F1(I)*FAC1 END DO C CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) CALL SOLB (N-1,LDE1,E1(1,2),1,1,Z1(2),IP1) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*Z1(I) END DO Z1(1)=(Z1(1)+SUM)/E1(5,1) DO I=2,N Z1(I)=Z1(I)-Z1(1)*E1(5,I) END DO RETURN C C ----------------------------------------------------------- C 3 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO I=1,N S1=0.0D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) S1=S1-FMAS(I-J+MBDIAG,J)*F1(J) END DO Z1(I)=Z1(I)+S1*FAC1 END DO CALL SOL (N,LDE1,E1,Z1,IP1) RETURN C C ----------------------------------------------------------- C 4 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX DO I=1,N S1=0.0D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) S1=S1-FMAS(I-J+MBDIAG,J)*F1(J) END DO Z1(I)=Z1(I)+S1*FAC1 END DO CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE SLVRAR C C *********************************************************** C SUBROUTINE SLVRAI(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,ALPHN,BETAN,E2R,E2I,LDE1,Z2,Z3, & F2,F3,CONT,IP2,IPHES,IER,IJOB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1), & IP2(NM1),IPHES(N),Z2(N),Z3(N),F2(N),F3(N) DIMENSION E2R(LDE1,NM1),E2I(LDE1,NM1) COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX DO I=1,N S2=-F2(I) S3=-F3(I) Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX DO I=1,N S2=-F2(I) S3=-F3(I) Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOLBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,Z2(2),Z3(2),IP2) SUMR=0.0D0 SUMI=0.0D0 DO I=2,N FJAC4=FJAC(4,I) SUMR=SUMR+FJAC4*Z2(I) SUMI=SUMI+FJAC4*Z3(I) END DO CREAL=E2R(5,1) CIMAG=E2I(5,1) CNORM=CREAL**2+CIMAG**2 CREAL=CREAL/CNORM CIMAG=-CIMAG/CNORM Z2A=Z2(1)+SUMR Z3A=Z3(1)+SUMI Z2(1)=Z2A*CREAL-Z3A*CIMAG Z3(1)=Z3A*CREAL+Z2A*CIMAG DO I=2,N Z2(I)=Z2(I)-(Z2(1)*E2R(5,I)-Z3(1)*E2I(5,I)) Z3(I)=Z3(I)-(Z2(1)*E2I(5,I)+Z3(1)*E2R(5,I)) END DO RETURN C C ----------------------------------------------------------- C 3 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO I=1,N S2=0.0D0 S3=0.0D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) BB=FMAS(I-J+MBDIAG,J) S2=S2-BB*F2(J) S3=S3-BB*F3(J) END DO Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2) RETURN C C ----------------------------------------------------------- C 4 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX DO I=1,N S2=0.0D0 S3=0.0D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) BB=FMAS(I-J+MBDIAG,J) S2=S2-BB*F2(J) S3=S3-BB*F3(J) END DO Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOLBC(N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE SLVRAI C C *********************************************************** C SUBROUTINE SLVRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,FAC1,ALPHN,BETAN,E1,E2R,E2I,LDE1,Z1,Z2,Z3, & F1,F2,F3,CONT,IP1,IP2,IPHES,IER,IJOB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1), & E2R(LDE1,NM1),E2I(LDE1,NM1),IP1(NM1),IP2(NM1), & IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),F3(N) COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX DO I=1,N S2=-F2(I) S3=-F3(I) Z1(I)=Z1(I)-F1(I)*FAC1 Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOL (N,LDE1,E1,Z1,IP1) CALL SOLC (N,LDE1,E2R,E2I,Z2,Z3,IP2) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX DO I=1,N S2=-F2(I) S3=-F3(I) Z1(I)=Z1(I)-F1(I)*FAC1 Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO C CALL SOLB (N,LDE1,E1,MLE,MUE,Z1,IP1) C CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) CALL SOLB (N-1,LDE1,E1(1,2),1,1,Z1(2),IP1) CALL SOLBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,Z2(2),Z3(2),IP2) SUM=0.0D0 SUMR=0.0D0 SUMI=0.0D0 DO I=2,N FJAC4=FJAC(4,I) SUM=SUM+FJAC4*Z1(I) SUMR=SUMR+FJAC4*Z2(I) SUMI=SUMI+FJAC4*Z3(I) END DO Z1(1)=(Z1(1)+SUM)/E1(5,1) CREAL=E2R(5,1) CIMAG=E2I(5,1) CNORM=CREAL**2+CIMAG**2 CREAL=CREAL/CNORM CIMAG=-CIMAG/CNORM Z2A=Z2(1)+SUMR Z3A=Z3(1)+SUMI Z2(1)=Z2A*CREAL-Z3A*CIMAG Z3(1)=Z3A*CREAL+Z2A*CIMAG DO I=2,N Z1(I)=Z1(I)-Z1(1)*E1(5,I) Z2(I)=Z2(I)-(Z2(1)*E2R(5,I)-Z3(1)*E2I(5,I)) Z3(I)=Z3(I)-(Z2(1)*E2I(5,I)+Z3(1)*E2R(5,I)) END DO RETURN C C ----------------------------------------------------------- C 3 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO I=1,N S1=0.0D0 S2=0.0D0 S3=0.0D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) BB=FMAS(I-J+MBDIAG,J) S1=S1-BB*F1(J) S2=S2-BB*F2(J) S3=S3-BB*F3(J) END DO Z1(I)=Z1(I)+S1*FAC1 Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOL (N,LDE1,E1,Z1,IP1) CALL SOLC(N,LDE1,E2R,E2I,Z2,Z3,IP2) RETURN C C ----------------------------------------------------------- C 4 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX DO I=1,N BB=FMAS(1,I) S1=-BB*F1(I) S2=-BB*F2(I) S3=-BB*F3(I) Z1(I)=Z1(I)+S1*FAC1 Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN END DO CALL SOLB (N-1,LDE1,E1(1,2),1,1,Z1(2),IP1) CALL SOLBC (N-1,LDE1,E2R(1,2),E2I(1,2),1,1,Z2(2),Z3(2),IP2) SUM=0.0D0 SUMR=0.0D0 SUMI=0.0D0 DO I=2,N FJAC4=FJAC(4,I) SUM=SUM+FJAC4*Z1(I) SUMR=SUMR+FJAC4*Z2(I) SUMI=SUMI+FJAC4*Z3(I) END DO Z1(1)=(Z1(1)+SUM)/E1(5,1) CREAL=E2R(5,1) CIMAG=E2I(5,1) CNORM=CREAL**2+CIMAG**2 CREAL=CREAL/CNORM CIMAG=-CIMAG/CNORM Z2A=Z2(1)+SUMR Z3A=Z3(1)+SUMI Z2(1)=Z2A*CREAL-Z3A*CIMAG Z3(1)=Z3A*CREAL+Z2A*CIMAG DO I=2,N Z1(I)=Z1(I)-Z1(1)*E1(5,I) Z2(I)=Z2(I)-(Z2(1)*E2R(5,I)-Z3(1)*E2I(5,I)) Z3(I)=Z3(I)-(Z2(1)*E2I(5,I)+Z3(1)*E2R(5,I)) END DO RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE SLVRAD C C *********************************************************** C SUBROUTINE ESTRAD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & H,DD1,DD2,DD3,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1, & E1,LDE1,Z1,Z2,Z3,CONT,F1,F2,IP1,IPHES,SCAL,ERR, & FIRST,REJECT,FAC1,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1), & SCAL(N),IPHES(N),Z1(N),Z2(N),Z3(N),F1(N),F2(N),Y0(N),Y(N) DIMENSION CONT(N),RPAR(1),IPAR(1) LOGICAL FIRST,REJECT COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG HEE1=DD1/H HEE2=DD2/H HEE3=DD3/H GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C 1 CONTINUE C ------ B=IDENTITY, JACOBIAN A FULL MATRIX DO I=1,N F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) CONT(I)=F2(I)+Y0(I) END DO CALL SOL (N,LDE1,E1,CONT,IP1) GOTO 77 C 2 CONTINUE C ------ B=IDENTITY, JACOBIAN A BANDED MATRIX DO I=1,N F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) CONT(I)=F2(I)+Y0(I) END DO C CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) CALL SOLB (N-1,LDE1,E1(1,2),1,1,CONT(2),IP1) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*CONT(I) END DO CONT(1)=(CONT(1)+SUM)/E1(5,1) DO I=2,N CONT(I)=CONT(I)-CONT(1)*E1(5,I) END DO GOTO 77 C 3 CONTINUE C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO I=1,N F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) END DO DO I=1,N SUM=0.D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) SUM=SUM+FMAS(I-J+MBDIAG,J)*F1(J) END DO F2(I)=SUM CONT(I)=SUM+Y0(I) END DO CALL SOL (N,LDE1,E1,CONT,IP1) GOTO 77 C 4 CONTINUE C ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX DO I=1,N F1(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I) END DO DO I=1,N SUM=FMAS(1,I)*F1(I) F2(I)=SUM CONT(I)=SUM+Y0(I) END DO CALL SOLB (N-1,LDE1,E1(1,2),1,1,CONT(2),IP1) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*CONT(I) END DO CONT(1)=(CONT(1)+SUM)/E1(5,1) DO I=2,N CONT(I)=CONT(I)-CONT(1)*E1(5,I) END DO GOTO 77 C C -------------------------------------- C 77 CONTINUE ERR=0.D0 DO I=1,N ERR=ERR+(CONT(I)/SCAL(I))**2 END DO ERR=MAX(SQRT(ERR/N),1.D-10) C IF (ERR.LT.1.D0) RETURN IF (FIRST.OR.REJECT) THEN DO I=1,N CONT(I)=Y(I)+CONT(I) END DO CALL FCN(N,X,CONT,F1,RPAR,IPAR) NFCN=NFCN+1 DO I=1,N CONT(I)=F1(I)+F2(I) END DO GOTO (31,32,31,32,31,32,55,55,55,55,31,32,31,32,31), IJOB C ------ FULL MATRIX OPTION 31 CONTINUE CALL SOL(N,LDE1,E1,CONT,IP1) GOTO 88 C ------ BANDED MATRIX OPTION 32 CONTINUE C CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) CALL SOLB (N-1,LDE1,E1(1,2),1,1,CONT(2),IP1) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*CONT(I) END DO CONT(1)=(CONT(1)+SUM)/E1(5,1) DO I=2,N CONT(I)=CONT(I)-CONT(1)*E1(5,I) END DO GOTO 88 C ----------------------------------- 88 CONTINUE ERR=0.D0 DO I=1,N ERR=ERR+(CONT(I)/SCAL(I))**2 END DO ERR=MAX(SQRT(ERR/N),1.D-10) END IF RETURN C ----------------------------------------------------------- 55 CONTINUE RETURN END C C END OF SUBROUTINE ESTRAD C C *********************************************************** C SUBROUTINE ESTRAV(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & H,DD,FCN,NFCN,Y0,Y,IJOB,X,M1,M2,NM1,NS,NNS, & E1,LDE1,ZZ,CONT,FF,IP1,IPHES,SCAL,ERR, & FIRST,REJECT,FAC1,RPAR,IPAR) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E1(LDE1,NM1),IP1(NM1), & SCAL(N),IPHES(N),ZZ(NNS),FF(NNS),Y0(N),Y(N) DIMENSION DD(NS),CONT(N),RPAR(1),IPAR(1) LOGICAL FIRST,REJECT COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C 1 CONTINUE C ------ B=IDENTITY, JACOBIAN A FULL MATRIX DO I=1,N SUM=0.D0 DO K=1,NS SUM=SUM+DD(K)*ZZ(I+(K-1)*N) END DO FF(I+N)=SUM/H CONT(I)=FF(I+N)+Y0(I) END DO CALL SOL (N,LDE1,E1,CONT,IP1) GOTO 77 C 2 CONTINUE C ------ B=IDENTITY, JACOBIAN A BANDED MATRIX DO I=1,N SUM=0.D0 DO K=1,NS SUM=SUM+DD(K)*ZZ(I+(K-1)*N) END DO FF(I+N)=SUM/H CONT(I)=FF(I+N)+Y0(I) END DO CALL SOLB (N-1,LDE1,E1(1,2),1,1,CONT(2),IP1) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*CONT(I) END DO CONT(1)=(CONT(1)+SUM)/E1(5,1) DO I=2,N CONT(I)=CONT(I)-CONT(1)*E1(5,I) END DO GOTO 77 C 3 CONTINUE C ------ B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX DO I=1,N SUM=0.D0 DO K=1,NS SUM=SUM+DD(K)*ZZ(I+(K-1)*N) END DO FF(I)=SUM/H END DO DO I=1,N SUM=0.D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J) END DO FF(I+N)=SUM CONT(I)=SUM+Y0(I) END DO CALL SOL (N,LDE1,E1,CONT,IP1) GOTO 77 C 4 CONTINUE C ------ B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX DO I=1,N SUM=0.D0 DO K=1,NS SUM=SUM+DD(K)*ZZ(I+(K-1)*N) END DO FF(I)=SUM/H END DO DO I=1,N SUM=0.D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) SUM=SUM+FMAS(I-J+MBDIAG,J)*FF(J) END DO FF(I+N)=SUM CONT(I)=SUM+Y0(I) END DO CALL SOLB (N,LDE1,E1,MLE,MUE,CONT,IP1) GOTO 77 C C -------------------------------------- C 77 CONTINUE ERR=0.D0 DO I=1,N ERR=ERR+(CONT(I)/SCAL(I))**2 END DO ERR=MAX(SQRT(ERR/N),1.D-10) C IF (ERR.LT.1.D0) RETURN IF (FIRST.OR.REJECT) THEN DO I=1,N CONT(I)=Y(I)+CONT(I) END DO CALL FCN(N,X,CONT,FF,RPAR,IPAR) NFCN=NFCN+1 DO I=1,N CONT(I)=FF(I)+FF(I+N) END DO GOTO (31,32,31,32,31,32,55,55,55,55,31,32,31,32,31), IJOB C ------ FULL MATRIX OPTION 31 CONTINUE CALL SOL (N,LDE1,E1,CONT,IP1) GOTO 88 C ------ BANDED MATRIX OPTION 32 CONTINUE CALL SOLB (N-1,LDE1,E1(1,2),1,1,CONT(2),IP1) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*CONT(I) END DO CONT(1)=(CONT(1)+SUM)/E1(5,1) DO I=2,N CONT(I)=CONT(I)-CONT(1)*E1(5,I) END DO GOTO 88 C ----------------------------------- 88 CONTINUE ERR=0.D0 DO I=1,N ERR=ERR+(CONT(I)/SCAL(I))**2 END DO ERR=MAX(SQRT(ERR/N),1.D-10) END IF RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE ESTRAV C C *********************************************************** C SUBROUTINE SLVROD(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,FAC1,E,LDE,IP,DY,AK,FX,YNEW,HD,IJOB,STAGE1) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1), & IP(NM1),DY(N),AK(N),FX(N),YNEW(N) LOGICAL STAGE1 COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C IF (HD.EQ.0.D0) THEN DO I=1,N AK(I)=DY(I) END DO ELSE DO I=1,N AK(I)=DY(I)+HD*FX(I) END DO END IF C GOTO (1,2,3,4,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX IF (STAGE1) THEN DO I=1,N AK(I)=AK(I)+YNEW(I) END DO END IF CALL SOL (N,LDE,E,AK,IP) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX C IF (STAGE1) THEN C DO I=1,N C AK(I)=AK(I)+YNEW(I) C END DO C END IF C CALL SOLB (N,LDE,E,MLE,MUE,AK,IP) IF (STAGE1) THEN DO I=1,N AK(I)=AK(I)+YNEW(I) END DO END IF CALL SOLB (N-1,LDE,E(1,2),1,1,AK(2),IP) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*AK(I) END DO AK(1)=(AK(1)+SUM)/E(5,1) DO I=2,N AK(I)=AK(I)-AK(1)*E(5,I) END DO RETURN C C ----------------------------------------------------------- C 3 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A FULL MATRIX IF (STAGE1) THEN DO I=1,N SUM=0.D0 DO J=MAX(1,I-MLMAS),MIN(N,I+MUMAS) SUM=SUM+FMAS(I-J+MBDIAG,J)*YNEW(J) END DO AK(I)=AK(I)+SUM END DO END IF CALL SOL (N,LDE,E,AK,IP) RETURN C C ----------------------------------------------------------- C 4 CONTINUE C --- B IS A BANDED MATRIX, JACOBIAN A BANDED MATRIX IF (STAGE1) THEN DO I=1,N AK(I)=AK(I)+FMAS(1,I)*YNEW(I) END DO END IF CALL SOLB (N-1,LDE,E(1,2),1,1,AK(2),IP) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*AK(I) END DO AK(1)=(AK(1)+SUM)/E(5,1) DO I=2,N AK(I)=AK(I)-AK(1)*E(5,I) END DO RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE SLVROD C C C *********************************************************** C SUBROUTINE SLVSEU(N,FJAC,LDJAC,MLJAC,MUJAC,FMAS,LDMAS,MLMAS,MUMAS, & M1,M2,NM1,FAC1,E,LDE,IP,IPHES,DEL,IJOB) IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION FJAC(LDJAC,N),FMAS(LDMAS,NM1),E(LDE,NM1),DEL(N) DIMENSION IP(NM1),IPHES(N) COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,1,2,55,55,55,55,55,55,55,55,55,55,55), IJOB C C ----------------------------------------------------------- C 1 CONTINUE C --- B=IDENTITY, JACOBIAN A FULL MATRIX CALL SOL (N,LDE,E,DEL,IP) RETURN C C ----------------------------------------------------------- C 2 CONTINUE C --- B=IDENTITY, JACOBIAN A BANDED MATRIX C CALL SOLB (N,LDE,E,MLE,MUE,DEL,IP) CALL SOLB (N-1,LDE,E(1,2),1,1,DEL(2),IP) SUM=0.0D0 DO I=2,N SUM=SUM+FJAC(4,I)*DEL(I) END DO DEL(1)=(DEL(1)+SUM)/E(5,1) DO I=2,N DEL(I)=DEL(I)-DEL(1)*E(5,I) END DO RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE SLVSEU C