C ******************************************
C     VERSION OF March 5, 2024
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,3,2,55,55,55,55,55,55,2,2,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
   3  CONTINUE
C ---  B IS A DIAGONAL MATRIX, 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*FMAS(MBDIAG,J)
      END DO
      CALL DEC (N,LDE1,E1,IP1,IER)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B = IDENTITY OR DIAGONAL, SUMEXP LINEAR ALGEBRA
      ND=MUE+1
      TERM = FAC1
      DO J=1,ND
         DO  I=1,ND
            E1(I,J)=-FJAC(I,J)
         END DO
         IF (IJOB.EQ.4) TERM = FAC1*FMAS(MBDIAG,J)
         E1(J,J)=E1(J,J)+TERM
      END DO
      SUM=0.0D0
      LEXP=N-ND-2
      II=ND+2
      DO I=1,LEXP
         FACT=FJAC(3,II)
         II=II+1
         IF (I.EQ.1.OR.FACT.EQ.0.0D0) THEN
            CC=1.0D0
         ELSE
            CC=FACT*U1
         END IF
         U1=CC/(FAC1-FJAC(2,II))
         SUM=SUM+FJAC(1,II)*U1
      END DO
      DO J=1,ND
         DO  I=1,ND
            E1(I,J)=E1(I,J)-SUM*FJAC(I,ND+1)*FJAC(J,ND+2)
         END DO
      END DO
      CALL DEC (ND,LDE1,E1,IP1,IER)
      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),
     &          E2R(LDE1,NM1),E2I(LDE1,NM1),IP2(NM1)
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,3,2,55,55,55,55,55,55,2,2,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
   3  CONTINUE
C ---  B IS A DIAGONAL 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
         BB=FMAS(MBDIAG,J)
         E2R(J,J)=E2R(J,J)+ALPHN*BB
         E2I(J,J)=BETAN*BB
      END DO
      CALL DECC(N,LDE1,E2R,E2I,IP2,IER)
      RETURN
C
C -----------------------------------------------------------
C
   2  CONTINUE
C ---  B = IDENTITY OR DIAGONAL, SUMEXP LINEAR ALGEBRA
      ND=MUE+1
      BB=1.0D0
      DO J=1,ND
         DO  I=1,ND
            E2R(I,J)=-FJAC(I,J)
            E2I(I,J)=0.D0
         END DO
         IF (IJOB.EQ.4) BB = FMAS(MBDIAG,J)
         E2R(J,J)=E2R(J,J)+ALPHN*BB
         E2I(J,J)=BETAN*BB
      END DO
      LEXP=N-ND-2
      SUMR=0.0D0
      SUMI=0.0D0
      II=ND+2
      DO I=1,LEXP
         II=II+1
         FACT=FJAC(3,II-1)
         IF (I.EQ.1.OR.FACT.EQ.0.0D0) THEN
            ALGAM=ALPHN-FJAC(2,II)
            DENOM=ALGAM**2+BETAN**2
            UR=ALGAM/DENOM
            UI=-BETAN/DENOM
         ELSE
            USAVE=(UR*ALGAM+UI*BETAN)/DENOM
            UI=FACT*(UI*ALGAM-UR*BETAN)/DENOM
            UR=FACT*USAVE
         END IF
         SUMR=SUMR+FJAC(1,II)*UR
         SUMI=SUMI+FJAC(1,II)*UI
      END DO
      DO J=1,ND
         DO  I=1,ND
            PROD=FJAC(I,ND+1)*FJAC(J,ND+2)
            E2R(I,J)=E2R(I,J)-SUMR*PROD
            E2I(I,J)=E2I(I,J)-SUMI*PROD
         END DO
      END DO
      CALL DECC (ND,LDE1,E2R,E2I,IP2,IER)
      RETURN
C
C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE DECOMC
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)
      DIMENSION CONT(4*N)
      COMMON/LINAL/MLE,MUE,MBJAC,MBB,MDIAG,MDIFF,MBDIAG
