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) LOGICAL CALHES COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,55,55,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 RETURN C 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), & E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1) COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,55,55,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 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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,55,55,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) DO I=1,N DATA(2*I-1)=z1(I) DATA(2*I)=0.d0 END DO N2=2*NS CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N z1(I)=DATA(2*I-1) END DO 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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,55,55,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 C CALL SOLBC (N,LDE1,E2R,E2I,MLE,MUE,Z2,Z3,IP2) DO I=1,N DATA(2*I-1)=z2(I) DATA(2*I)=z3(I) END DO CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=ALPHN+4.0D0*ALF*NSSQ N2=2*NS DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERMR=(TI-TCOS(J))*NSSQ TERMI=BETAN*NSSQ TNORM=TERMR**2+TERMI**2 TERMR=TERMR/TNORM TERMI=-TERMI/TNORM IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATASA=DATA(IJ-1)*TERMR-DATA(IJ)*TERMI DATA(IJ)=DATA(IJ)*TERMR+DATA(IJ-1)*TERMI DATA(IJ-1)=DATASA DATASA=DATA(IJ2-1)*TERMR-DATA(IJ2)*TERMI DATA(IJ2)=DATA(IJ2)*TERMR+DATA(IJ2-1)*TERMI DATA(IJ2-1)=DATASA END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N z2(I)=DATA(2*I-1) z3(I)=DATA(2*I) END DO 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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG C GOTO (1,2,55,55,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) DO I=1,N DATA(2*I-1)=z1(I) DATA(2*I)=0.d0 END DO N2=2*NS CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N z1(I)=DATA(2*I-1) END DO DO I=1,N DATA(2*I-1)=z2(I) DATA(2*I)=z3(I) END DO CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=ALPHN+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERMR=(TI-TCOS(J))*NSSQ TERMI=BETAN*NSSQ TNORM=TERMR**2+TERMI**2 TERMR=TERMR/TNORM TERMI=-TERMI/TNORM IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATASA=DATA(IJ-1)*TERMR-DATA(IJ)*TERMI DATA(IJ)=DATA(IJ)*TERMR+DATA(IJ-1)*TERMI DATA(IJ-1)=DATASA DATASA=DATA(IJ2-1)*TERMR-DATA(IJ2)*TERMI DATA(IJ2)=DATA(IJ2)*TERMR+DATA(IJ2-1)*TERMI DATA(IJ2-1)=DATASA END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N z2(I)=DATA(2*I-1) z3(I)=DATA(2*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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG HEE1=DD1/H HEE2=DD2/H HEE3=DD3/H GOTO (1,2,55,55,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 DO I=1,N DATA(2*I-1)=CONT(I) DATA(2*I)=0.d0 END DO N2=2*NS CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N CONT(I)=DATA(2*I-1) 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 DO I=1,N DATA(2*I-1)=CONT(I) DATA(2*I)=0.d0 END DO N2=2*NS CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N CONT(I)=DATA(2*I-1) 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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NP,NSSQ COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG GOTO (1,2,55,55,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 DO I=1,N DATA(2*I-1)=CONT(I) DATA(2*I)=0.d0 END DO N2=2*NP CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NP TI=FAC1N-TCOS(I) II=I*2 DO J=1,NP TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N CONT(I)=DATA(2*I-1) 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,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 DO I=1,N DATA(2*I-1)=CONT(I) DATA(2*I)=0.d0 END DO N2=2*NP CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NP TI=FAC1N-TCOS(I) II=I*2 DO J=1,NP TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N CONT(I)=DATA(2*I-1) 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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ 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,55,55,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 DO I=1,N DATA(2*I-1)=AK(I) DATA(2*I)=0.d0 END DO N2=2*NS CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N AK(I)=DATA(2*I-1) 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/datafft/ data(256*256*4) COMMON /FOURIER/TCOS(512),NF(2),ALF,NDIM,NS,NSSQ 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 DO I=1,N DATA(2*I-1)=DEL(I) DATA(2*I)=0.d0 END DO N2=2*NS CALL FOURN (DATA,NF,NDIM,-1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,-1) FAC1N=FAC1+4.0D0*ALF*NSSQ DO I=1,NS TI=FAC1N-TCOS(I) II=I*2 DO J=1,NS TERM=(TI-TCOS(J))*NSSQ IJ=II+(J-1)*N2 IJ2=2*NSSQ+IJ DATA(IJ-1)=DATA(IJ-1)/TERM DATA(IJ)=DATA(IJ)/TERM DATA(IJ2-1)=DATA(IJ2-1)/TERM DATA(IJ2)=DATA(IJ2)/TERM END DO END DO CALL FOURN (DATA,NF,NDIM,1) CALL FOURN (DATA(2*NSSQ+1),NF,NDIM,1) DO I=1,N DEL(I)=DATA(2*I-1) END DO RETURN C C ----------------------------------------------------------- C 55 CONTINUE RETURN END C C END OF SUBROUTINE SLVSEU C C SUBROUTINE FOURN (DATA,NN,NDIM,ISIGN) C --- explication in "Numerical Recipies", page 451 IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION NN(NDIM),DATA(*) NTOT=1 DO IDIM=1,NDIM NTOT=NTOT*NN(IDIM) END DO NPREV=1 DO IDIM=1,NDIM N=NN(IDIM) NREM=NTOT/(N*NPREV) IP1=2*NPREV IP2=IP1*N IP3=IP2*NREM I2REV=1 DO I2=1,IP2,IP1 IF (I2.LT.I2REV) THEN DO I1=I2,I2+IP1-2,2 DO I3=I1,IP3,IP2 I3REV=I2REV+I3-I2 TEMPR=DATA(I3) TEMPI=DATA(I3+1) DATA(I3)=DATA(I3REV) DATA(I3+1)=DATA(I3REV+1) DATA(I3REV)=TEMPR DATA(I3REV+1)=TEMPI END DO END DO END IF IBIT=IP2/2 1 IF ((IBIT.GE.IP1).AND.(I2REV.GT.IBIT)) THEN I2REV=I2REV-IBIT IBIT=IBIT/2 GO TO 1 END IF I2REV=I2REV+IBIT END DO IFP1=IP1 2 IF (IFP1.LT.IP2) THEN IFP2=2*IFP1 THETA=ISIGN*6.2831853071795864769252D0/(IFP2/IP1) WPR=-2.D0*SIN(0.5D0*THETA)**2 WPI=SIN(THETA) WR=1.0D0 WI=0.0D0 DO I3=1,IFP1,IP1 DO I1=I3,I3+IP1-2,2 DO I2=I1,IP3,IFP2 K1=I2 K2=K1+IFP1 TEMPR=WR*DATA(K2)-WI*DATA(K2+1) TEMPI=WR*DATA(K2+1)+WI*DATA(K2) DATA(K2)=DATA(K1)-TEMPR DATA(K2+1)=DATA(K1+1)-TEMPI DATA(K1)=DATA(K1)+TEMPR DATA(K1+1)=DATA(K1+1)+TEMPI END DO END DO WTEMP=WR WR=WR*WPR-WI*WPI+WR WI=WI*WPR+WTEMP*WPI+WI END DO IFP1=IFP2 GO TO 2 END IF NPREV=N*NPREV END DO RETURN END