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