C
      GOTO (1,2,3,2,55,55,55,55,55,55,2,2,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
   3  CONTINUE
C ---  B IS A DIAGONAL MATRIX, JACOBIAN A FULL MATRIX
      DO I=1,N
         BB=FMAS(MBDIAG,I)
         S2=-BB*F2(I)
         S3=-BB*F3(I)
         Z1(I)=Z1(I)-BB*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, SUMEXP LINEAR ALGEBRA
      BB=1.0D0
      DO I=1,N
         IF (IJOB.EQ.4) BB=FMAS(MBDIAG,I)
         S2=-BB*F2(I)
         S3=-BB*F3(I)
         Z1(I)=Z1(I)-BB*F1(I)*FAC1
         Z2(I)=Z2(I)+S2*ALPHN-S3*BETAN
         Z3(I)=Z3(I)+S3*ALPHN+S2*BETAN
      END DO
      CALL SOLEXP(N,FJAC,LDJAC,MUE,NM1,LDE1,E1,Z1,IP1,FAC1,IJOB)
      CALL SOLEXPC(N,FJAC,LDJAC,MUE,NM1,LDE1,E2R,E2I,Z2,Z3,
     &             IP1,IP2,ALPHN,BETAN,IJOB)
C
      RETURN

C -----------------------------------------------------------
C
  55  CONTINUE
      RETURN
      END
C
C     END OF SUBROUTINE SLVRAD
C
C ***********************************************************
C
      SUBROUTINE SOLEXP(N,FJAC,LDJAC,MUJAC,NM1,LDE1,E1,Z1,IP1,
     &                  FAC1,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),E1(LDE1,NM1),IP1(NM1),Z1(N)
      ND=MUJAC+1
      LEXP=N-ND-2
      SUM=0.0D0
      II=ND+1
      DO I=1,LEXP  
         FACT=FJAC(3,II+1)
         IF (I.EQ.1.OR.FACT.EQ.0.0D0) THEN
            FJINV=1.0D0/(FAC1-FJAC(2,II+2))
            U1=Z1(II)*FJINV
         ELSE
            U1=(Z1(II)+FACT*U1)*FJINV
         END IF
         SUM=SUM+FJAC(1,II+2)*U1
         II=II+1
      END DO
      ND1=ND+1
      DO I=1,ND
         Z1(I)=Z1(I)+SUM*FJAC(I,ND1)
      END DO
      CALL SOL (ND,LDE1,E1,Z1,IP1)
      ND2=ND+2
      PROD=0.0D0
      DO IPRO=1,ND
         PROD=PROD+FJAC(IPRO,ND2)*Z1(IPRO)
      END DO
      II=ND+1
      DO I=1,LEXP  
         FACT=FJAC(3,II+1)
         IF (I.EQ.1.OR.FACT.EQ.0.0D0) THEN
            FJINV=1.0D0/(FAC1-FJAC(2,II+2))
            Z1(II)=(Z1(II)+PROD)*FJINV
         ELSE
            Z1(II)=(Z1(II)+FACT*Z1(II-1))*FJINV
         END IF
         II=II+1
      END DO
      Z1(N-1)=0.0D0
      Z1(N)=0.0D0
      RETURN
      END
C
C     END OF SUBROUTINE SOLEXP
C
C ***********************************************************
C
      SUBROUTINE SOLEXPC(N,FJAC,LDJAC,MUJAC,NM1,LDE1,E2R,E2I,Z2,Z3,
     &             IP1,IP2,ALPHN,BETAN,IJOB)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DIMENSION FJAC(LDJAC,N),E2R(LDE1,NM1),E2I(LDE1,NM1)
      DIMENSION IP1(NM1),IP2(NM1),Z2(N),Z3(N)
      ND=MUJAC+1
      LEXP=N-ND-2
      II=ND+1
      SUMR=0.0D0
      SUMI=0.0D0
      DO I=1,LEXP
         FACT=FJAC(3,II+1)
         IF (I.EQ.1.OR.FACT.EQ.0.0D0) THEN
            ALGAM=ALPHN-FJAC(2,II+2)
            DENOM=ALGAM**2+BETAN**2
            UR=(Z2(II)*ALGAM+Z3(II)*BETAN)/DENOM
            UI=(Z3(II)*ALGAM-Z2(II)*BETAN)/DENOM
         ELSE
            UR=FACT*UR
            UI=FACT*UI
            USAVE=((Z2(II)+UR)*ALGAM+(Z3(II)+UI)*BETAN)/DENOM
            UI=((Z3(II)+UI)*ALGAM-(Z2(II)+UR)*BETAN)/DENOM
            UR=USAVE
         END IF
         SUMR=SUMR+FJAC(1,II+2)*UR
         SUMI=SUMI+FJAC(1,II+2)*UI
         II=II+1
      END DO
      DO I=1,ND
         Z2(I)=Z2(I)+SUMR*FJAC(I,ND+1)
         Z3(I)=Z3(I)+SUMI*FJAC(I,ND+1)
      END DO
      CALL SOLC (ND,LDE1,E2R,E2I,Z2,Z3,IP2)
      ND2=ND+2
      PRODR=0.0D0
      PRODI=0.0D0
      DO IPRO=1,ND
         PRODR=PRODR+FJAC(IPRO,ND2)*Z2(IPRO)
         PRODI=PRODI+FJAC(IPRO,ND2)*Z3(IPRO)
      END DO
      II=ND+1
      DO I=1,LEXP
         FACT=FJAC(3,II+1)
         IF (I.EQ.1.OR.FACT.EQ.0.0D0) THEN
            ALGAM=ALPHN-FJAC(2,II+2)
            DENOM=ALGAM**2+BETAN**2
            UR=((Z2(II)+PRODR)*ALGAM+(Z3(II)+PRODI)*BETAN)/DENOM
            UI=((Z3(II)+PRODI)*ALGAM-(Z2(II)+PRODR)*BETAN)/DENOM
         ELSE
            UR=FACT*UR
            UI=FACT*UI
            USAVE=((Z2(II)+UR)*ALGAM+(Z3(II)+UI)*BETAN)/DENOM
            UI=((Z3(II)+UI)*ALGAM-(Z2(II)+UR)*BETAN)/DENOM
            UR=USAVE
         END IF
         Z2(II)=UR
         Z3(II)=UI
         II=II+1
      END DO
      Z2(N-1)=0.0D0
      Z2(N)=0.0D0
      Z3(N-1)=0.0D0
      Z3(N)=0.0D0
      RETURN
      END
C
C     END OF SUBROUTINE SOLEXPC
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,2,55,55,55,55,55,55,2,2,55,55,55), IJOB
C
   1  CONTINUE
C ------  B=IDENTITY, JACOBIAN A FULL MATRIX
      NL=N
      DO  I=1,NL
         F2(I)=HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I)
         CONT(I)=F2(I)+Y0(I)
      END DO
      CALL SOL (NL,LDE1,E1,CONT,IP1) 
      GOTO 77
C
   3  CONTINUE
C ------  B IS A DIAGONAL MATRIX, JACOBIAN A FULL MATRIX
      NL=N
      DO I=1,NL
         F2(I)=FMAS(MBDIAG,I)*(HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I))
         CONT(I)=F2(I)+Y0(I)
      END DO
      CALL SOL (NL,LDE1,E1,CONT,IP1) 
      GOTO 77
C
   2  CONTINUE
C ---  B = IDENTITY, SUMEXP LINEAR ALGEBRA (only for mi =1)
      NL=N-2
      BB=1.0D0
      DO I=1,NL
         IF (IJOB.EQ.4) BB=FMAS(MBDIAG,I)
         F2(I)=BB*(HEE1*Z1(I)+HEE2*Z2(I)+HEE3*Z3(I))
         CONT(I)=F2(I)+Y0(I)
      END DO
      CALL SOLEXP(N,FJAC,LDJAC,MUJAC,NM1,LDE1,E1,CONT,IP1,FAC1,IJOB)
      GOTO 77
C
C --------------------------------------
C
  77  CONTINUE
      ERR=0.D0
      DO  I=1,NL
         ERR=ERR+(CONT(I)/SCAL(I))**2
      END DO
      ERR=MAX(SQRT(ERR/NL),1.D-10)
C
      IF (ERR.LT.1.D0) RETURN
      IF (FIRST.OR.REJECT) THEN
          DO I=1,NL
             CONT(I)=Y(I)+CONT(I)
          END DO
          CALL FCN(N,X,CONT,F1,RPAR,IPAR)
          NFCN=NFCN+1
          DO I=1,NL
             CONT(I)=F1(I)+F2(I)
          END DO
          GOTO (31,32,31,32,31,32,55,55,55,55,32,32,31,32,31), IJOB
C ------ FULL MATRIX OPTION
  31      CONTINUE
          CALL SOL(N,LDE1,E1,CONT,IP1) 
          GOTO 88
C ---  B=IDENTITY, SUMEXP LINEAR ALGEBRA (only for mi =2)
  32      CONTINUE
          CALL SOLEXP(N,FJAC,LDJAC,MUJAC,NM1,LDE1,E1,CONT,IP1,
     &                 FAC1,IJOB)
          GOTO 88
C -----------------------------------
   88     CONTINUE
          ERR=0.D0 
          DO I=1,NL
             ERR=ERR+(CONT(I)/SCAL(I))**2
          END DO
          ERR=MAX(SQRT(ERR/NL),1.D-10)
       END IF
       RETURN
C -----------------------------------------------------------
  55   CONTINUE
       RETURN
       END
C
C     END OF SUBROUTINE ESTRAD
C
