      SUBROUTINE MULCON (FY,DFY,BC,IVPSOL,N,M,T,X,XW,TAU,TAUMIN,TAUMAX,
     &   UMAX,EPS,INFO,RWORK,LRWORK,IWORK,LIWORK)
C*    Begin Prologue MULCON
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION RWORK(LRWORK),IWORK(LIWORK),T(M),X(N,M),XW(N,M),INFO(9)
      EXTERNAL FY,DFY,BC,IVPSOL
C
C---------------------------------------------------------------------
C
C*  Title
C
C     (Mul)tiple shooting for parameter dependent two-point boundary 
C     value problems with (Con)tinuation method.
C
C*  Written by        P. Deuflhard, P. Kunkel
C*  Purpose           Solution of parameter dependent two-point boundary 
C                     value problems.
C*  Method            Numerical pathfollowing with automatic steplength 
C                     control
C*  Category          I1b2b - Differential and integral equations
C                     Parameter-dependent general BVP's
C*  Keywords          Numerical pathfollowing, Homotopy Method,
C                     Nonlinear boundary value problems,
C                     Multiple shooting
C*  Version           0.9
C*  Revision          July 1985
C*  Latest Change     January 1991
C*  Library           CodeLib
C*  Code              Fortran 77, Double Precision
C*  Environment       Standard Fortran 77 environment on PC's,
C                     workstations and hosts.
C*  Copyright     (c) Konrad-Zuse-Zentrum fuer
C                     Informationstechnik Berlin (ZIB)
C                     Takustrasse 7, D-14195 Berlin-Dahlem
C                     phone : + 49/30/84185-0
C                     fax   : + 49/30/84185-125
C*  Contact           Bodo Erdmann
C                     ZIB, Division Scientific Computing, 
C                          Department Numerical Analysis and Modelling
C                     phone : + 49/30/84185-185
C                     fax   : + 49/30/84185-107
C                     e-mail: erdmann@zib.de
C
C*    Reference:
C
C     /1/ P. Deuflhard, B. Fiedler, P. Kunkel:
C         Efficient Numerical Pathfollowing Beyond Critical Points
C         Univ. Heidelberg, Sfb 123, Tech. Rep. 278 (1984)
C
C     /2/ In: P. Deuflhard, B. Engquist (Ed.) :
C         Large Scale Scientific Computing
C         Birkhauser/Boston, Series "Progress In Scientific Computing",
C         Vol. 7 (1987), P. 97-113
C
C  ---------------------------------------------------------------
C
C* Licence
C    You may use or modify this code for your own non commercial
C    purposes for an unlimited time.
C    In any case you should not deliver this code without a special
C    permission of ZIB.
C    In case you intend to use the code commercially, we oblige you
C    to sign an according licence agreement with ZIB.
C
C* Warranty
C    This code has been tested up to a certain level. Defects and
C    weaknesses, which may be included in the code, do not establish
C    any warranties by ZIB. ZIB does not take over any liabilities
C    which may follow from aquisition or application of this code.
C
C* Software status
C    This code is not under special care of ZIB and belongs to ZIB
C    software class 3.
C
C     ------------------------------------------------------------
C
C*   Summary
C
C     Continuation method for two-point boundary value problems
C                         x(t)'=f(t,x(t),tau)
C                             0=r(x(t(1)),x(t(m)),tau)
C
C     Numerical pathfollowing with automatic steplength control
C     Optional computation of turning and (simple) bifurcation points
C
C     ------------------------------------------------------------
C
C
C  INPUT PARAMETERS  (* MARKS INOUT PARAMETERS)
C
C     - FY (T,X,TAU,F)
C               EXTERNAL SUBROUTINE DEFINING THE SYSTEM OF
C               ALGEBRAIC EQUATIONS
C                 T          TIME VARIABLE                INPUT
C                 X(N)       VECTOR OF STATE VARIABLES    INPUT
C                 TAU        HOMOTOPY PARAMETER           INPUT
C                 F(N)       VALUES OF F(T,X,TAU)         OUTPUT
C
C    - DFY (T,X,TAU,DFDX,DFDTAU)
C              EXTERNAL SUBROUTINE REALIZING THE JACOBIAN MATRIX
C              CORRESPONDING TO TO ALGEBRAIC SYSTEM DEFINED BY FY
C                 T          TIME VARIABLE                INPUT
C                 X(N)       VECTOR OF STATE VARIABLES    INPUT
C                 TAU        HOMOTOPY PARAMETER           INPUT
C                 DFDX(N,N)  JACOBIAN MATRIX DF/DX        OUTPUT
C                 DFDTAU(N)  GRADIENT DF/DTAU             OUTPUT
C
C    - BC (XA,XB,TAU,R)
C              EXTERNAL SUBROUTINE DEFINING THE BOUNDARY CONDITIONS
C                 XA(N)      STATE VARIABLES ON LEFT BOUNDARY   INPUT
C                 XB(N)      STATE VARIABLES ON RIGHT BOUNDARY  INPUT
C                 TAU        HOMOTOPY PARAMETER                 INPUT
C                 R(N)       VALUES OF R(X(T(1),X(T(M)),TAU)    OUTPUT
C
C    - IVPSOL
C               EXTERNAL SUBROUTINE SUPPLYING AN INITIAL VALUE PROBLEM
C               INTEGRATOR.
C               USE OF THE SUPPLIED INTEGRATOR  MUDIFX  IS RECOMMENDED.
C                  PARAMETERS: SEE MUDIFX.
C
C    - N        NUMBER OF ALGEBRAIC EQUATIONS
C
C    - M        NUMBER OF MULTIPLE SHOOTING POINTS
C
C    - T(M)     VECTOR OF THE MULTIPLE SHOOTING POINTS
C
C    - X(N)   * ESTIMATE OF SOLUTION ON CONTINUATION PATH
C               FOR INITIAL TAU
C
C    - XW(N)    INITIAL SCALING QUANTITIES FOR X(N)
C               (MULCON USES ADAPTIVE INTERNAL SCALING
C               ALONG THE CONTINUATION PATH)
C
C    - TAU    * INITIAL VALUE OF PARAMETER TAU
C
C    - TAUMIN   MINIMUM ALLOWED VALUE FOR TAU
C
C    - TAUMAX   MAXIMUM ALLOWED VALUE FOR TAU
C
C    - EPS      REQUIRED RELATIVE ACCURACY
C
C    - INFO(9)  ARRAY USED FOR COMMUNICATION BETWEEN
C               THE PARTICULAR ROUTINES
C                 INFO(1)  PRINT PARAMETER
C                            0  NO PRINT
C                            1  PRINT OF X(N), TAU, AND INFORMATION
C                               ILLUSTRATING THE CONTINUATION PROCESS
C                            2  ADDITIONALLY INFORMATION ABOUT THE
C                               CONVERGENCE BEHAVIOR OF THE
C                               GAUSS-NEWTON METHOD
C                 INFO(2)  MAXIMUM PERMITTED NUMBER OF
C                          CONTINUATION STEPS
C                 INFO(3)  INITIAL DIRECTION
C                           +1  PATHFOLLOWING ALONG A SINGLE BRANCH
C                               STARTING IN POSITIVE TAU-DIRECTION
C                           -1  PATHFOLLOWING ALONG A SINGLE BRANCH
C                               STARTING IN NEGATIVE TAU-DIRECTION
C                 INFO(4)  CRITICAL POINT OPTIONS
C                            0  NO CRITICAL POINT DETERMINATION
C                            1  TURNING POINT COMPUTATION
C                            2  BIFURCATION POINT COMPUTATION
C                            3  TURNING AND BIFURCATION POINT
C                               COMPUTATION  (RECOMMENDED OPTION)
C                 INFO(5)  ERROR EXIT PARAMETER  (SEE OUTPUT PARAMETERS)
C                 INFO(6)  INTERNALLY USED
C                 INFO(7)  INTERNALLY USED
C                 INFO(8)  INTERNALLY USED
C                 INFO(9)  INTERNALLY USED
C
C    - RWORK    REAL WORKSPACE
C
C    - LRWORK   LENGTH OF REAL WORKSPACE
C               MUST BE AT LEAST   3 * N**2 + 25 * N + 19
C
C    - IWORK    INTEGER WORKSPACE
C
C    - LIWORK   LENGTH OF INTEGER WORKSPACE
C               MUST BE AT LEAST   N + 1
C
C    REMARK: FOR EASE OF IMPLEMENTATION ONLY REGULAR SOLUTIONS ARE
C            ASSUMED NEAR INITIAL TAU, TAUMIN, AND TAUMAX. COMPUTATIONS
C            ARE PERFORMED STRICTLY BETWEEN TAUMIN AND TAUMAX.
C
C
C  OUTPUT PARAMETERS
C
C    - X(N)     FINAL SOLUTION VALUES
C
C    - TAU      FINAL PARAMETER VALUE
C
C    - INFO(5)  ERROR EXIT PARAMETER
C                 0  NO ERROR OCCURRED
C                 1  MORE THAN INFO(2) CONTINUATION STEPS
C                 2  STEPLENGTH IN CONTINUATION PROCESS TOO SMALL
C                    (RELATIVE DIFFERENCE LESS THAN 10*EPS)
C                 3  RANK DEFICIENT SOLUTION OBTAINED
C                 4  INITIAL GUESS OF X(N) TOO BAD
C                 5  RANK DEFICIENT JACOBIAN
C                 6  WORKSPACE TOO SMALL
C
C
C  EXTERNAL UNITS
C
C      COMMON /MUUNIT/ UPR,UDIAG
C
C      - UPR     PRINT UNIT
C                (STANDARD UNIT 6)
C      - UDIAG   PLOT INFORMATION UNIT, INPUT TO PLOT ROUTINE PLTHM
C                DISC FILE WITH CARD IMAGE  TO BE DECLARED BY THE USER
C                SET UDIAG=0, IF NO PLOT IS DESIRED
C                (STANDARD UNIT 2)
C
C*    End Prologue MULCON
      INTEGER UPR,UDIAG
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /MUUNIT/ UPR,UDIAG
C  MACHINE DEPENDENT CONSTANTS
C            EPMACH  UNIT ROUND OFF
C            SMALL   SQRT OF SMALLEST MACHINE NUMBER
      CALL ZIBCONST(EPMACH,SMALL)
C  WORKSPACE REQUIREMENTS
      N1=N+1
      NN1=N*N1
      N1Q=N1*N1
      M1=M-1
      NDIS=N*M+1
      NQ=N*N
      LRWREQ=3*N**2+25*N+19
      LIWREQ=N+1
C.    IF (LRWREQ.GT.LRWORK .OR. LIWREQ.GT.LIWORK) GOTO 9000
C  REAL WORKSPACE DISTRIBUTION
      LY=1
      LYQ=LY+NDIS
      LYD=LYQ+NDIS
      LYDH=LYD+NDIS
      LYH=LYDH+NDIS
      LDYT=LYH+NDIS
      LDYTA=LDYT+NDIS
      LDYTH=LDYTA+NDIS
      LV=LDYTH+NDIS
      LD=LV+NDIS
      LYW=LD+N1
      LG=LYW+NDIS
      LP=LG+NQ*M1
      LA=LP+N*M
      LB=LA+NQ
      LE=LB+NQ
      LEH=LE+NN1
      LES=LEH+N1Q
      LYA=LES+NN1
      LDY=LYA+NDIS
      LDYQ=LDY+NDIS
      LW=LDYQ+NDIS
      LXU=LW+NDIS
      LHH=LXU+N*M1
      LX1=LHH+N*M1
      LXM=LX1+N
      LDX1=LXM+N
      LWM=LDX1+N1
      LR=LWM+NN1
      LDE=LR+N
      LF=LDE+N
      LFH=LF+N
      LU=LFH+N
      LUS=LU+N
      LDU=LUS+N
      LDR=LDU+N
      LDHH=LDR+N
      LDDX=LDHH+N*M1
      LRF=LDDX+N*M
      LW1=LRF+M
      LW2=LW1+NDIS
C  INTEGER WORKSPACE DISTRIBUTION
      LIPIV=1
      LIROW=LIPIV+N1
      LICOL=LIROW+N
      LICOLB=LICOL+N
      LIA=LICOLB+N
      LIB=LIA+NQ
C  REORGANIZE INITIAL VALUE
      DO 1000 J=1,M
      DO 1000 I=1,N
      RWORK(I+(J-1)*N)=X(I,J)
1000  CONTINUE
      RWORK(NDIS)=TAU
C  EXTERNAL INITIAL SCALING
      DO 2000 J=1,M
      DO 2000 I=1,N
      RWORK(LYW+I+(J-1)*N-1)=XW(I,J)
2000  CONTINUE
      RWORK(LYW+NDIS-1)=DMAX1(DABS(TAU),0.01D0*DABS(TAUMIN),
     &   0.01D0*DABS(TAUMAX))
C  CALL CONTINUATION ROUTINE
      CALL MUHOMQ(FY,BC,IVPSOL,N,N1,M,M1,NDIS,TAUMIN,TAUMAX,UMAX,EPS,
     &   INFO,T,RWORK(LY),RWORK(LYQ),RWORK(LYD),RWORK(LYDH),RWORK(LYH),
     &   RWORK(LDYT),RWORK(LDYTA),RWORK(LDYTH),RWORK(LV),
     &   RWORK(LD),RWORK(LYW),RWORK(LG),RWORK(LP),RWORK(LA),RWORK(LB),
     &   RWORK(LE),RWORK(LEH),RWORK(LES),
     &   RWORK(LYA),RWORK(LDY),RWORK(LDYQ),RWORK(LW),
     &   RWORK(LXU),RWORK(LHH),RWORK(LX1),RWORK(LXM),
     &   RWORK(LDX1),RWORK(LWM),RWORK(LR),RWORK(LDE),
     &   RWORK(LF),RWORK(LFH),RWORK(LU),
     &   RWORK(LW1),RWORK(LW2),IWORK(LIPIV),
     &   IWORK(LIROW),IWORK(LICOL),IWORK(LICOLB),IWORK(LIA),IWORK(LIB),
     &   RWORK(LUS),RWORK(LDU),RWORK(LDHH),RWORK(LDR),RWORK(LDDX),
     &   RWORK(LRF),DFY)
C  RESTORE FINAL VALUES
      DO 3000 J=1,M
      DO 3000 I=1,N
      X(I,J)=RWORK(I+(J-1)*N)
3000  CONTINUE
      TAU=RWORK(NDIS)
      RETURN
C  FAIL EXIT
9000  CONTINUE
      INFO(5)=6
      KPRINT=INFO(1)
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0 .AND. LRWREQ.GT.LRWORK) WRITE(UPR,60002) LRWREQ
      IF (KPRINT.GT.0 .AND. LIWREQ.GT.LIWORK) WRITE(UPR,60003) LIWREQ
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      RETURN
60001 FORMAT(///)
60002 FORMAT('   REAL  WORKSPACE NOT SUFFICIENT',4X,'AT LEAST',I5,
     &   ' LOCATIONS NEEDED')
60003 FORMAT(' INTEGER WORKSPACE NOT SUFFICIENT',4X,'AT LEAST',I5,
     &   ' LOCATIONS NEEDED')
      END
C
C
      SUBROUTINE MUHOMQ (FY,BC,IVPSOL,N,N1,M,M1,NDIS,
     &   TAUMIN,TAUMAX,UMAX,EPS,INFO,TQ,Y,YQ,YD,YDH,YH,
     &   DYT,DYTA,DYTH,V,D,YW,G,P,A,B,E,EH,ES,YA,DY,DYQ,W,XU,HH,X1,XM,
     &   DX1,WM,RQ,DE,F,FH,U,W1,W2,IPIV,IROW,ICOL,ICOLB,IA,IB,US,DU,
     &   DHH,DR,DDX,RF,DFY)
C
C  SUBROUTINE  MUHOMQ TO BE USED WITH ROUTINE MULCON
C
C  DRIVER ROUTINE FOR MULCON
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION TQ(M),Y(NDIS),YQ(NDIS),YD(NDIS),YDH(NDIS),
     $YH(NDIS),DYT(NDIS),DYTA(NDIS),DYTH(NDIS),V(NDIS),D(N1),YW(NDIS),
     $G(N,N,M1),P(N,M),A(N,N),B(N,N),E(N,N1),EH(N1,N1),ES(N,N1),
     $YA(NDIS),DY(NDIS),DYQ(NDIS),W(NDIS),
     $XU(N,M1),HH(N,M1),X1(N),XM(N),DX1(N1),WM(N,N1),RQ(N),DE(N),
     $F(N),FH(N),U(N),W1(NDIS),W2(NDIS),
     $US(N),DU(N),DHH(N,M1),DR(N),DDX(N,M),RF(M),
     $HOF(999),HOJ(999)
      DOUBLE PRECISION MUSCPR
      INTEGER KOF(999),KOJ(999)
      INTEGER INFO(9),IPIV(N1),IROW(N),ICOL(N),ICOLB(N),IA(N,N),IB(N,N)
      INTEGER UPR,UDIAG
      EXTERNAL FY,DFY,BC,IVPSOL
      COMMON /COUNT/ IFCTEV,ITER,IDECS,ISOLS,ITRAJ
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /MUUNIT/ UPR,UDIAG
      COMMON /MUVALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      COMMON /MUMETH/ FACTOR,IEXTR,IJACM
C-----------------------------------------------------------------------
C  PREPARATIONS
      ITER=0
      IHALT=0
      IRED=0
      KPRINT=INFO(1)
      IFCTEV=0
      IDECS=0
      ISOLS=0
      ITRAJ=0
      INFO(6)=IGNMAX
      DO 1010 I=1,N
      DO 1010 K=1,N
      IA(I,K)=0
      IB(I,K)=0
      G(I,K,1)=0.D0
1010  CONTINUE
      NM=N*M
      TOL=SMALL1*EPS
      TOLF=SMALL1*EPS
      TOLJ=DSQRT(TOLF)
      IF (TOLJ.GT.1.D-3) TOLJ=1.D-3
      RELDIF=SMALL1**2
      XTHR=SMALL
      HSTART=SMALL1*(TQ(2)-TQ(1))
      NYMAX=M1
      INITF=0
      INITJ=0
C  CHECK FOR CONSISTENCY OF INITIAL VALUES
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60201)
      IF (KPRINT.GT.2) WRITE(UPR,60500) Y(NDIS),(Y(I),I=1,NM)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INFO(8)=1
      CALL MUGNHO(N,N1,M,M1,NDIS,Y,FY,BC,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,E,EH,ES,YA,DY,DYQ,W,
     &   IROW,ICOL,ICOLB,IA,IB,TQ,IVPSOL,XU,HH,X1,XM,DX1,WM,RQ,G,A,B,P,
     &   DE,F,FH,U,HSTART,XTHR,RELDIF,NE,NB,US,DU,DHH,DR,DDX,RF,NEA,
     &   NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,TOLF,TOLJ)
      IERR=INFO(9)
      IF (IERR.NE.0) GOTO 9940
      ITS=INFO(7)
1500  CONTINUE
      CALL MUTAND(N,N1,NDIS,FY,Y,F,FH,YW,E,COND,DETFA,IPIV,EH,
     &   V,D,IRANK,BC,IVPSOL,M,M1,TQ,X1,XM,RQ,XU,HH,W,HSTART,TOLJ,WM,
     &  INFO,A,B,P,RELDIF,NE,NB,IA,IB,DE,IROW,ICOL,ICOLB,XTHR,G,ES,UNRM,
     &  NEA,DFY,HOJ,KOJ,INITJ)
      IERR=INFO(9)
      IF (IERR.NE.0) GOTO 9940
      NEA=NE
      NE1=NE+1
      IPIVS=IPIV(NE1)
      IF (IPIVS.NE.NE1) IPIVSQ=ICOL(IPIVS)
      IF (IPIVS.EQ.NE1) IPIVSQ=NDIS
      DETHMA=DETFA
      IF (IPIVS.NE.NE1) CALL MUDETH(N,N1,NE,NE1,E,WM,D,NE1,IPIV,DETHMA)
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
      IF (KPRINT.GT.0) WRITE(UPR,60502) Y(NDIS),UNRM
      DO 1020 J=1,M
      IF (KPRINT.GT.0) WRITE(UPR,60503) J,TQ(J),(Y(I+(J-1)*N),I=1,N)
1020  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
C.......................................................................
C  START OF BRANCH FOLLOWING
      INIT=0
      ITER=1
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60202) ITER
      DO 1110 I=1,NDIS
      YQ(I)=Y(I)
1110  CONTINUE
C.......................................................................
C  INITIAL ESTIMATE FOR STEPLENGTH
      IRANKB=NB
      DO 1201 I=1,NE
      US(I)=0.D0
1201  DYT(IPIV(I))=V(I)
      DYT(IPIV(NE1))=-1.D0
      IF (IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUREFN(N,N1,NE,NE1,ES,DYT,DU,US,IRANKB,W,IRANK,D,IPIV,
     &   WM,EPH,NE1,EH,EPS,E,2,INFO)
      IF (INFO(9).NE.0) GOTO 9940
      DO 1202 I=1,NE1
1202  V(I)=DYT(IPIV(I))
      DO 1205 I=1,N
      DYT(I)=0.D0
1205  CONTINUE
      DO 1210 I=1,NE
      IF (IPIV(I).NE.NE1)
     &   DYT(ICOL(IPIV(I)))=-V(I)/D(NE1)*YW(ICOL(IPIV(I)))
      IF (IPIV(I).EQ.NE1) DYT(NDIS)=-V(I)/D(NE1)*YW(NDIS)
1210  CONTINUE
      IF (IPIVS.NE.NE1) DYT(ICOL(IPIVS))=1.D0/D(NE1)*YW(ICOL(IPIVS))
      IF (IPIVS.EQ.NE1) DYT(NDIS)=1.D0/D(NE1)*YW(NDIS)
      DO 40000 I=1,N
      DO 40000 J=1,M1
40000 HH(I,J)=0.D0
      CALL MURECU(N,M,M1,1,HH,G,P,DYT,DYT(NDIS),V,W)
      IF (NYMAX.GT.0 .AND. NE.GT.0)
     &   CALL MUSWP(N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,DYT,DHH,
     &   RQ,A,B,DR,DU,DE,W,WM,IROW,DX1,NE,NE1,IRANK,IRANKB,IPIV,D,E,EH,
     &   ICOL,YW,YW(NDIS),DDX,RF,DYT(NDIS),3,RELDIF,TOLF,TOLJ,INFO)
      IF (INFO(9).NE.0) GOTO 9940
      IF (NB.EQ.0) GOTO 1213
      DO 1212 I=1,NB
      DYT(N*M1+ICOLB(I))=0.D0
1212  CONTINUE
1213  CONTINUE
      S=MUSCPR(N,M,DYT,DYT,DYT(NDIS),DYT(NDIS),YW,YW(NDIS),TQ)
      S=DSQRT(S)
      DO 1211 I=1,NDIS
      DYT(I)=DYT(I)/(S*YW(I))
1211  CONTINUE
      SIGMA=SMALL1
      SIGNUM=DBLE(INFO(3))
      IF (DYT(NDIS)*SIGNUM.GT.0.D0) GOTO 1221
      DO 1220 I=1,NDIS
      DYT(I)=-DYT(I)
1220  CONTINUE
1221  CONTINUE
      IPIVA=IPIVS
      IPIVAQ=IPIVSQ
      IF (UDIAG.GT.0) WRITE(UDIAG,20102) IPIVSQ,UNRM
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) Y,DYT,YW
      GOTO 2200
C-----------------------------------------------------------------------
C  STEPLENGTH PREDICTOR
2000  CONTINUE
      ITER=ITER+1
      IF (ITER.GT.INFO(2)) GOTO 9910
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60202) ITER
      IF (ITS.EQ.0) GOTO 2020
      TH=DSQRT(SUMYQH/SUMYH)
      VAL1=0.25D0*THMAX*DSQRT(SUMYH)/DABS(VAL1)
      DO 2010 I=1,NDIS
      W(I)=YDH(I)-YQ(I)
2010  CONTINUE
      VAL2=TH*DSQRT(MUSCPR(N,M,W,W,W(NDIS),W(NDIS),YW,YW(NDIS),TQ))
      IF (VAL2.LT.THR*VAL1) VAL2=THR*VAL1
      R=DSQRT(VAL1/VAL2)
      SIGMA=R*SIGMA
      GOTO 2030
C  EMPIRICAL STEPLENGTH INCREASE IN NEARLY LINEAR CASE
2020  CONTINUE
      SIGMA=SIGMA/DSQRT(THR)
C  EXTRAPOLATED STEPLENGTH BOUND
2030  CONTINUE
      IF (IHALT.EQ.0) GOTO 2040
      SIGMAQ=FACTOR*(YHALT-YQ(IPIVAQ))/(DYT(IPIVAQ)*YW(IPIVAQ))
      IF (SIGMAQ.LE.10.D0*EPS .OR. SIGMA.LE.SIGMAQ) GOTO 2040
      SIGMA=SIGMAQ
      IF (KPRINT.GT.1) WRITE(UPR,60206)
2040  CONTINUE
      IF (INIT.NE.0) IPIVA=IPIVS
      IF (INIT.NE.0) IPIVAQ=IPIVSQ
      GOTO 2200
C.......................................................................
C  STEPLENGTH CORRECTOR
2110  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60401)
      TH=DSQRT(SUMYQH/SUMYH)
      R=DSQRT(0.25D0*THMAX/TH)
      IF (R.GT.0.7D0) R=0.7D0
      IF (R.LT.0.1D0) R=0.1D0
      GOTO 2190
2120  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60402) ICON
      R=0.7D0
      GOTO 2190
2130  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60403)
      R=0.7D0
2190  CONTINUE
      IF (SIGMA.LT.EPS) GOTO 9920
      SIGMA=R*SIGMA
      INIT=1
C.......................................................................
C  TRIAL VALUE
2200  CONTINUE
      IF (DYT(NDIS).EQ.0.D0) DYT(NDIS)=SMALL
      SIGMAX=(TAUMAX-YQ(NDIS))/(DYT(NDIS)*YW(NDIS))
      IF (SIGMAX.GT.0.D0 .AND. SIGMA.GT.SIGMAX) SIGMA=SIGMAX
      SIGMIN=(TAUMIN-YQ(NDIS))/(DYT(NDIS)*YW(NDIS))
      IF (SIGMIN.GT.0.D0 .AND. SIGMA.GT.SIGMIN) SIGMA=SIGMIN
      DO 2220 I=1,NDIS
      YD(I)=YQ(I)+SIGMA*DYT(I)*YW(I)
      YDH(I)=YD(I)
2220  CONTINUE
C.......................................................................
C  RESCALING
      IF (IERR.NE.0) GOTO 2490
      CALL MUSCLE(N,M,M1,YQ,DYT,DYT(NDIS),DYTH,DYTH(NDIS),
     &   YW,YW(NDIS),XTHR,YQ(NDIS),TAUMIN,TAUMAX,TQ)
      DO 2420 I=1,NDIS
      Y(I)=YQ(I)
2420  CONTINUE
2490  CONTINUE
      IF (SIGMA.EQ.SIGMIN .OR. SIGMA.EQ.SIGMAX) GOTO 3100
C-----------------------------------------------------------------------
C  ITERATION BACK TO CONTINUATION PATH BY GAUSS-NEWTON METHOD
      IF (KPRINT.GT.0) WRITE(UPR,60205) SIGMA,IPIVA
      IF (KPRINT.GT.2) WRITE(UPR,60500) YD(NDIS),(YD(I),I=1,NM)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INFO(8)=2
      CALL MUGNHO(N,N1,M,M1,NDIS,YD,FY,BC,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,E,EH,ES,YA,DY,DYQ,W,
     &   IROW,ICOL,ICOLB,IA,IB,TQ,IVPSOL,XU,HH,X1,XM,DX1,WM,RQ,G,A,B,P,
     &   DE,F,FH,U,HSTART,XTHR,RELDIF,NE,NB,US,DU,DHH,DR,DDX,RF,NEA,
     &   NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,TOLF,TOLJ)
      IERR=INFO(9)
      ITS=INFO(7)
      IF (IERR.NE.0 .OR. ITS.NE.0) INIT=1
      IF (INIT.EQ.0) GOTO 2020
      GOTO (2110,2110,9950,3110,2110,9970,5300),IERR
      DO 3010 I=1,NDIS
      W(I)=DYT(I)*YW(I)
3010  CONTINUE
      VAL1=MUSCPR(N,M,W,V,W(NDIS),V(NDIS),YW,YW(NDIS),TQ)
      CALL MUTAND(N,N1,NDIS,FY,YD,F,FH,YW,E,COND,DETF,IPIV,EH,
     &   V,D,IRANK,BC,IVPSOL,M,M1,TQ,X1,XM,RQ,XU,HH,W,HSTART,TOLJ,WM,
     &  INFO,A,B,P,RELDIF,NE,NB,IA,IB,DE,IROW,ICOL,ICOLB,XTHR,G,ES,UNRM,
     &  NEA,DFY,HOJ,KOJ,INITJ)
      IERR=INFO(9)
      GOTO (9950,2110,9970,5300),IERR
      IPIVS=IPIV(NE1)
      IF (IPIVS.NE.NE1) IPIVSQ=ICOL(IPIVS)
      IF (IPIVS.EQ.NE1) IPIVSQ=NDIS
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
      GOTO 4000
C.......................................................................
C  HIT FINAL VALUE
3100  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60205) SIGMA,IPIVA
      IF (SIGMA.EQ.SIGMIN) YD(NDIS)=TAUMIN
      IF (SIGMA.EQ.SIGMAX) YD(NDIS)=TAUMAX
      GOTO 3120
3110  CONTINUE
      IF (YD(NDIS).LT.TAUMIN) YD(NDIS)=TAUMIN
      IF (YD(NDIS).GT.TAUMAX) YD(NDIS)=TAUMAX
3120  CONTINUE
      IF (KPRINT.GT.2) WRITE(UPR,60500) YD(NDIS),(YD(I),I=1,NM)
      IF (KPRINT.GT.2) WRITE(UPR,60001)
      INIT=2
      INFO(8)=3
      CALL MUGNHO(N,N1,M,M1,NDIS,YD,FY,BC,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,E,EH,ES,YA,DY,DYQ,W,
     &   IROW,ICOL,ICOLB,IA,IB,TQ,IVPSOL,XU,HH,X1,XM,DX1,WM,RQ,G,A,B,P,
     &   DE,F,FH,U,HSTART,XTHR,RELDIF,NE,NB,US,DU,DHH,DR,DDX,RF,NEA,
     &   NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,TOLF,TOLJ)
      IERR=INFO(9)
      ITS=INFO(7)
      GOTO (2110,2110,9930,2110,2110,9970,5300),IERR
      CALL MUTAND(N,N1,NDIS,FY,YD,F,FH,YW,E,COND,DETF,IPIV,EH,
     &   V,D,IRANK,BC,IVPSOL,M,M1,TQ,X1,XM,RQ,XU,HH,W,HSTART,TOLJ,WM,
     &  INFO,A,B,P,RELDIF,NE,NB,IA,IB,DE,IROW,ICOL,ICOLB,XTHR,G,ES,UNRM,
     &  NEA,DFY,HOJ,KOJ,INITJ)
      IERR=INFO(9)
      GOTO (9930,2110,9970,5300),IERR
      IPIVS=IPIV(NE1)
      IF (IPIVS.NE.NE1) IPIVSQ=ICOL(IPIVS)
      IF (IPIVS.EQ.NE1) IPIVSQ=NDIS
      IF (KPRINT.GT.0) WRITE(UPR,60300) ITS
C-----------------------------------------------------------------------
C  DETERMINANTS
4000  CONTINUE
      DETHM=DETF
      IF (IPIVS.NE.NE1) CALL MUDETH(N,N1,NE,NE1,E,WM,D,NE1,IPIV,DETHM)
      DETFH=DETF
      IF (IPIVS.NE.IPIVA)
     &   CALL MUDETH(N,N1,NE,NE1,E,WM,D,IPIVA,IPIV,DETFH)
      IF (KPRINT.GT.1) WRITE(UPR,60501) DETFA,DETFH,DETHMA,DETHM
C.......................................................................
C  NEW NORMALIZED TANGENT
      IRANKB=NB
      DO 1203 I=1,NE
      US(I)=0.D0
1203  W(IPIV(I))=V(I)
      W(IPIV(NE1))=-1.D0
      IF (IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUREFN(N,N1,NE,NE1,ES,W,DU,US,IRANKB,V,IRANK,D,IPIV,
     &   WM,EPH,NE1,EH,EPS,E,2,INFO)
      IF (INFO(9).NE.0) GOTO 2130
      DO 1204 I=1,NE1
1204  V(I)=W(IPIV(I))
      DO 4105 I=1,N
      W(I)=0.D0
4105  CONTINUE
      DO 4110 I=1,NE
      IF (IPIV(I).NE.NE1)
     &   W(ICOL(IPIV(I)))=-V(I)/D(NE1)*YW(ICOL(IPIV(I)))
      IF (IPIV(I).EQ.NE1) W(NDIS)=-V(I)/D(NE1)*YW(NDIS)
4110  CONTINUE
      IF (IPIVS.NE.NE1) W(ICOL(IPIVS))=1.D0/D(NE1)*YW(ICOL(IPIVS))
      IF (IPIVS.EQ.NE1) W(NDIS)=1.D0/D(NE1)*YW(NDIS)
      DO 40001 I=1,N
      DO 40001 J=1,M1
40001 HH(I,J)=0.D0
      CALL MURECU(N,M,M1,1,HH,G,P,W,W(NDIS),V,WM)
      IF (NYMAX.GT.0 .AND. NE.GT.0)
     &   CALL MUSWP(N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,W,DHH,
     &   RQ,A,B,DR,DU,DE,V,WM,IROW,DX1,NE,NE1,IRANK,IRANKB,IPIV,D,E,EH,
     &   ICOL,YW,YW(NDIS),DDX,RF,W(NDIS),3,RELDIF,TOLF,TOLJ,INFO)
      IF (INFO(9).NE.0) GOTO 2130
      IF (NB.EQ.0) GOTO 4113
      DO 4112 I=1,NB
      W(N*M1+ICOLB(I))=0.D0
4112  CONTINUE
4113  CONTINUE
      VN1=MUSCPR(N,M,DYQ,DYQ,DYQ(NDIS),DYQ(NDIS),YW,YW(NDIS),TQ)
      VN2=MUSCPR(N,M,W,W,W(NDIS),W(NDIS),YW,YW(NDIS),TQ)
      WINKEL=MUSCPR(N,M,DYQ,W,DYQ(NDIS),W(NDIS),YW,YW(NDIS),TQ)
      WINKEL=WINKEL/DSQRT(VN1*VN2)
      VN2=DSQRT(VN2)
      DO 4111 I=1,NDIS
      W(I)=W(I)/(VN2*YW(I))
4111  CONTINUE
C  CHECK SIGN OF TANGENT
      SIGNUM=DSIGN(1.D0,YD(IPIVSQ)-Y(IPIVSQ))
      IF (W(IPIVSQ)*SIGNUM.GT.0.D0) GOTO 4121
      DO 4120 I=1,NDIS
      W(I)=-W(I)
4120  CONTINUE
4121  CONTINUE
C.......................................................................
C  CHECK FOR CONFLICT IN SIGNS
      IF (SIGMA.LE.10.D0*EPS) GOTO 4210
      ICON=1
      IF (DYT(NDIS)*W(NDIS).LE.0.D0 .AND.
     &   (DETHMA*DETHM.GT.0.D0 .OR. IPIVSQ+IPIVAQ.EQ.2*NDIS)) GOTO 2120
      ICON=2
      IF (DYT(NDIS)*W(NDIS).GT.0.D0 .AND. DETHMA*DETHM.LT.0.D0 .AND.
     &   DETFH*DETFA.GT.0.D0) GOTO 2120
      ICON=3
      IF (DETFH*DETFA.LE.0.D0 .AND. DETHM*DETHMA.GT.0.D0) GOTO 2120
      ICON=4
      IF (DYT(NDIS)*W(NDIS).LE.0.D0 .AND.
     &   DETHMA*DETHM.LE.0.D0 .AND. DETFH*DETFA.LE.0.D0) GOTO 2120
      ICON=5
      IF (DYT(IPIVSQ)*W(IPIVSQ).LT.0.D0) GOTO 2120
      ICON=6
      IF (DYT(IPIVAQ)*W(IPIVAQ).LT.0.D0) GOTO 2120
4210  CONTINUE
C.......................................................................
C  SAVE VALUES FOR NEXT ITERATE
      DO 4310 I=1,NDIS
      DYTA(I)=DYT(I)
      DYT(I)=W(I)
      YQ(I)=YD(I)
4310  CONTINUE
C-----------------------------------------------------------------------
C  CHECK FOR TURNING POINT
      IF (DYTA(NDIS)*DYT(NDIS).GT.0.D0) GOTO 5010
      CALL MUTURN(N,N1,IPIVA,EPS,INFO,Y,W1,YD,YH,DYT,DYTH,W2,V,
     &   YW,IPIV,FY,TAUMIN,TAUMAX,D,E,EH,YA,DY,DYQ,W,F,FH,U,
     &  M,M1,NDIS,BC,ES,IROW,ICOL,ICOLB,IA,IB,TQ,IVPSOL,XU,HH,X1,XM,DX1,
     &   WM,RQ,G,A,B,P,DE,HSTART,XTHR,RELDIF,NE,NE1,NB,IPIVAQ,UNRMQ,
     &   US,DU,DHH,DR,DDX,RF,NEA,NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,
     &   TOLF,TOLJ)
      IERRQ=INFO(9)
      IF (IERRQ.EQ.4) GOTO 9970
      IF (IERRQ.EQ.5) GOTO 5300
      IF (IERRQ.GT.0 .AND. SIGMA.GE.10.D0*EPS) GOTO 5200
      IF (IPIV(NE1).NE.NE1) IPIVH=ICOL(IPIV(NE1))
      IF (IPIV(NE1).EQ.NE1) IPIVH=NDIS
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20103) IPIVH,UNRMQ
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20000) YH,W,YW
5010  CONTINUE
C.......................................................................
C  CHECK FOR BIFURCATION POINT
      IF (DETFH*DETFA.GT.0.D0) GOTO 6000
      DO 5111 I=1,N1
      YD(I)=YQ(I)
5111  CONTINUE
      CALL MUBIFC(N,N1,IPIVA,EPS,INFO,Y,W1,YD,YH,V,
     &  YW,IPIV,FY,TAUMIN,TAUMAX,D,E,EH,YA,DY,DYQ,W,F,FH,U,DETFH,DETFA,
     &  M,M1,NDIS,BC,ES,IROW,ICOL,ICOLB,IA,IB,TQ,IVPSOL,XU,HH,X1,XM,DX1,
     &  WM,RQ,G,A,B,P,DE,HSTART,XTHR,RELDIF,NE,NE1,NB,IPIVAQ,UNRMQ,
     &  US,DU,DHH,DR,DDX,RF,NEA,NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,
     &  TOLF,TOLJ)
      IERRQ=INFO(9)
      IF (IERRQ.EQ.4) GOTO 9970
      IF (IERRQ.EQ.5) GOTO 5300
      IF (IERRQ.GT.0) GOTO 5200
      IF (IPIV(NE1).NE.NE1) IPIVH=ICOL(IPIV(NE1))
      IF (IPIV(NE1).EQ.NE1) IPIVH=NDIS
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20101) IPIVH,UNRMQ
      IF (UDIAG.GT.0 .AND. IERRQ.EQ.0) WRITE(UDIAG,20000) YH,W,YW
      GOTO 6000
C.......................................................................
C  RESTORE FORMER VALUES
5200  CONTINUE
      DO 5210 I=1,NDIS
      YQ(I)=Y(I)
      DYT(I)=DYTA(I)
5210  CONTINUE
      INIT=1
      GOTO 2130
C.......................................................................
C  RESTART DUE TO CHANGE IN PROJECTIONS
5300  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60601)
      INIT=0
      GOTO 1500
C-----------------------------------------------------------------------
C  FINISH OUTPUT FOR THIS STEP
6000  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60502) YQ(NDIS),UNRM
      DO 6001 J=1,M
      IF (KPRINT.GT.0) WRITE(UPR,60503) J,TQ(J),(YQ(I+(J-1)*N),I=1,N)
6001  CONTINUE
      IF (UDIAG.GT.0) WRITE(UDIAG,20100) IPIVSQ,UNRM
      IF (UDIAG.GT.0) WRITE(UDIAG,20000) YQ,DYT,YW
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INIT.EQ.2 .OR. UNRM.GE.UMAX) GOTO 9900
C  EXTRAPOLATED STEPLENGTH BOUND
      IHALT=0
      IF (DETFA*DETFH.LE.0.D0 .OR. DABS(DETFA).LE.DABS(DETFH) .OR.
     &   IEXTR.NE.1) GOTO 6010
      IHALT=1
      S=DETFH/DETFA
      YHALT=(YQ(IPIVAQ)-S*Y(IPIVAQ))/(1.D0-S)
6010  CONTINUE
      DETFA=DETF
      DETHMA=DETHM
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9900  CONTINUE
      INFO(5)=0
      IF (KPRINT.GT.0) WRITE(UPR,60900) IFCTEV,ITRAJ
      DO 9901 I=1,NDIS
      Y(I)=YQ(I)
9901  CONTINUE
      RETURN
C  FAIL EXIT
C  MORE THAN ITMAX STEPS
9910  CONTINUE
      INFO(5)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      RETURN
C  STEPLENGTH TOO SMALL
9920  CONTINUE
      INFO(5)=2
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      RETURN
C  RANK DEFICIENT SOLUTION OBTAINED
9930  CONTINUE
      INFO(5)=3
      IF (KPRINT.GT.0) WRITE(UPR,60502) YD(NDIS),UNRM
      DO 9931 J=1,M
      IF (KPRINT.GT.0) WRITE(UPR,60503) J,TQ(J),(YD(I+(J-1)*N),I=1,N)
9931  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60903)
      RETURN
C  STARTING VALUES TOO BAD
9940  CONTINUE
      INFO(5)=4
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60904)
      RETURN
C  RANK DEFICIENCY OF JACOBIAN
9950  CONTINUE
      INFO(5)=5
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60905)
      RETURN
C  NON-ZERO RESIDUAL IN SEPARATED BOUNDARY CONDITION
9970  CONTINUE
      INFO(5)=7
      IF (KPRINT.GT.0) WRITE(UPR,60100)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (KPRINT.GT.0) WRITE(UPR,60907)
      RETURN
C-----------------------------------------------------------------------
20000 FORMAT(4D18.10)
20100 FORMAT('0',I4,D18.10)
20101 FORMAT('1',I4,D18.10)
20102 FORMAT('2',I4,D18.10)
20103 FORMAT('3',I4,D18.10)
60001 FORMAT(/)
60002 FORMAT(///)
60100 FORMAT(' ',132('*'))
60201 FORMAT(' STEP  0       CHECK FOR CONSISTENCY OF INITIAL VALUES'/)
60202 FORMAT(' STEP',I4/)
60205 FORMAT(' SCALED STEPLENGTH',D11.4/
     &   ' CURRENT CONTINUATION PARAMETER',I4/)
60206 FORMAT(' STEPLENGTH REDUCED BY EXTRAPOLATED STEPLENGTH BOUND')
60300 FORMAT(/' GAUSS-NEWTON METHOD REQUIRED',I3,' ITERATIONS'//)
60401 FORMAT(/' GAUSS-NEWTON METHOD FAILED'/' ',132('.')/)
60402 FORMAT(/' CONFLICT NO.',I2,' IN SIGNS'/' ',132('.')/)
60403 FORMAT(/' STEP RETRIED WITH REDUCED STEPLENGTH'/' ',132('.')/)
60500 FORMAT(' TAU=',D14.6,'      X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60501 FORMAT(/' DET WRT ACTUAL CON PARM',2D17.7/
     &   ' DET WRT  REAL  CON PARM',2D17.7/)
60502 FORMAT(' TAU=',D14.6,'      NORM=',D14.6/)
60503 FORMAT('      NODE',I3,'     T=',D14.6/
     &   25X,'X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60601 FORMAT(' RESTART DUE TO CHANGE IN PROJECTIONS'/' ',132('.')/)
60900 FORMAT(' MULCON REQUIRED',I7,' FUNCTION EVALUATIONS AND',I7,
     &   ' TRAJECTORY EVALUATIONS'//)
60901 FORMAT(' TERMINATION AFTER INFO(2) CONTINUATION STEPS'//)
60902 FORMAT(' TERMINATION SINCE STEPLENGTH TOO SMALL'//)
60903 FORMAT(' RANK DEFICIENT SOLUTION OBTAINED'//)
60904 FORMAT(' TERMINATION SINCE INITIAL GUESS TOO FAR AWAY FROM',
     &   ' A REGULAR POINT OF THE CONTINUATION PATH'/
     &   ' USE NONLINEAR EQUATION SOLVER FOR BETTER INITIAL DATA'//)
60905 FORMAT(' TERMINATION SINCE JACOBIAN RANK DEFICIENT'//)
60907 FORMAT(' TERMINATION SINCE NON-ZERO RESIDUAL IN SEPARATED',
     &   ' BOUNDARY CONDITION WAS FOUND'//)
      END
C
C
      SUBROUTINE MUGNHO (N,N1,M,M1,NDIS,Y,FY,BC,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,PIVOT,V,D,YW,E,EH,ES,YA,DY,DYQ,W,
     &   IROW,ICOL,ICOLB,IA,IB,T,IVPSOL,XU,HH,X1,XM,DX1,WM,
     &   R,G,A,B,P,DE,F,FH,U,HSTART,XTHR,RELDIF,NE,NB,US,DU,
     &   DHH,DR,DDX,RF,NEA,NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,TOLF,
     &   TOLJ)
C
C  SUBROUTINE  MUGNHO TO BE USED WITH ROUTINE MULCON
C
C  GAUSS-NEWTON METHOD AS CORRECTOR IN PATHFOLLOWING PROCEDURE
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(NDIS),V(NDIS),D(N1),YW(NDIS)
      DOUBLE PRECISION E(N,N1),EH(N1,N1),ES(N,N1)
      DOUBLE PRECISION YA(NDIS),DY(NDIS),DYQ(NDIS),W(NDIS)
      DOUBLE PRECISION T(M),XU(N,M1),HH(N,M1),X1(N),XM(N),DX1(N1),WM(1)
      DOUBLE PRECISION R(N),G(N,N,M1),A(N,N),B(N,N),P(N,M),DE(N)
      DOUBLE PRECISION F(N),FH(N),U(N)
      DOUBLE PRECISION US(N),DU(N),DHH(N,M1),DR(N),DDX(N,M),RF(M)
      DOUBLE PRECISION HOF(M1),HOJ(M1),VSAVE(999),VH(999)
      DOUBLE PRECISION MUSCPR
      INTEGER PIVOT(N1),IROW(N),ICOL(N),ICOLB(N),IA(N,N),IB(N,N),INFO(9)
      INTEGER KOF(M1),KOJ(M1)
      INTEGER UPR,UDIAG
      EXTERNAL FY,DFY,BC,IVPSOL
      COMMON /COUNT/ IFCTEV,IDUMQ,IDECS,ISOLS,ITRAJ
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /MUUNIT/ UPR,UDIAG
      COMMON /MUVALS/ SMALL1,THMAX,THR,IDUM,IGNMAX
      COMMON /MUMETH/ FACTOR,IEXTR,IJACM
C  INTERNAL PARAMETERS
      SIGMA=0.5D0/THMAX
C  INITIAL PREPARATIONS TO START MUGNHO
      THMAX2=1.D0
      IF (INFO(8).NE.1) THMAX2=THMAX**2
      KPRINT=INFO(1)
      ITMAX=INFO(6)
      ITER=0
      LEVEL=0
      NEW=0
      LYM=N*M1+1
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y
C-----------------------------------------------------------------------
C  COMPUTATION OF RESIDUAL VECTOR
1     CALL MURESI(FY,BC,IVPSOL,N,N1,M,M1,T,Y,Y(NDIS),X1,XM,R,
     &   XU,HH,W,HSTART,TOLF,INFO,UNRM,YW,YW(NDIS),TFAIL,G,P,0,W,DFY,
     &   HOF,KOF,INITF)
      IF (INFO(9).NE.0) GOTO 95
      IF (LEVEL.EQ.1) GOTO 43
C-----------------------------------------------------------------------
C  DIFFERENCE APPROXIMATION OF THE SCALED JACOBIAN MATRIX A
C  INCLUDING FEED-BACK DEVICE
2     NEW=0
      CALL MUDERA(BC,N,M,Y,Y(LYM),Y(NDIS),R,W,A,B,P(1,M),
     &   YW,YW(NDIS),RELDIF)
      CALL MUCHAB(N,N1,M,NE,NB,Y,A,B,IA,IB,P(1,M),R,DE,YW,YW(NDIS),
     &   IROW,ICOL,ICOLB,XTHR,INFO)
      IF (INFO(9).NE.0) GOTO 96
      IF (INFO(8).EQ.1) NEA=NE
      IF (NE.NE.NEA) GOTO 97
      NE1=NE+1
      NQ=NE
      IF (INFO(8).EQ.2) NQ=NE1
      CALL MUDERG(FY,N,N1,M,M1,T,Y,Y(NDIS),XU,YW,YW(NDIS),W,TFAIL,G,P,
     &   IVPSOL,HSTART,TOLJ,INFO,WM,DFY,HOJ,KOJ,INITJ)
      GOTO 4
C-----------------------------------------------------------------------
C  RANK-1 APPROXIMATION OF THE JACOBIAN DUE TO BROYDEN
C  (SCALED VERSION)
3     NEW=NEW+1
301   CALL MUDERA(BC,N,M,Y,Y(LYM),Y(NDIS),R,W,A,B,P(1,M),
     &   YW,YW(NDIS),RELDIF)
      CALL MUCHAB(N,N1,M,NE,NB,Y,A,B,IA,IB,P(1,M),R,DE,YW,YW(NDIS),
     &   IROW,ICOL,ICOLB,XTHR,INFO)
      IF (INFO(9).NE.0) GOTO 96
      IF (INFO(8).EQ.1) NEA=NE
      IF (NE.NE.NEA .AND. INFO(8).NE.1) GOTO 97
      NE1=NE+1
      NQ=NE
      IF (INFO(8).EQ.2) NQ=NE1
      CALL MURK1G(N,N1,M,M1,YW,YW(NDIS),DY,DY(NDIS),HH,W,G,P)
C-----------------------------------------------------------------------
C  SOLUTION OF THE LINEAR SYSTEM
4     CALL MUCDNS(N,N1,M,M1,NE,NE1,G,A,B,YW,DE,E,ES,W,IROW,ICOL,
     &   P,YW(NDIS),WM,DYQ)
C  SAVE MATRIX E ON ES FOR REFINEMENT
401   DO 7410 K=1,NE1
      DO 7410 I=1,NE
7410  ES(I,K)=E(I,K)
C  HOUSEHOLDER TRIANGULARIZATION
      COND=1.D0/EPMACH
      IRANK=NE
      CALL MUDECC(E,N,N1,NB,NE,NQ,IRANK,COND,D,PIVOT,0,EH,V)
      IDECS=IDECS+1
      D1=DABS(D(1))
C  LINEAR LEAST SQUARES SOLUTION
43    CALL MURHS1(N,NE,M1,1,HH,R,B,G,F,DE,W,WM,IROW)
      DO 431 K=1,NE
      US(K)=F(K)
431   U(K)=F(K)
      IRANKB=NB
      CALL MUSOLC(E,N,N1,IRANKB,NE,NQ,DX1,U,IRANK,D,PIVOT,0,EH,W)
      ISOLS=ISOLS+1
      IF (IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUREFN(N,N1,NE,NE1,ES,DX1,DU,US,IRANKB,W,IRANK,D,PIVOT,
     &   WM,EPH,NQ,EH,EPS,E,LEVEL,INFO)
      IF (INFO(9).NE.0) GOTO 95
C  DESCALING
      DO 4314 I=1,N
4314  DYQ(I)=0.D0
      DO 4311 I=1,NE
      IQ=ICOL(I)
4311  DYQ(IQ)=DX1(I)*YW(IQ)
      DYQ(NDIS)=0.D0
      IF (INFO(8).EQ.2) DYQ(NDIS)=DX1(NE1)*YW(NDIS)
      CALL MURECU(N,M,M1,1,HH,G,P,DYQ,DYQ(NDIS),WM,W)
      IF (NYMAX.GT.0 .AND. IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUSWP(N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,DYQ,DHH,
     &   R,A,B,DR,DU,DE,W,WM,IROW,DX1,NE,NQ,IRANK,IRANKB,PIVOT,D,E,EH,
     &   ICOL,YW,YW(NDIS),DDX,RF,DYQ(NDIS),LEVEL,RELDIF,TOLF,TOLJ,INFO)
      IF (INFO(9).NE.0) GOTO 95
      IF (NB.EQ.0) GOTO 4313
      DO 4312 I=1,NB
      DYQ(N*M1+ICOLB(I))=0.D0
4312  CONTINUE
4313  CONTINUE
      IF (INFO(8).NE.2) GOTO 99899
      IF (LEVEL.EQ.1) GOTO 98999
      IPIVS=PIVOT(NE1)
      DO 1203 I=1,NE
1203  W(PIVOT(I))=V(I)
      W(PIVOT(NE1))=-1.D0
      DO 1208 I=1,NE
1208  US(I)=0.D0
      IF (IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUREFN(N,N1,NE,NE1,ES,W,DU,US,IRANKB,V,IRANK,D,PIVOT,
     &   WM,EPH,NE1,EH,EPS,E,2,INFO)
      IF (INFO(9).NE.0) GOTO 95
      DO 1204 I=1,NE1
1204  V(I)=W(PIVOT(I))
      DO 4105 I=1,N
      W(I)=0.D0
4105  CONTINUE
      DO 4110 I=1,NE
      IF (PIVOT(I).NE.NE1)
     &   W(ICOL(PIVOT(I)))=-V(I)/D(NE1)*YW(ICOL(PIVOT(I)))
      IF (PIVOT(I).EQ.NE1) W(NDIS)=-V(I)/D(NE1)*YW(NDIS)
4110  CONTINUE
      IF (IPIVS.NE.NE1) W(ICOL(IPIVS))=1.D0/D(NE1)*YW(ICOL(IPIVS))
      IF (IPIVS.EQ.NE1) W(NDIS)=1.D0/D(NE1)*YW(NDIS)
      DO 40001 I=1,N
      R(I)=0.D0
      DO 40001 J=1,M1
40001 HH(I,J)=0.D0
      CALL MURECU(N,M,M1,1,HH,G,P,W,W(NDIS),V,WM)
      IF (NYMAX.GT.0 .AND. NE.GT.0)
     &   CALL MUSWP(N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,W,DHH,
     &   R,A,B,DR,DU,DE,V,WM,IROW,DX1,NE,NE1,IRANK,IRANKB,PIVOT,D,E,EH,
     &   ICOL,YW,YW(NDIS),DDX,RF,W(NDIS),2,RELDIF,TOLF,TOLJ,INFO)
      IF (INFO(9).NE.0) GOTO 95
      IF (NB.EQ.0) GOTO 4113
      DO 4112 I=1,NB
      W(N*M1+ICOLB(I))=0.D0
4112  CONTINUE
4113  CONTINUE
      S=MUSCPR(N,M,W,W,W(NDIS),W(NDIS),YW,YW(NDIS),T)
      S=DSQRT(S)
      DO 4114 I=1,NDIS
4114  VH(I)=W(I)/S
      IF (ITER.NE.0 .OR. LEVEL.EQ.1 .OR. INFO(8).NE.2) GOTO 41162
      DO 41152 I=1,NDIS
41152 VSAVE(I)=VH(I)
41162 CONTINUE
98999 CONTINUE
      SK1=MUSCPR(N,M,DYQ,VH,DYQ(NDIS),VH(NDIS),YW,YW(NDIS),T)
      SK2=MUSCPR(N,M,VH,VH,VH(NDIS),VH(NDIS),YW,YW(NDIS),T)
      SK1=SK1/SK2
      DO 4116 I=1,NDIS
      DYQ(I)=DYQ(I)-SK1*VH(I)
4116  CONTINUE
99899 CONTINUE
      IF (LEVEL.EQ.1) GOTO 44
      CONV=0.D0
      DO 433 L=1,NDIS
      S=DYQ(L)
      YA(L)=Y(L)
      Y(L)=Y(L)+S
433   CONV=DMAX1(CONV,DABS(S)/YW(L))
C  TEST OF ACCURACY
      IF (CONV.LE.EPS) GOTO 9
434   SUMY=MUSCPR(N,M,DYQ,DYQ,DYQ(NDIS),DYQ(NDIS),YW,YW(NDIS),T)
      DO 435 L=1,NDIS
      S=DYQ(L)
435   DY(L)=S
      LEVEL=1
      GOTO 1
44    SUMYQ=MUSCPR(N,M,DYQ,DYQ,DYQ(NDIS),DYQ(NDIS),YW,YW(NDIS),T)
441   CONTINUE
C-----------------------------------------------------------------------
C  RESTRICTED NATURAL MONOTONICITY TEST
      SUMYQ=DMAX1(SMALL*SUMY,SUMYQ)
      IF (ITER.NE.0) GOTO 51
      SUMYH=SUMY
      SUMYQH=SUMYQ
51    CONTINUE
      IF (SUMYQ.GT.THMAX2*SUMY) GOTO 69
C  A-POSTERIORI ESTIMATE OF RELAXATION FACTOR
      FCH=0.5D0*DSQRT(SUMY/SUMYQ)
C-----------------------------------------------------------------------
C  PREPARATIONS TO START THE FOLLOWING ITERATION STEP
      LEVEL=0
      ITER=ITER+1
      IF (KPRINT.GT.1 .AND. ITER.EQ.1) WRITE(UPR,60003)
      IF (KPRINT.GT.1) WRITE(UPR,60002)
     &   ITER,SUMY,SUMYQ,NEW,IRANK,COND,D1,(PIVOT(I),I=1,NQ)
      IF (KPRINT.GT.2) WRITE(UPR,60001) Y
      IF (Y(NDIS).LT.TAUMIN .OR. Y(NDIS).GT.TAUMAX) GOTO 94
      IF (ITER.GE.ITMAX) GOTO 91
      IF (FCH.LT.SIGMA .OR. IRANK.LT.NE .OR. IJACM.EQ.2) GOTO 2
      IF (IJACM.EQ.1) GOTO 3
      NEW=NEW-1
      GOTO 43
69    CONTINUE
      IF (INFO(8).NE.1 .OR. NEW.EQ.0) GOTO 92
      DO 691 I=1,NQ
691   Y(I)=YA(I)
      LEVEL=0
      GOTO 1
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9     INFO(7)=ITER
      IF (INFO(8).NE.2) GOTO 998
      DO 999 I=1,NDIS
999   V(I)=VSAVE(I)
998   CONTINUE
      IF (IRANK.LT.NE) GOTO 93
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  TERMINATION AFTER ITMAX ITERATIONS
91    INFO(9)=1
      IF (KPRINT.GT.1) WRITE(UPR,60901)
      RETURN
C  RESTRICTED MONOTONICITY TEST VIOLATED
92    INFO(9)=2
      IF (KPRINT.GT.1) WRITE(UPR,60902)
      SUMYH=SUMY
      SUMYQH=SUMYQ
      RETURN
C  RANK DEFICIENCY OF JACOBIAN
93    INFO(9)=3
      IF (KPRINT.GT.1) WRITE(UPR,60903)
      RETURN
C  TAUMIN OR TAUMAX CROSSED
94    INFO(9)=4
      IF (KPRINT.GT.1) WRITE(UPR,60904)
      RETURN
C  DISCRETIZATION OF BOUNDARY VALUE PROBLEM FAILED
95    INFO(9)=5
      IF (KPRINT.GT.1) WRITE(UPR,60905)
      RETURN
C  NON-ZERO RESIDUAL IN SEPARATED BOUNDARY CONDITION
96    INFO(9)=6
      IF (KPRINT.GT.1) WRITE(UPR,60906)
      RETURN
C  CHANGE IN PROJECTIONS
97    INFO(9)=7
      IF (KPRINT.GT.1) WRITE(UPR,60907)
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60002 FORMAT(5X,I2,4X,D10.4,3X,D10.4,4X,I2,4X,I2,
     &   4X,D8.2,3X,D8.2,3X,20I3)
60003 FORMAT(4X,'ITER',5X,'LEVELX',6X,'LEVELXQ',6X,'NEW',2X,'RANK',
     &   5X,'COND',7X,'SENS',9X,'PIVOTS')
60901 FORMAT(' TERMINATION AFTER ITMAX ITERATIONS')
60902 FORMAT(' TERMINATION SINCE RESTRICTED MONOTONICITY TEST WAS ',
     &   'VIOLATED')
60903 FORMAT(' RANK DEFICIENT SOLUTION OBTAINED')
60904 FORMAT(' TERMINATION SINCE TAUMIN OR TAUMAX WAS CROSSED')
60905 FORMAT(' TERMINATION SINCE DISCRETIZATION OF BOUNDARY VALUE',
     &   ' PROBLEM FAILED')
60906 FORMAT(' TERMINATION SINCE NON-ZERO RESIDUAL IN SEPARATED',
     &   ' BOUNDARY CONDITION WAS FOUND')
60907 FORMAT(' TERMINATION DUE TO CHANGE IN PROJECTIONS')
      END
C
C
      SUBROUTINE MUTURN (N,N1,IPIVA,EPS,INFO,Y,Y0,Y1,YH,DYT,DYT0,DYT1,V,
     &   YW,IPIV,FY,TAUMIN,TAUMAX,D,E,EH,YA,DY,DYQ,W,F,FH,U,
     &  M,M1,NDIS,BC,ES,IROW,ICOL,ICOLB,IA,IB,T,IVPSOL,XU,HH,X1Q,XM,DX1,
     &  WM,RQ,G,AQ,BQ,P,DE,HSTART,XTHR,RELDIF,NE,NE1,NB,IPIVAQ,UNRM,
     &  US,DU,DHH,DR,DDX,RF,NEA,NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,
     &  TOLF,TOLJ)
C
C  SUBROUTINE  MUTURN TO BE USED WITH ROUTINE MULCON
C
C  DETERMINATION OF TURNING POINTS
C  INTERVAL METHOD WITH CUBIC HERMITE INTERPOLATION AS PREDICTOR
C  AND GAUSS-NEWTON METHOD AS CORRECTOR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(NDIS),Y0(NDIS),Y1(NDIS),YH(NDIS)
      DOUBLE PRECISION DYT(NDIS),DYT0(NDIS),DYT1(NDIS)
      DOUBLE PRECISION V(NDIS),YW(NDIS)
      DOUBLE PRECISION D(N1),E(N,N1),EH(N1,N1)
      DOUBLE PRECISION YA(NDIS),DY(NDIS),DYQ(NDIS),W(NDIS)
      DOUBLE PRECISION F(N),FH(N),U(N)
      DOUBLE PRECISION ES(N,N1),T(M),XU(N,M1),HH(N,M1),X1Q(N),XM(N),
     $DX1(N1),WM(N,N),RQ(N),G(N,N,M1),AQ(N,N),BQ(N,N),P(N,M),DE(N)
      DOUBLE PRECISION US(N),DU(N),DHH(N,M1),DR(N),DDX(N,M),RF(M)
      DOUBLE PRECISION HOF(M1),HOJ(M1)
      INTEGER KOF(M1),KOJ(M1)
      INTEGER INFO(9),IPIV(N1)
      INTEGER IROW(N),ICOL(N),ICOLB(N),IA(N,N),IB(N,N)
      INTEGER UPR,UDIAG
      EXTERNAL FY,DFY,BC,IVPSOL
      COMMON /MUUNIT/ UPR,UDIAG
      COMMON /MUVALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C  PREPARATIONS
      KPRINT=INFO(1)
      INFO(9)=-1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INFO(4).EQ.0 .OR. INFO(4).EQ.2) GOTO 9990
      IPIVH=IPIVAQ
      IF (IPIVA.EQ.NE1) IPIVH=ICOL(IPIV(NE1))
      DO 1000 I=1,NDIS
      Y0(I)=Y(I)
      DYT1(I)=DYT(I)
1000  CONTINUE
      F0P=DYT0(NDIS)/DYT0(IPIVH)
      F1P=DYT1(NDIS)/DYT1(IPIVH)
      INFO(1)=0
      ITER=0
      X0=Y0(IPIVH)/YW(IPIVH)
      X1=Y1(IPIVH)/YW(IPIVH)
      H=X1-X0
      DEL=DABS(H)
C-----------------------------------------------------------------------
C  ITERATION LOOP
2000  CONTINUE
C  HERMITE INTERPOLATION
      F0=Y0(NDIS)/YW(NDIS)
      F1=Y1(NDIS)/YW(NDIS)
      A=6.D0*(F0-F1)+3.D0*H*(F0P+F1P)
      B=-6.D0*(F0-F1)-2.D0*H*(2.D0*F0P+F1P)
      C=H*F0P
      Z=B**2-4.D0*A*C
      Z=DSIGN(DSQRT(Z),-B)
      Z=(-B+Z)/(2.D0*A)
      IF (Z.LT.0.D0 .OR. Z.GT.1.D0) Z=(C/A)/Z
      YH(IPIVH)=Y0(IPIVH)+H*Z*YW(IPIVH)
      DO 2020 I=1,NDIS
      IF (I.EQ.IPIVH) GOTO 2010
      F0H=Y0(I)/YW(I)
      F1H=Y1(I)/YW(I)
      F0PH=DYT0(I)/DYT0(IPIVH)
      F1PH=DYT1(I)/DYT1(IPIVH)
      A=2.D0*(F0H-F1H)+H*(F0PH+F1PH)
      B=-3.D0*(F0H-F1H)-H*(2.D0*F0PH+F1PH)
      C=H*F0PH
      R=F0H+Z*(C+Z*(B+Z*A))
      YH(I)=R*YW(I)
2010  CONTINUE
2020  CONTINUE
C  ITERATION BACK TO CONTINUATION PATH
      ITER=ITER+1
      IF (ITER.GT.ITMAX) GOTO 9910
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      INFO(8)=2
      CALL MUGNHO(N,N1,M,M1,NDIS,YH,FY,BC,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,E,EH,ES,YA,DY,DYQ,W,
     &  IROW,ICOL,ICOLB,IA,IB,T,IVPSOL,XU,HH,X1Q,XM,DX1,WM,RQ,G,AQ,BQ,P,
     &   DE,F,FH,U,HSTART,XTHR,RELDIF,NE,NB,US,DU,DHH,DR,DDX,RF,NEA,
     &   NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,TOLF,TOLJ)
      IERR=INFO(9)
      GOTO (9920,9920,9920,9920,9920,9940,9950),IERR
      CALL MUTAND(N,N1,NDIS,FY,YH,F,FH,YW,E,COND,DETF,IPIV,EH,
     &   V,D,IRANK,BC,IVPSOL,M,M1,T,X1Q,XM,RQ,XU,HH,W,HSTART,TOLJ,WM,
     &INFO,AQ,BQ,P,RELDIF,NE,NB,IA,IB,DE,IROW,ICOL,ICOLB,XTHR,G,ES,UNRM,
     &   NEA,DFY,HOJ,KOJ,INITJ)
      IERR=INFO(9)
      GOTO (9920,9920,9940,9950),IERR
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      IF (DABS(YH(IPIVH)-Y0(IPIVH)).GT.DEL*YW(IPIVH) .OR.
     &    DABS(YH(IPIVH)-Y1(IPIVH)).GT.DEL*YW(IPIVH)) GOTO 9930
      IPIVS=IPIV(NE1)
      IRANKB=NB
      DO 1201 I=1,NE
1201  W(IPIV(I))=V(I)
      W(IPIV(NE1))=-1.D0
      IF (IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUREFN(N,N1,NE,NE1,ES,W,DU,US,IRANKB,V,IRANK,D,IPIV,
     &   WM,EPH,NE1,EH,EPS,E,2,INFO)
      IF (INFO(9).NE.0) GOTO 9920
      DO 1202 I=1,NE1
1202  V(I)=W(IPIV(I))
      DO 2030 I=1,N
      W(I)=0.D0
2030  CONTINUE
      DO 2031 I=1,NE
      IF (IPIV(I).NE.NE1)
     &   W(ICOL(IPIV(I)))=-V(I)/D(NE1)*YW(ICOL(IPIV(I)))
      IF (IPIV(I).EQ.NE1) W(NDIS)=-V(I)/D(NE1)*YW(NDIS)
2031  CONTINUE
      IF (IPIVS.NE.NE1) W(ICOL(IPIVS))=1.D0/D(NE1)*YW(ICOL(IPIVS))
      IF (IPIVS.EQ.NE1) W(NDIS)=1.D0/D(NE1)*YW(NDIS)
      DO 40000 I=1,N
      DO 40000 J=1,M1
40000 HH(I,J)=0.D0
      CALL MURECU(N,M,M1,1,HH,G,P,W,W(NDIS),V,WM)
      IF (NYMAX.GT.0 .AND. NE.GT.0)
     &   CALL MUSWP(N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,W,DHH,
     & RQ,AQ,BQ,DR,DU,DE,V,WM,IROW,DX1,NE,NE1,IRANK,IRANKB,IPIV,D,E,EH,
     &   ICOL,YW,YW(NDIS),DDX,RF,W(NDIS),2,RELDIF,TOLF,TOLJ,INFO)
      IF (INFO(9).NE.0) GOTO 9920
      IF (NB.EQ.0) GOTO 1213
      DO 1212 I=1,NB
      W(N*M1+ICOLB(I))=0.D0
1212  CONTINUE
1213  CONTINUE
      DO 2032 I=1,NDIS
      W(I)=W(I)/YW(I)
2032  CONTINUE
      FQP=W(NDIS)/W(IPIVH)
C  NEW INCLUSION INTERVAL
      IF (F0P*FQP.GT.0.D0) GOTO 2200
      DQ=DABS(YH(IPIVH)-Y1(IPIVH))
      F1P=FQP
      DO 2110 I=1,NDIS
      Y1(I)=YH(I)
      DYT1(I)=W(I)
2110  CONTINUE
      GOTO 2300
2200  CONTINUE
      DQ=DABS(YH(IPIVH)-Y0(IPIVH))
      F0P=FQP
      DO 2210 I=1,NDIS
      Y0(I)=YH(I)
      DYT0(I)=W(I)
2210  CONTINUE
2300  CONTINUE
      X0=Y0(IPIVH)/YW(IPIVH)
      X1=Y1(IPIVH)/YW(IPIVH)
      H=X1-X0
      DEL=DABS(H)
      IF (KPRINT.GT.0 .AND. ITER.EQ.1) WRITE(UPR,60101)
      IF (KPRINT.GT.0) WRITE(UPR,60102) ITER,YH(IPIVH),DEL,FQP,INFO(7)
      IF (FQP.EQ.0.D0) GOTO 9000
      IF (DEL.LT.EPS) GOTO 9000
      IF (DQ/YW(IPIVH).LT.EPS) GOTO 9000
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9000  CONTINUE
      INFO(9)=0
      IF (KPRINT.GT.0) WRITE(UPR,60103) ITER
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(NDIS),UNRM
      DO 9001 J=1,M
      IF (KPRINT.GT.0) WRITE(UPR,60105) J,T(J),(YH(I+(J-1)*N),I=1,N)
9001  CONTINUE
      GOTO 9990
C  FAIL EXIT
C  MORE THAN ITMAX ITERATIONS
9910  CONTINUE
      INFO(9)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      GOTO 9990
C  GAUSS-NEWTON METHOD FAILED
9920  CONTINUE
      INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      GOTO 9990
C  ITERATION STRATEGY FAILED
9930  CONTINUE
      INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60903)
      GOTO 9990
C  NON-ZERO RESIDUAL IN SEPARATED BOUNADRY CONDITIONS
9940  CONTINUE
      INFO(9)=4
      IF (KPRINT.GT.0) WRITE(UPR,60904)
      GOTO 9990
C  CHANGE IN PROJECTIONS
9950  CONTINUE
      INFO(9)=5
      IF (KPRINT.GT.0) WRITE(UPR,60905)
C  RETURN TO CALLING ROUTINE
9990  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      INFO(1)=KPRINT
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(' ',132('-'))
60002 FORMAT(' POSSIBLY TURNING POINT DETECTED')
60003 FORMAT(/)
60101 FORMAT(/4X,'IGNC',7X,'YHIT',8X,'DELTAX',8X,'DERIV',6X,'ITER')
60102 FORMAT(5X,I2,4X,D11.4,3X,D10.4,3X,D11.4,4X,I2)
60103 FORMAT(/' MUTURN REQUIRED',I3,' GAUSS-NEWTON-CALLS'/)
60104 FORMAT(' TAU=',D14.6,'      NORM=',D14.6/)
60105 FORMAT('      NODE',I3,'     T=',D14.6/
     &   25X,'X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60201 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,4(D14.6,6X)/))
60901 FORMAT(/' MUTURN TERMINATED AFTER ITMAX ITERATIONS')
60902 FORMAT(/' MUTURN TERMINATED SINCE GAUSS-NEWTON METHOD FAILED')
60903 FORMAT(/' MUTURN TERMINATED SINCE ITERATION STRATEGY FAILED')
60904 FORMAT(/' MUTURN TERMINATED DUE TO NON-ZERO RESIDUAL IN',
     &   ' SEPARATED BOUNDARY CONDITIONS')
60905 FORMAT(/' MUTURN TERMINATED DUE TO CHANGE IN PROJECTIONS')
      END
C
C
      SUBROUTINE MUBIFC (N,N1,IPIVA,EPS,INFO,Y,Y0,Y1,YH,V,
     &  YW,IPIV,FY,TAUMIN,TAUMAX,D,E,EH,YA,DY,DYQ,W,F,FH,U,DETFS,DETFAH,
     &  M,M1,NDIS,BC,ES,IROW,ICOL,ICOLB,IA,IB,T,IVPSOL,XU,HH,X1Q,XM,DX1,
     &  WM,RQ,G,AQ,BQ,P,DE,HSTART,XTHR,RELDIF,NE,NE1,NB,IPIVAQ,UNRM,
     &  US,DU,DHH,DR,DDX,RF,NEA,NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,
     &  TOLF,TOLJ)
C
C  SUBROUTINE  MUBIFC  TO BE USED WITH ROUTINE MULCON
C
C  DETERMINATION OF BIFURCATION POINTS
C  INTERVAL METHOD WITH BRENT ALGORITHM IN CONNECTION WITH
C  LINEAR INTERPOLATION AS PREDICTOR
C  AND GAUSS-NEWTON METHOD AS CORRECTOR
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(NDIS),Y0(NDIS),Y1(NDIS),YH(NDIS)
      DOUBLE PRECISION V(NDIS),YW(NDIS)
      DOUBLE PRECISION D(N1),E(N,N1),EH(N1,N1)
      DOUBLE PRECISION YA(NDIS),DY(NDIS),DYQ(NDIS),W(NDIS)
      DOUBLE PRECISION F(N),FH(N),U(N)
      DOUBLE PRECISION ES(N,N1),T(M),XU(N,M1),HH(N,M1),X1Q(N),XM(N),
     $DX1(N1),WM(N,N1),RQ(N),G(N,N,M1),AQ(N,N),BQ(N,N),P(N,M),DE(N)
      DOUBLE PRECISION US(N),DU(N),DHH(N,M1),DR(N),DDX(N,M),RF(M)
      DOUBLE PRECISION HOF(M1),HOJ(M1)
      INTEGER KOF(M1),KOJ(M1)
      INTEGER INFO(9),IPIV(N1)
      INTEGER IROW(N),ICOL(N),ICOLB(N),IA(N,N),IB(N,N)
      INTEGER UPR,UDIAG
      EXTERNAL FY,DFY,BC,IVPSOL
      COMMON /MUUNIT/ UPR,UDIAG
      COMMON /MUVALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C  PREPARATIONS
      KPRINT=INFO(1)
      INFO(9)=-1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      IF (INFO(4).LT.2) GOTO 9990
      DETF=DETFS
      DETFA=DETFAH
      DO 1000 I=1,NDIS
      Y0(I)=Y(I)
1000  CONTINUE
      INFO(1)=0
      ITER=0
      X0=Y0(IPIVAQ)/YW(IPIVAQ)
      X1=Y1(IPIVAQ)/YW(IPIVAQ)
      H=X1-X0
      DEL=DABS(H)
      AZ=X0
      FA=DETFA
      BZ=X1
      FB=DETF
C-----------------------------------------------------------------------
C  ITERATION LOOP
C  BRENT ALGORITHM
2000  CONTINUE
      CZ=AZ
      FC=FA
      DZ=BZ-CZ
      EZ=DZ
2010  IF (DABS(FC).GE.DABS(FB)) GOTO 2015
      AZ=BZ
      BZ=CZ
      CZ=AZ
      FA=FB
      FB=FC
      FC=FA
2015  CONTINUE
      RM=(CZ-BZ)*0.5D0
      IF (DABS(CZ-BZ).LE.EPS .AND. ITER.GT.0) GOTO 9000
      IF (DABS(EZ).LT.EPS) GOTO 2030
      IF (DABS(FA).LE.DABS(FB)) GOTO 2030
      SZ=FB/FA
      IF (AZ.NE.CZ) GOTO 2020
      PZ=(CZ-BZ)*SZ
      QZ=1.D0-SZ
      GOTO 2025
2020  QZ=FA/FC
      RZ=FB/FC
      RZ1=RZ-1.D0
      PZ=SZ*((CZ-BZ)*QZ*(QZ-RZ)-(BZ-AZ)*RZ1)
      QZ=(QZ-1.D0)*RZ1*(SZ-1.D0)
2025  IF (PZ.GT.0.D0) QZ=-QZ
      IF (PZ.LT.0.D0) PZ=-PZ
      SZ=EZ
      EZ=DZ
      IF (PZ+PZ.GE.3.D0*RM*QZ) GOTO 2030
      IF (PZ+PZ.GE.DABS(SZ*QZ)) GOTO 2030
      DZ=PZ/QZ
      GOTO 2035
2030  EZ=RM
      DZ=EZ
2035  AZ=BZ
      FA=FB
      TEMP=DZ
      IF (DABS(TEMP).LE.0.5D0*EPS) TEMP=DSIGN(0.5D0*EPS,RM)
      BZ=BZ+TEMP
      SZ=BZ
      Z=(SZ-AZ)/(CZ-AZ)
      IF (FA.EQ.DETF) Z=1.D0-Z
C  LINEAR INTERPOLATION
      DO 2080 I=1,NDIS
      F0H=Y0(I)/YW(I)
      F1H=Y1(I)/YW(I)
      R=F0H+Z*(F1H-F0H)
      YH(I)=R*YW(I)
2080  CONTINUE
C  ITERATION BACK TO CONTINUATION PATH
      ITER=ITER+1
      IF (ITER.GT.ITMAX) GOTO 9910
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      INFO(8)=2
      CALL MUGNHO(N,N1,M,M1,NDIS,YH,FY,BC,EPS,INFO,SUMYH,SUMYQH,
     &   TAUMIN,TAUMAX,IPIV,V,D,YW,E,EH,ES,YA,DY,DYQ,W,
     &  IROW,ICOL,ICOLB,IA,IB,T,IVPSOL,XU,HH,X1Q,XM,DX1,WM,RQ,G,AQ,BQ,P,
     &   DE,F,FH,U,HSTART,XTHR,RELDIF,NE,NB,US,DU,DHH,DR,DDX,RF,NEA,
     &   NYMAX,DFY,HOF,HOJ,KOF,KOJ,INITF,INITJ,TOLF,TOLJ)
      IERR=INFO(9)
      INFO(9)=-1
      GOTO (9920,9920,9010,9920,9920,9940,9950),IERR
      CALL MUTAND(N,N1,NDIS,FY,YH,F,FH,YW,E,COND,DETFH,IPIV,EH,
     &   V,D,IRANK,BC,IVPSOL,M,M1,T,X1Q,XM,RQ,XU,HH,W,HSTART,TOLJ,WM,
     &INFO,AQ,BQ,P,RELDIF,NE,NB,IA,IB,DE,IROW,ICOL,ICOLB,XTHR,G,ES,UNRM,
     &   NEA,DFY,HOJ,KOJ,INITJ)
      IERR=INFO(9)
      GOTO (9920,9920,9940,9950),IERR
      IPIVS=IPIV(NE1)
      IF (KPRINT.GT.2) WRITE(UPR,60201) YH
      IF (DABS(YH(IPIVAQ)-Y0(IPIVAQ)).GT.DEL*YW(IPIVAQ) .OR.
     &    DABS(YH(IPIVAQ)-Y1(IPIVAQ)).GT.DEL*YW(IPIVAQ)) GOTO 9930
      IRANKB=NB
      DO 1201 I=1,NE
1201  W(IPIV(I))=V(I)
      W(IPIV(NE1))=-1.D0
      IF (IRANK.EQ.NE .AND. NE.GT.0)
     &   CALL MUREFN(N,N1,NE,NE1,ES,W,DU,US,IRANKB,V,IRANK,D,IPIV,
     &   WM,EPH,NE1,EH,EPS,E,2,INFO)
      IF (INFO(9).NE.0) GOTO 9920
      DO 1202 I=1,NE1
1202  V(I)=W(IPIV(I))
      DO 2033 I=1,N
      W(I)=0.D0
2033  CONTINUE
      DO 2031 I=1,NE
      IF (IPIV(I).NE.NE1)
     &   W(ICOL(IPIV(I)))=-V(I)/D(NE1)*YW(ICOL(IPIV(I)))
      IF (IPIV(I).EQ.NE1) W(NDIS)=-V(I)/D(NE1)*YW(NDIS)
2031  CONTINUE
      IF (IPIVS.NE.NE1) W(ICOL(IPIVS))=1.D0/D(NE1)*YW(ICOL(IPIVS))
      IF (IPIVS.EQ.NE1) W(NDIS)=1.D0/D(NE1)*YW(NDIS)
      CALL MURECU(N,M,M1,1,HH,G,P,W,W(NDIS),V,WM)
      DO 40000 I=1,N
      DO 40000 J=1,M1
40000 HH(I,J)=0.D0
      IF (NYMAX.GT.0 .AND. NE.GT.0)
     &   CALL MUSWP(N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,W,DHH,
     & RQ,AQ,BQ,DR,DU,DE,V,WM,IROW,DX1,NE,NE1,IRANK,IRANKB,IPIV,D,E,EH,
     &   ICOL,YW,YW(NDIS),DDX,RF,W(NDIS),2,RELDIF,TOLF,TOLJ,INFO)
      IF (INFO(9).NE.0) GOTO 9920
      IF (NB.EQ.0) GOTO 1213
      DO 1212 I=1,NB
      W(N*M1+ICOLB(I))=0.D0
1212  CONTINUE
1213  CONTINUE
      DO 2032 I=1,NDIS
      W(I)=W(I)/YW(I)
2032  CONTINUE
      IF (IPIV(NE1).NE.IPIVA)
     &   CALL MUDETH(N,N1,NE,NE1,E,WM,D,IPIVA,IPIV,DETFH)
C  NEW INCLUSION INTERVAL
      IF (DETFA*DETFH.GT.0.D0) GOTO 2200
      DETF=DETFH
      DO 2110 I=1,NDIS
      Y1(I)=YH(I)
2110  CONTINUE
      GOTO 2300
2200  CONTINUE
      DETFA=DETFH
      DO 2210 I=1,NDIS
      Y0(I)=YH(I)
2210  CONTINUE
2300  CONTINUE
      X0=Y0(IPIVAQ)/YW(IPIVAQ)
      X1=Y1(IPIVAQ)/YW(IPIVAQ)
      H=X1-X0
      DEL=DABS(H)
      IF (KPRINT.GT.0 .AND. ITER.EQ.1) WRITE(UPR,60101)
      IF (KPRINT.GT.0)
     &   WRITE(UPR,60102) ITER,YH(IPIVAQ),DEL,DETFH,INFO(7)
      FB=DETFH
      IF (FB*FC.LE.0.D0) GOTO 2010
      GOTO 2000
C-----------------------------------------------------------------------
C  SOLUTION EXIT
9000  CONTINUE
      INFO(9)=0
9010  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60103) ITER
      IF (KPRINT.GT.0) WRITE(UPR,60104) YH(NDIS),UNRM
      DO 9001 J=1,M
      IF (KPRINT.GT.0) WRITE(UPR,60105) J,T(J),(YH(I+(J-1)*N),I=1,N)
9001  CONTINUE
      IF (KPRINT.GT.0 .AND. IRANK.LT.N) WRITE(UPR,60109)
      GOTO 9990
C  FAIL EXIT
C  MORE THAN ITMAX ITERATIONS
9910  CONTINUE
      INFO(9)=1
      IF (KPRINT.GT.0) WRITE(UPR,60901)
      GOTO 9990
C  GAUSS-NEWTON METHOD FAILED
9920  CONTINUE
      INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60902)
      GOTO 9990
C  ITERATION STRATEGY FAILED
9930  CONTINUE
      INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60903)
      GOTO 9990
C  NON-ZERO RESIDUAL IN SEPARATED BOUNADRY CONDITIONS
9940  CONTINUE
      INFO(9)=4
      IF (KPRINT.GT.0) WRITE(UPR,60904)
      GOTO 9990
C  CHANGE IN PROJECTIONS
9950  CONTINUE
      INFO(9)=5
      IF (KPRINT.GT.0) WRITE(UPR,60905)
C  RETURN TO CALLING ROUTINE
9990  CONTINUE
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      INFO(1)=KPRINT
      RETURN
C-----------------------------------------------------------------------
60001 FORMAT(' ',132('-'))
60002 FORMAT(' POSSIBLY BIFURCATION POINT DETECTED')
60003 FORMAT(/)
60101 FORMAT(/4X,'IGNC',7X,'YHIT',8X,'DELTAX',8X,'DETER',6X,'ITER')
60102 FORMAT(5X,I2,4X,D11.4,3X,D10.4,3X,D11.4,4X,I2)
60103 FORMAT(/' MUBIFC REQUIRED',I3,' GAUSS-NEWTON-CALLS'/)
60104 FORMAT(' TAU=',D14.6,'      NORM=',D14.6/)
60105 FORMAT('      NODE',I3,'     T=',D14.6/
     &   25X,'X=',4(D14.6,4X)/9(27X,4(D14.6,4X)/))
60109 FORMAT(/' NO TANGENT AT SOLUTION AVAILABLE')
60201 FORMAT(18X,'   Y=',5(D14.6,6X)/9(23X,5(D14.6,6X)/))
60901 FORMAT(/' MUBIFC TERMINATED AFTER ITMAX ITERATIONS')
60902 FORMAT(/' MUBIFC TERMINATED SINCE GAUSS-NEWTON METHOD FAILED')
60903 FORMAT(/' MUBIFC TERMINATED SINCE ITERATION STRATEGY FAILED')
60904 FORMAT(/' MUBIFC TERMINATED DUE TO NON-ZERO RESIDUAL IN',
     &   ' SEPARATED BOUNDARY CONDITIONS')
60905 FORMAT(/' MUBIFC TERMINATED DUE TO CHANGE IN PROJECTIONS')
      END
C
C
      SUBROUTINE MUDETH (N,N1,NE,NE1,A,AH,D,IPIVA,IPIV,DET)
C
C  SUBROUTINE  MUDETH  TO BE USED WITH ROUTINE MULCON
C
C  COMPUTATION OF DETERMINANT OF A HESSENBERG MATRIX
C  BY GIVENS ROTATIONS
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(N,N1),AH(N,N1),D(N1),IPIV(N1)
C  SEARCH OLD PIVOT
      DO 1010 I=1,NE1
      IQ=I
      IF (IPIV(I).EQ.IPIVA) GOTO 1020
1010  CONTINUE
1020  CONTINUE
C  PRODUCT OF FIRST IQ-1 DIAGONAL ELEMENTS
      DET=1.D0
      IF (IQ.EQ.1) GOTO 2011
      IEND=IQ-1
      DO 2010 I=1,IEND
      DET=DET*D(I)
2010  CONTINUE
2011  CONTINUE
C  N-IQ GIVENS TRANSFORMATIONS
      IF (IQ.EQ.NE1) GOTO 3050
      IBEG=IQ+1
      X1=A(IBEG-1,IBEG)
      IF (IBEG.GT.NE) GOTO 3040
      DO 3010 I=IBEG,NE1
      AH(IQ,I)=A(IQ,I)
3010  CONTINUE
      DO 3030 I=IBEG,NE
      Y1=D(I)
      R=DSQRT(X1**2+Y1**2)
      C=X1/R
      S=Y1/R
      DET=DET*(C*X1+S*Y1)
      I1=I+1
      DO 3020 J=I1,NE1
      X2=AH(I-1,J)
      Y2=A(I,J)
      AH(I,J)=-S*X2+C*Y2
3020  CONTINUE
      X1=AH(I,I1)
3030  CONTINUE
3040  CONTINUE
      DET=DET*X1
3050  CONTINUE
C  SIGN OF PIVOT VECTOR
      DO 4020 I=1,NE
      IF (I.EQ.IQ) GOTO 4020
      I1=I+1
      DO 4010 J=I1,NE1
      IF (J.EQ.IQ) GOTO 4010
      IF (IPIV(J).LT.IPIV(I)) DET=-DET
4010  CONTINUE
4020  CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MUTAND (N,N1,NDIS,FY,Y,F,FH,YW,E,COND,DET,PIVOT,EH,
     &   V,D,IRANK,BC,IVPSOL,M,M1,T,X1,XM,R,XU,HH,W,HSTART,TOLJ,WM,
     &  INFO,A,B,P,RELDIF,NE,NB,IA,IB,DE,IROW,ICOL,ICOLB,XTHR,G,ES,UNRM,
     &  NEA,DFY,HOJ,KOJ,INITJ)
C
C  SUBROUTINE  MUTAND  TO BE USED WITH ROUTINE MULCON
C
C  COMPUTATION OF TANGENT AND DETERMINANT AT GIVEN POINT OF
C  CONTINUATION PATH BY QR DECOMPOSITION
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION F(N),FH(N),E(N,N1)
      DOUBLE PRECISION EH(N1,N1),V(N1),D(N1),WM(N,N)
      DOUBLE PRECISION T(M),Y(NDIS),X1(N),XM(N),R(N),XU(N,M1),HH(N,M1),
     $W(NDIS),A(N,N),B(N,N),P(N,M),YW(NDIS),G(N,N,M1),ES(N,N1),DE(N)
      DOUBLE PRECISION HOJ(M1)
      INTEGER KOJ(M1)
      INTEGER INFO(9),IROW(N),ICOL(N),ICOLB(N),IA(N,N),IB(N,N)
      INTEGER PIVOT(N1)
      EXTERNAL FY,DFY,BC,IVPSOL
      COMMON /COUNT/ IFCTEV,ITER,IDECS,ISOLS,ITRAJ
      COMMON /MACHIN/ EPMACH,SMALL
C  COMPUTE NEW TANGENT AND DETERMINANT
      CALL MURESI(FY,BC,IVPSOL,N,N1,M,M1,T,Y,Y(NDIS),X1,XM,R,
     &   XU,HH,W,HSTART,TOLJ,INFO,UNRM,YW,YW(NDIS),TFAIL,G,P,1,W,DFY,
     &   HOJ,KOJ,INITJ)
      IF (INFO(9).NE.0) GOTO 92
      LYM=N*M1+1
      CALL MUDERA(BC,N,M,Y,Y(LYM),Y(NDIS),R,W,A,B,P(1,M),
     &   YW,YW(NDIS),RELDIF)
      CALL MUCHAB(N,N1,M,NE,NB,Y,A,B,IA,IB,P(1,M),R,DE,YW,YW(NDIS),
     &   IROW,ICOL,ICOLB,XTHR,INFO)
      IF (INFO(9).NE.0) GOTO 93
      IF (INFO(8).EQ.1) NEA=NE
      IF (NE.NE.NEA) GOTO 94
      NE1=NE+1
      CALL MUCDNS(N,N1,M,M1,NE,NE1,G,A,B,YW,DE,E,ES,W,IROW,ICOL,
     &   P,YW(NDIS),WM,V)
      DO 7410 I=1,NE
      DO 7410 J=1,NE1
7410  ES(I,J)=E(I,J)
      COND=1.D0/EPMACH
      IRANK=NE
      CALL MUDECC(E,N,N1,NB,NE,NE1,IRANK,COND,D,PIVOT,0,EH,V)
      IDECS=IDECS+1
      IF (IRANK.LT.NE) GOTO 91
      DET=1.D0
      DO 1020 I=1,NE
      DET=DET*D(I)
1020  CONTINUE
      NM1=NE-1
      IF (NM1.EQ.0) GOTO 1031
      DO 1030 I=1,NM1
      I1=I+1
      DO 1030 J=I1,NE
      IF (PIVOT(J).LT.PIVOT(I)) DET=-DET
1030  CONTINUE
1031  CONTINUE
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  RANK DEFICIENCY
91    CONTINUE
      INFO(9)=1
      RETURN
C  BOUNDARY VALUE DISCRETIZATION FAILED
92    CONTINUE
      INFO(9)=2
      RETURN
C  NON-ZERO RESIDUAL IN SEPARATED BOUNDARY CONDITIONS
93    CONTINUE
      INFO(9)=3
      RETURN
C  CHANGE IN PROJECTIONS
94    CONTINUE
      INFO(9)=4
      RETURN
      END
C
C
      BLOCK DATA MUBLOK
C
C  BLOCK DATA  TO BE USED WITH ROUTINE MULCON
C
C    COMMON /MUUNIT/ UPR,UDIAG,UBIF
C      INPUT/OUTPUT UNITS
C      - UPR     PRINT UNIT                            OUTPUT
C      - UDIAG   UNIT FOR PLOT INFORMATION             OUTPUT
C                (IN CONNECTION WITH ROUTINE PLTHM)
C
C    COMMON /MUVALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C      LIMIT VALUES IN CONTINUATION PROCESS
C      - SMALL1  USED AS INITIAL STEPLENGTH ESTIMATE AND AS MAXIMUM
C                PERMITTED STEPSIZE IN NUMERICAL DIFFERENTIATION
C      - THMAX   MAXIMUM PERMITTED RATIO IN RESTRICTED MONOTONICITY TEST
C      - THR     THRESHOLD IN NEARLY LINEAR CASE
C      - ITMAX   MAXIMUM PERMITTED NUMBER OF ITERATIONS
C                IN INTERVAL METHODS
C      - IGNMAX  MAXIMUM PERMITTED NUMBER OF GAUSS-NEWTON ITERATIONS
C                PER CORRECTOR CALL
C
C    COMMON /MUMETH/ FACTOR,IEXTR,IJACM
C      METHOD DESCRIBING PARAMETERS
C      - FACTOR  FACTOR FOR EXTRAPOLATED STEPLENGTH BOUND
C      - IEXTR   METHOD FOR EXTRAPOLATED STEPLENGTH BOUND
C                  0  NO EXTRAPOLATED STEPLENGTH BOUND
C                  1  WITH EXTRAPOLATED STEPLENGTH BOUND
C      - IJACM   METHOD FOR UPDATING JACOBIAN IN GAUSS-NEWTON METHOD
C                  0  KEEP FIRST JACOBIAN
C                  1  BROYDEN UPDATES OF JACOBIAN
C                  2  NEW JACOBIAN IN EACH ITERATION
C                  3  AS 1  BUT IN A COMPUTING TIME SAVING WAY
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER UPR,UDIAG
      COMMON /MUUNIT/ UPR,UDIAG
      COMMON /MUVALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
      COMMON /MUMETH/ FACTOR,IEXTR,IJACM
      DATA UPR/6/,UDIAG/2/
      DATA SMALL1/0.01D0/,THMAX/0.25D0/,THR/0.1D0/,ITMAX/20/,IGNMAX/10/
      DATA FACTOR/1.1D0/,IEXTR/0/,IJACM/1/
      END
C
C
      SUBROUTINE MURHS1 (N,NE,M1,JIN,HH,R,B,G,U,DE,V,BG,IROW)
C  COMPUTATION OF CONDENSED RIGHT-HAND SIDE U(NE)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION G(N,N,M1),HH(N,M1),B(N,N),BG(N,N),U(N),DE(N),V(N),R(N)
      INTEGER IROW(N)
      DO 100 I=1,NE
      IR=IROW(I)
100   U(I)=DE(IR)*R(IR)
      IF(JIN.GT.M1) RETURN
      DO 110 I=1,NE
      IR=IROW(I)
      S=U(I)
      DO 111 K=1,N
      TH=DE(IR)*B(IR,K)
      BG(I,K)=TH
111   S=S+TH*HH(K,M1)
110   U(I)=S
      IF(M1.EQ.1.OR.JIN.EQ.M1) RETURN
C
      M2=M1-1
      DO 200 JJ=JIN,M2
      J=M2+JIN-JJ
      J1=J+1
      DO 200 I=1,NE
      DO 210 K=1,N
      S=0.D0
      DO 211 L=1,N
211   S=S+BG(I,L)*G(L,K,J1)
210   V(K)=S
      S=U(I)
      DO 220 K=1,N
      S=S+V(K)*HH(K,J)
220   BG(I,K)=V(K)
200   U(I)=S
      RETURN
      END
C
C
      SUBROUTINE MURECU (N,M,M1,JIN,HH,G,P,DX,DP,U,V)
C  RECURSIVE SOLUTION OF M1 LINEAR (N,N)-SYSTEMS
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION G(N,N,M1),DX(N,M),HH(N,M1),P(N,M1),U(N),V(N)
      DO 10 I=1,N
10    U(I)=DX(I,1)
      DO 100 J=1,M1
      J1=J+1
      DO 110 I=1,N
      S=0.D0
      IF (J.GE.JIN) S=HH(I,J)
      S=S+P(I,J)*DP
      DO 111 K=1,N
111   S=S+G(I,K,J)*U(K)
      V(I)=S
110   DX(I,J1)=S
      DO 120 I=1,N
120   U(I)=V(I)
100   CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MUDERA (BC,N,M,X1,XM,TAU,R,RH,A,B,P,XW,TAUW,RELDIF)
C
C  DIFFERENCE APPROX. OF BOUNDARY DERIVATIVE MATRICES A(N,N) AND B(N,N)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION A(N,N),B(N,N),P(N),R(N),RH(N),X1(N),XM(N),XW(N,M)
      EXTERNAL BC
      DO 1 K=1,N
      XH=X1(K)
      S=RELDIF*XW(K,1)
      IF (XH.LT.0.D0) S=-S
      X1(K)=XH+S
      CALL BC(X1,XM,TAU,RH)
      X1(K)=XH
      S=1.D0/S
      DO 11 I=1,N
11    A(I,K)=(RH(I)-R(I))*S
      XH=XM(K)
      S=RELDIF*XW(K,M)
      IF (XH.LT.0.D0) S=-S
      XM(K)=XH+S
      CALL BC(X1,XM,TAU,RH)
      XM(K)=XH
      S=1.D0/S
      DO 12 I=1,N
12    B(I,K)=(RH(I)-R(I))*S
1     CONTINUE
      XH=TAU
      S=RELDIF*TAUW
      IF (XH.LT.0.D0) S=-S
      TAU=XH+S
      CALL BC(X1,XM,TAU,RH)
      TAU=XH
      S=1.D0/S
      DO 13 I=1,N
13    P(I)=(RH(I)-R(I))*S
      RETURN
      END
C
C
      SUBROUTINE MUDERG (FY,N,N1,M,M1,T,X,TAU,XU,XW,TAUW,XJ,TJ,G,P,
     &   IVPSOL,HSTART,TOL,INFO,W,DFY,HOS,KOS,INIT)
C
C  DIFFERENCE APPROXIMATION OF WRONSKIAN MATRICES G(1),..,G(M1)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION T(M),X(N,M),G(N,N,M1),P(N,M1),XW(N,M),XU(N,M1),XJ(N1)
      DIMENSION W(N),HOS(M1)
      INTEGER INFO(9),KOS(M1)
      INTEGER UPR,UDIAG
      COMMON /COUNT/ IFCTEV,ITER,IDECS,ISOLS,ITRAJ
      COMMON /MUUNIT/ UPR,UDIAG
      EXTERNAL FY,DFY,IVPSOL
      KPRINT=INFO(1)
      HSAVE=HSTART
      J=1
1     J1=J+1
      TJA=T(J)
      TJ1=T(J1)
      H=HSAVE
      HMAX=DABS(TJ1-TJA)
      TJ=TJA
      DO 11 K=1,N
11    XJ(K)=X(K,J)
      XJ(N1)=TAU
      UNRM=0.D0
      NVAR=N**2+2*N+1
      CALL IVPSOL(N,NVAR,FY,DFY,TJ,XJ,TJ1,XW(1,J),XW(1,J1),TAUW,
     &   TOL,HMAX,H,UNRM,G(1,1,J),P(1,J),1,W,HOS(J),KOS(J),INIT)
      IF (H.EQ.0.D0) GOTO 9
      HSAVE=H
      J=J1
      IF (J.LT.M) GOTO 1
      ITRAJ=ITRAJ+1
      INIT=1
      INFO(9)=0
      RETURN
C  ERROR RETURN
C  SINGULAR TRAJECTORY
9     INFO(9)=J
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      RETURN
60001 FORMAT(' TERMINATION AT SINGULAR TRAJECTORY')
      END
C
C
      SUBROUTINE MURK1G (N,N1,M,M1,XW,TAUW,DX,DTAU,HH,DXJ,G,P)
C  RANK-1 UPDATES OF WRONSKIAN MATRICES G(1),...,G(M1)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION G(N,N,M1),P(N,M),DX(N,M),XW(N,M),HH(N,M1),DXJ(N1)
      DO 100 J=1,M1
      DNM=0.D0
      DO 110 I=1,N
      T=DX(I,J)/XW(I,J)
      DXJ(I)=T/XW(I,J)
110   DNM=DNM+T*T
      T=DTAU/TAUW
      DXJ(N1)=T/TAUW
      DNM=DNM+T*T
      IF (DNM.EQ.0.D0) GOTO 100
      DO 120 K=1,N
      T=DXJ(K)/DNM
      DO 120 I=1,N
      S=G(I,K,J)
      IF (S.NE.0.D0) G(I,K,J)=S+T*HH(I,J)
120   CONTINUE
      T=DXJ(N1)/DNM
      DO 121 I=1,N
      S=P(I,J)
      IF (S.NE.0.D0) P(I,J)=S+T*HH(I,J)
121   CONTINUE
100   CONTINUE
      RETURN
      END
C
C
      SUBROUTINE MUSCLE (N,M,M1,X,DXT,DTAUT,DXTH,DTAUTH,XW,TAUW,XTHR,
     &   TAU,TAUMIN,TAUMAX,TQ)
C  PROVIDES SCALING XW(N,M) OF VARIABLES X(N,M)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION X(N,M),XW(N,M),DXT(N,M),DXTH(N,M),TQ(M)
      COMMON /MUVALS/ SMALL1,THMAX,THR,ITMAX,IGNMAX
C  THRESHOLD
      DO 222 I=1,N
      XMAX=.5D0*(DABS((TQ(2)-TQ(1))*X(I,1))+DABS((TQ(M)-TQ(M1))*X(I,M)))
      IF (M.EQ.2) GOTO 2231
      DO 223 J=2,M1
      XMAX=XMAX+DABS((TQ(J+1)-TQ(J-1))*X(I,J))
223   CONTINUE
2231  XMAX=XMAX/DABS(TQ(M)-TQ(1))
      IF (XMAX.LT.XTHR) XMAX=XTHR
      DO 224 J=1,M
      T=XMAX
      IF (T.LT.XW(I,J)) T=XW(I,J)
      DXTH(I,J)=DXT(I,J)*XW(I,J)/T
      XW(I,J)=T
224   CONTINUE
222   CONTINUE
      T=DABS(TAU)
      TSCAL=0.01D0*DMAX1(DABS(TAUMIN),DABS(TAUMAX))
      IF (T.LT.TSCAL) T=TSCAL
      IF (T.LT.TAUW) T=TAUW
      DTAUTH=DTAUT*TAUW/T
      TAUW=T
      RETURN
      END
C
C
      SUBROUTINE MURESI(FY,BC,IVPSOL,N,N1,M,M1,T,X,TAU,X1,XM,R,
     &   XU,HH,YTEMP,HSTART,TOL,INFO,UNRM,XW,TAUW,TFAIL,G,P,IVAR,W,DFY,
     &   HOS,KOS,INIT)
C  COMPUTATION OF TRAJECTORIES
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION T(M),X(N,M),X1(N),XM(N),R(N)
      DIMENSION XU(N,M1),HH(N,M1),YTEMP(N1),INFO(9)
      DIMENSION XW(N,M),G(N,N,M1),P(N,M1),W(N),HOS(M1),KOS(M1)
      INTEGER UPR,UDIAG
      EXTERNAL FY,DFY,BC,IVPSOL
      COMMON /COUNT/ IFCTEV,IDUM,IDECS,ISOLS,ITRAJ
      COMMON /MUUNIT/ UPR,UDIAG
      KPRINT=INFO(1)
      UNRM=0.D0
      NVAR=N**2+2*N+1
      ITRAJ=ITRAJ+1
2000  J=1
      HSAVE=HSTART
2100  J1=J+1
      TJ=T(J)
      TJ1=T(J1)
      H=HSAVE
      HMAX=DABS(TJ1-TJ)
      DO 2110 K=1,N
2110  YTEMP(K)=X(K,J)
      YTEMP(N1)=TAU
      CALL IVPSOL(N,NVAR,FY,DFY,TJ,YTEMP,TJ1,XW(1,J),XW(1,J1),TAUW,
     &   TOL,HMAX,H,UNRM,G(1,1,J),P(1,J),IVAR,W,HOS(J),KOS(J),INIT)
      HFIN=HSAVE
      HSAVE=H
      IF (H.EQ.0.D0) GOTO 9999
C  CONTINUITY CONDITIONS
2200  DO 2210 K=1,N
      TH=YTEMP(K)
      XU(K,J)=TH
2210  HH(K,J)=TH-X(K,J1)
      J=J1
      IF (J.LT.M) GOTO 2100
C  TWO-POINT BOUNDARY CONDITIONS
      DO 2300 I=1,N
      XM(I)=X(I,M)
2300  X1(I)=X(I,1)
      CALL BC(X1,XM,TAU,R)
      INIT=1
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  SINGULAR TRAJECTORY
9999  CONTINUE
      INFO(9)=J
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      RETURN
60001 FORMAT(' TERMINATION AT SINGULAR TRAJECTORY')
      END
C
C
      SUBROUTINE MUCHAB(N,N1,M,NE,NB,X,A,B,IA,IB,P,R,DE,XW,TAUW,
     &   IROW,ICOL,ICOLB,XTHR,INFO)
C  DETERMINATION OF SPARSE STRUCTURE OF MATRICES A AND B
C  AND DETERMINATION OF INTERNAL ROW SCALING OF SENSITIVITY MATRIX E
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION X(N,M),A(N,N),B(N,N),IA(N,N),IB(N,N),P(N),R(N),DE(N)
      DIMENSION XW(N,M),IROW(N),ICOL(N),ICOLB(N),INFO(9)
      INTEGER UPR,UDIAG
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /MUUNIT/ UPR,UDIAG
      KPRINT=INFO(1)
      ISUM=0
      DO 6100 I=1,N
      S=DABS(P(I))*TAUW
      DO 6110 K=1,N
      TH=DABS(A(I,K))*XW(K,1)
      IF (S.LT.TH) S=TH
      TH=DABS(B(I,K))*XW(K,M)
      IF (S.LT.TH) S=TH
      IF (IA(I,K).GT.0) GOTO 6111
      IF (A(I,K).EQ.0.D0) GOTO 6111
      IA(I,K)=1
      ISUM=1
6111  IF (IB(I,K).GT.0) GOTO 6110
      IF (B(I,K).EQ.0.D0) GOTO 6110
      IB(I,K)=1
      ISUM=1
6110  CONTINUE
      IF (S.LT.XTHR) S=XTHR
6100  DE(I)=1.D0/S
      IF (ISUM.EQ.0) GOTO 6400
C  DETERMINATION OF ROW AND COLUMN PERMUTATION VECTORS
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      DO 6210 I=1,N
      ICOL(I)=I
      ICOLB(I)=I
6210  IROW(I)=I
C  SEARCH FOR SEPARABLE LINEAR BOUNDARY CONDITIONS AT T(1)
      NE=N
      DO 6220 I=1,N
      IF (P(I).NE.0.D0) GOTO 6220
      DO 6221 K=1,N
      IF (IB(I,K).NE.0) GOTO 6220
6221  CONTINUE
      ISUM=0
      DO 6222 K=1,N
      IF (IA(I,K).EQ.0) GOTO 6222
      ISUM=ISUM+1
      ICA=K
6222  CONTINUE
      IF (ISUM.GT.1) GOTO 6220
      DO 6223 IS=1,N
      IH=ICOL(IS)
      IF (IH.EQ.ICA) ICOL(IS)=ICOL(NE)
      IH=IROW(IS)
      IF (IH.EQ.I) IROW(IS)=IROW(NE)
6223  CONTINUE
      ICOL(NE)=ICA
      IROW(NE)=I
      NE=NE-1
      IF (KPRINT.GT.0) WRITE(UPR,60002) I
      IF (DABS(R(I)).GT.10.D0*EPMACH*DABS(X(ICA,1))) GOTO 9999
6220  CONTINUE
C  SEARCH FOR SEPARABLE LINEAR BOUNDARY CONDITIONS AT T(M)
      NB=0
      IF (NE.EQ.0) GOTO 6400
      DO 6230 I=1,NE
      IR=IROW(I)
      IF (P(IR).NE.0.D0) GOTO 6230
      DO 6231 K=1,N
      IF (IA(IR,K).NE.0) GOTO 6230
6231  CONTINUE
      ISUM=0
      DO 6232 K=1,N
      IF (IB(IR,K).EQ.0) GOTO 6232
      ISUM=ISUM+1
      ICB=K
6232  CONTINUE
      IF (ISUM.GT.1) GOTO 6230
      NB=NB+1
      DO 6233 IS=1,N
      IH=ICOLB(IS)
      IF (IH.EQ.ICB) ICOLB(IS)=ICOLB(NB)
6233  CONTINUE
      ICOLB(NB)=ICB
      IROW(I)=IROW(NB)
      IROW(NB)=IR
      IF (KPRINT.GT.0) WRITE(UPR,60003) IR
      IF (DABS(R(IR)).GT.10.D0*EPMACH*DABS(X(ICB,M))) GOTO 9999
6230  CONTINUE
      IF (KPRINT.GT.0 .AND. NE.EQ.N .AND. NB.EQ.0) WRITE(UPR,60005)
      IF (KPRINT.GT.0) WRITE(UPR,60004)
6400  CONTINUE
      INFO(9)=0
      RETURN
C  FAIL EXIT
C  NON-ZERO RESIDUAL IN SEPARATED BOUNDARY CONDITION
9999  CONTINUE
      INFO(9)=1
      RETURN
60001 FORMAT(/' CHECK FOR SEPARATED BOUNDARY CONDITIONS INITIATED'/)
60002 FORMAT(' BOUNDARY CONDITION',I3,
     &   ' ASSUMED TO BE SEPARATED AT T(1)')
60003 FORMAT(' BOUNDARY CONDITION',I3,
     &   ' ASSUMED TO BE SEPARATED AT T(M)')
60004 FORMAT(' ')
60005 FORMAT(' NO BOUNDARY CONDITIONS ASSUMED TO BE SEPARATED')
      END
C
C
      SUBROUTINE MUCDNS(N,N1,M,M1,NE,NE1,G,A,B,XW,DE,E,ES,W,IROW,ICOL,
     &   P,TAUW,BG,V)
C  COMPUTATION OF SENSITIVITY MATRIX E=-(A+B*G(M1)*...*G(1))
C  (PROJECTIONS INCLUDED)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DIMENSION G(N,N,M1),A(N,N),B(N,N),XW(N,M),DE(N),E(N,N1),ES(N,N1)
      DIMENSION W(N),IROW(N),ICOL(N)
      DIMENSION P(N,M),BG(N,N),V(N)
      DO 7100 I=1,NE
      IR=IROW(I)
      DO 7100 K=1,N
7100  E(I,K)=B(IR,K)*DE(IR)
      DO 7200 JJ=1,M1
      J=M-JJ
      DO 7200 I=1,NE
      DO 7210 K=1,N
      S=0.D0
      DO 7211 L=1,N
7211  S=S+E(I,L)*G(L,K,J)
7210  W(K)=S
      DO 7220 K=1,N
7220  E(I,K)=W(K)
7200  CONTINUE
C  INTERNAL ROW AND COLUMN SCALING AND PERMUTATION OF MATRIX E
      DO 7400 K=1,NE
      KC=ICOL(K)
      S=XW(KC,1)
      DO 7400 I=1,NE
      IR=IROW(I)
7400  E(I,K)=-(A(IR,KC)*DE(IR)+E(I,KC))*S
C  CONDENSED PARAMETER COLUMN
      CALL MURHS1(N,NE,M1,1,P,P(1,M),B,G,V,DE,W,BG,IROW)
      DO 7401 I=1,NE
7401  E(I,NE1)=-V(I)*TAUW
      RETURN
      END
C
C
      SUBROUTINE MUREFN (N,N1,NE,NE1,ES,DX1,DU,US,IRANKB,W,IRANK,
     &   D,PIVOT,WM,EPH,NQ,EH,EPS,E,LEVEL,INFO)
C  ITERATIVE REFINEMENT OF DX1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION ES(N,N1),DX1(N1),DU(N),US(N),W(N1),D(N1),WM(N1),
     $EH(N1,N1),E(N,N1)
      INTEGER PIVOT(N1),INFO(9)
      INTEGER UPR,UDIAG
      COMMON /MUUNIT/ UPR,UDIAG
      EPH=EPS
      DO 3140 I=1,NE
      S=0.D0
      DO 3141 K=1,NQ
3141  S=S+ES(I,K)*DX1(K)
3140  DU(I)=US(I)-S
      DO 3142 I=1,NQ
3142  WM(I)=0.D0
C  SOLUTION OF RESIDUAL EQUATION
      IF (LEVEL.EQ.2) GOTO 4000
      CALL MUSOLC(E,N,N1,IRANKB,NE,NQ,WM,DU,IRANK,D,PIVOT,0,EH,W)
      GOTO 4009
4000  CONTINUE
      MH = IRANKB
      IF (MH.EQ.0)  MH = NE
      DO  31  J=1,IRANK
      S = 0.D0
      DO  311  I=J,MH
311   S = S+E(I,J)*DU(I)
      S = S/(D(J)*E(J,J))
      DO  312  I=J,NE
312   DU(I) = DU(I)+E(I,J)*S
      IF (J.EQ.IRANKB)  MH = NE
 31   CONTINUE
      IRK1 = IRANK+1
      DO  41  II=1,IRANK
      I = IRK1-II
      I1 = I + 1
      S = DU(I)
      IF (I1.GT.IRANK)  GO TO 41
      DO  4111  JJ=I1,IRANK
4111  S = S-E(I,JJ)*W(JJ)
41    W(I) = S/D(I)
      DO  50  J=1,NE
      IH=PIVOT(J)
50    WM(IH) = W(J)
4009  CONTINUE
      EPX1=0.D0
      EPDX1=0.D0
      DO 3150 I=1,NQ
      S=DABS(WM(I))
      IF(EPDX1.LT.S) EPDX1=S
      DX1(I)=DX1(I)+WM(I)
      S=DABS(DX1(I))
      IF(EPX1.LT.S) EPX1=S
3150  CONTINUE
      EPX1H=0.D0
      IF (EPX1.NE.0.D0) EPX1H=EPDX1/EPX1
      EPH=10.D0*EPDX1
      IF(EPX1H.GT.EPS) GOTO 9800
3160  INFO(9)=0
      RETURN
9800  INFO(9)=1
      IF (INFO(1).GT.0) WRITE(UPR,60001)
      RETURN
60001 FORMAT(' GAUSSIAN BLOCK ELIMINATION FAILED BY ILL-CONDITIONED',
     &   ' CONDENSED LINEAR SYSTEM')
      END
C
C
      SUBROUTINE MUSWP (N,N1,M,M1,NY,NYMAX,EPS,EPH,HH,G,P,DXQ,DHH,
     &   R,A,B,DR,DU,DE,W,WM,IROW,DX1,NE,NQ,IRANK,IRANKB,PIVOT,D,E,EH,
     &   ICOL,XW,TAUW,DDX,RF,DTAU,LEVEL,RELDIF,TOLF,TOLJ,INFO)
C  ITERATIVE REFINEMENT SWEEPS  NY=1,..,NYMAX
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION HH(N,M1),G(N,N,M1),P(N,M),DXQ(N,M),DHH(N,M1),
     $R(N),A(N,N),B(N,N),DR(N),DU(N),DE(N),W(N1),WM(N,N),
     $DX1(N1),D(N1),E(N,N1),EH(N1,N1),XW(N,M),DDX(N,M),RF(M)
      INTEGER PIVOT(N1),IROW(N),ICOL(N),INFO(9)
      INTEGER UPR,UDIAG
      COMMON /MUUNIT/ UPR,UDIAG
      DATA TOLMIN/1.D-12/,REDH/1.D-2/
      KPRINT=INFO(1)
      NY=0
      SIGDEL=10.D0
      SIGDLH=0.D0
      IF (EPH.LT.EPS) EPH=EPS
C
C  COMPUTATION OF REQUIRED CONTINUITY RESIDUALS DHH(N,M1)
      JN=1
      JIN=M
      GOTO 3230
C
3200  DO 3220 J=JN,M1
      J1=J+1
      DO 3220 I=1,N
      S=HH(I,J)+P(I,J)*DTAU
      DO 3221 K=1,N
3221  S=S+G(I,K,J)*DXQ(K,J)
3220  DHH(I,J)=S-DXQ(I,J1)
C
C  COMPUTATION OF BOUNDARY RESIDUAL DR(N)
3230  DO 3240 I=1,N
      S=R(I)+P(I,M)*DTAU
      DO 3241 K=1,N
3241  S=S+A(I,K)*DXQ(K,1)+B(I,K)*DXQ(K,M)
3240  DR(I)=S
C
C  COMPUTATION OF CONDENSED RESIDUAL DU(NE)
      CALL MURHS1(N,NE,M1,JIN,DHH,DR,B,G,DU,DE,W,WM,IROW)
C
C  COMPUTATION OF CORRECTION DDX(N,1)
      DO 3250 I=1,N1
3250  DX1(I)=0.D0
      IF (LEVEL.GE.2) GOTO 40000
      CALL MUSOLC
     &   (E,N,N1,IRANKB,NE,NQ,DX1,DU,IRANK,D,PIVOT,0,EH,W)
      GOTO 40009
40000 CONTINUE
      MH = IRANKB
      IF (MH.EQ.0)  MH = NE
      DO  31  J=1,IRANK
      S = 0.D0
      DO  311  I=J,MH
311   S = S+E(I,J)*DU(I)
      S = S/(D(J)*E(J,J))
      DO  312  I=J,NE
312   DU(I) = DU(I)+E(I,J)*S
      IF (J.EQ.IRANKB)  MH = NE
 31   CONTINUE
      IRK1 = IRANK+1
      DO  41  II=1,IRANK
      I = IRK1-II
      I1 = I + 1
      S = DU(I)
      IF (I1.GT.IRANK)  GO TO 41
      DO  4111  JJ=I1,IRANK
4111  S = S-E(I,JJ)*W(JJ)
41    W(I) = S/D(I)
      DO  50  J=1,NE
      IH=PIVOT(J)
50    DX1(IH) = W(J)
40009 CONTINUE
      DDTAU=0.D0
      IF (NQ.GT.NE) DDTAU=DX1(NQ)
      IF (NQ.GT.NE) DX1(NQ)=0.D0
C
C  DESCALING OF DDX(N,1), REFINEMENT OF DXQ(N,1)
C
      CORR=0.D0
      DO 3260 L=1,N
      I=ICOL(L)
      S=DX1(L)
3260  DDX(I,1)=S
      DO 3266 I=1,N
      S=DDX(I,1)
      IF (CORR.LT.DABS(S)) CORR=DABS(S)
      S=S*XW(I,1)
      DDX(I,1)=S
3266  DXQ(I,1)=DXQ(I,1)+S
      S=DDTAU
      IF (CORR.LT.DABS(S)) CORR=DABS(S)
      S=S*TAUW
      DDTAU=S
      DTAU=DTAU+S
      IF (CORR.LT.EPH) GOTO 3269
      EPH=CORR
      GOTO 9800
3269  RF(1)=CORR
C
C  RECURSIVE COMPUTATION OF DDX(N,2),...,DDX(N,M)
      CALL MURECU(N,M,M1,JIN,DHH,G,P,DDX,DDTAU,W,WM)
C
C  REFINEMENT OF DXQ(N,2),...,DXQ(N,M)
      DO 3270 J=2,M
      CORR=0.D0
      DO 3271 I=1,N
      S=DDX(I,J)
      DXQ(I,J)=DXQ(I,J)+S
      S=DABS(S)/XW(I,J)
      IF(CORR.LT.S) CORR=S
3271  CONTINUE
3270  RF(J)=CORR
C
C  DETERMINATION OF MUSWP INDEX JN
      JA=JN
      DO 3280 J=1,M
      IF(RF(J).GT.EPH) GOTO 3290
3280  JN=J
C
3290  NY=NY+1
      IF(JN.LE.JA) GOTO 9600
      IF(JN.EQ.M) GOTO 3900
      JIN=JN
      IF(NY.GT.1 .OR. LEVEL.LE.2) GOTO 3200
C
C  DETERMINATION AND ADAPTATION OF PARAMETERS TOL AND RELDIF
3900  IF(LEVEL.LE.2 .OR. NY.GT.1) GOTO 3920
      DO 3910 J=1,M1
      S=0.D0
      IF (RF(J).NE.0.D0) S=RF(J+1)/RF(J)
      IF(SIGDLH.LT.S) SIGDLH=S
      RF(J)=S
3910  CONTINUE
      SIGDEL=DMAX1(SIGDLH,SIGDEL)
      TH=TOLF*SIGDEL
      IF(TH.GT.REDH) GOTO 9700
      IF(TH.GT.EPH) EPH=TH
      TOLH=EPS/SIGDEL
      IF(TOLH.LT.TOLMIN) TOLH=TOLMIN
      TOLF=TOLH
      TOLJ=DSQRT(TOLF)
      IF (TOLJ.GT.1.D-3) TOLJ=1.D-3
      RELDIF=DSQRT(TOLF/SIGDEL)
3920  IF(JN.NE.M) GOTO 3200
      INFO(9)=0
      RETURN
C  FAIL EXIT
9600  INFO(9)=1
      IF (KPRINT.GT.0) WRITE(UPR,60001)
      RETURN
9700  INFO(9)=2
      IF (KPRINT.GT.0) WRITE(UPR,60002)
      RETURN
9800  INFO(9)=3
      IF (KPRINT.GT.0) WRITE(UPR,60003)
      RETURN
60001 FORMAT(' ITERATIVE REFINEMENT FAILED TO CONVERGE')
60002 FORMAT(' RELIABLE RELATIVE ACCURACY NOT SUFFICIENT')
60003 FORMAT(' GAUSSIAN BLOCK ELIMINATION FAILED BY ILL-CONDITIONED',
     &   ' CONDENSED LINEAR SYSTEM')
      END
C
C
      FUNCTION MUSCPR(N,M,X1,X2,TAU1,TAU2,XW,TAUW,T)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION X1(N,M),X2(N,M),XW(N,M),T(M)
      DOUBLE PRECISION MUSCPR
      RL=T(M)-T(1)
      MUSCPR=0.D0
      DO 1000 J=1,M
      IF (J.EQ.1) S=T(2)-T(1)
      IF (J.NE.1 .AND. J.NE.M) S=T(J+1)-T(J-1)
      IF (J.EQ.M) S=T(M)-T(M-1)
      S=S/RL
      SUM=0.D0
      DO 1010 I=1,N
      S1=X1(I,J)/XW(I,J)
      S2=X2(I,J)/XW(I,J)
1010  SUM=SUM+S1*S2
1000  MUSCPR=MUSCPR+S*SUM
      MUSCPR=0.5D0*MUSCPR
      MUSCPR=MUSCPR+(TAU1/TAUW)*(TAU2/TAUW)
      RETURN
      END
C
C
      SUBROUTINE MUDIFX (NQ,NVAR,FY,DFY,X,YQ,XEND,XW1,XW2,TAUW,
     &   EPS,HMAX,H,UNRM,GQ,PQ,IVAR,WQ,HOS,KOS,INIT)
C
C
C  EXPLICIT EXTRAPOLATION INTEGRATOR
C  FOR NON-STIFF SYSTEMS OF ORDINARY DIFFERENTIAL EQUATIONS
C  (BASED ON THE EXPLICIT MID-POINT DISCRETIZATION)
C
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      INTEGER NJ(10),INCR(10),NRED( 9)
C     INTEGER NJ(JM),INCR(JM),NRED(KM)
      DOUBLE PRECISION YQ(1),Y(99),YL(99),YM(99),DY(99),DZ(99),
     $  DT(99,10),
C                      Y(N),YL(N) ,YM(N) ,DY(N) ,DZ( N),DT( N,JM),
     $  D(10,10),A(10),AL(10,10),GQ(NQ,NQ),PQ(NQ),XW1(NQ),XW2(NQ),WQ(NQ)
C       D(JM,JM),A(JM),AL(JM,JM)
      DOUBLE PRECISION B,B1,C,DABS,DBLE,DM,DMAX,DSQRT,EPH,EPMACH,
     1 EPS,ERR,FC,FCM,FCO,FIVE,FJ,FJ1,FMIN,FN,G,H,HALF,HMAX,HMAXU,HR,
     2 H1,OMJ,OMJO,ONE,ONE1,Q,QUART,RED,RO,SAFE,SMALL,TA,TEN,
     3 TWO,U,V,W,X,XEND,XEPS,XN,YH,ZERO,TAU,UNRM
      DOUBLE PRECISION MURNTG
      COMMON /MACHIN/ EPMACH,SMALL
      COMMON /COUNT/ IFCTEV,ITER,IDECS,ISOLS,ITRAJ
      EXTERNAL FY,DFY
      DATA ZERO/0.D0/,FMIN/1.D-2/,RO/0.25D0/,QUART/0.25D0/,HALF/0.5D0/,
     *     SAFE/0.5D0/,ONE/1.D0/,ONE1/1.01D0/,TWO/2.D0/,TEN/1.D1/,
     *     FIVE/5.D0/
      DATA DT/990*0.D0/
C     DATA DT/N*JM*0.D0/
C
C  STEPSIZE SEQUENCE HA (DUE TO /3/ )
      DATA NJ/2,4,6,8,10,12,14,16,18,20/
C
C  ASSOCIATED MAXIMUM COLUMN NUMBER (1.LE.KM.LE.9)
      KM=6
C
C  ASSOCIATED MAXIMUM ROW NUMBER (2.LE.JM.LE.10)
      JM=KM+1
C
C
C  INTERNAL PARAMETERS
C  STANDARD VALUES FIXED BELOW
C
C
C  MAXIMUM PERMITTED NUMBER OF INTEGRATION STEPS PER INTERVAL
      NSTMAX=10000
C
C  MAXIMUM PERMITTED NUMBER OF STEPSIZE REDUCTIONS
      JRMAX=5
C
C  INITIAL PREPARATIONS
      N=NQ+1
      IF (IVAR.EQ.1) N=NVAR
      DO 49 I=1,NQ
49    Y(I)=YQ(I)
      TAU=YQ(NQ+1)
      Y(NQ+1)=UNRM**2
      IF (IVAR.EQ.0) GOTO 50
      DO 51 I=1,NQ
      DO 52 J=1,NQ
52    Y(2*NQ+1+(J-1)*NQ+I)=0.D0
      Y(NQ+1+I)=0.D0
51    Y(2*NQ+1+(I-1)*NQ+I)=1.D0
50    CONTINUE
      NSTEP=0
      EPH=RO*EPS
      FJ1=DBLE(NJ(1))
      A(1)=FJ1+ONE
      DO 60 J=2,JM
      J1=J-1
      INCR(J1)=0
      NRED(J1)=0
      FJ=DBLE(NJ(J))
      V=A(J1)+FJ
      A(J)=V
      DO 61 K=1,J1
      W=FJ/DBLE(NJ(K))
61    D(J,K)=W*W
      IF(J.EQ.2) GOTO 60
      W=V-FJ1
      DO 62 K1=2,J1
      K=K1-1
      U=(A(K1)-V)/(W*DBLE(K+K1))
      U=EPH**U
62    AL(J1,K)=U
60    CONTINUE
      KOH=1
      JOH=2
65    IF(JOH.GE.JM) GOTO 66
      IF(A(JOH+1)*ONE1.GT.A(JOH)*AL(JOH,KOH)) GOTO 66
      KOH=JOH
      JOH=JOH+1
      GOTO 65
66    K=0
      KM=KOH
      JM=KM+1
      INCR(JM)=-1
      OMJO=ZERO
      HMAX=DABS(HMAX)
      XEPS=(DABS(X)+DABS(XEND))*EPMACH*TEN
      FN=DBLE(N)
      H1=XEND-X
      HMAXU=HMAX
      HR=HMAX
      DMAX=FIVE
      IF (INIT.EQ.0) GOTO 401
      KOH=KOS
      JOH=KOH+1
      H=HOS
C
C  BASIC INTEGRATION STEP
401   IF(DABS(H1).LE.XEPS) GOTO 403
      Q=H1/H
      IF(Q.LE.EPMACH) GOTO 403
      IF(Q.GE.ONE1) GOTO 402
      HR=H
      H=H1
402   JRED=0
      DO 405 K=1,KM
405   INCR(K)=INCR(K)+1
      HMAX=DABS(H1)
      IF(HMAXU.LT.HMAX) HMAX=HMAXU
C
C  EXPLICIT EULER STARTING STEP
      CALL MUFVAR(NQ,FY,DFY,X,Y,Y(NQ+1),Y(NQ+2),Y(2*NQ+2),TAU,
     &   DZ,DZ(NQ+1),DZ(NQ+2),DZ(2*NQ+2),IVAR,WQ)
      IFCTEV=IFCTEV+1
C
10    XN=X+H
      FCM=DABS(H)/HMAX
      IF(FCM.LT.FMIN) FCM=FMIN
C
      DO 260 J=1,JM
      M=NJ(J)
      G=H/DBLE(M)
      B=G+G
      DO 210 I=1,N
      YL(I)=Y(I)
210   YM(I)=Y(I)+G*DZ(I)
      M=M-1
C  EXPLICIT MID-POINT RULE
      DO 220 K=1,M
      CALL MUFVAR(NQ,FY,DFY,X+G*DBLE(K),YM,YM(NQ+1),YM(NQ+2),
     &   YM(2*NQ+2),TAU,DY,DY(NQ+1),DY(NQ+2),DY(2*NQ+2),IVAR,WQ)
      IFCTEV=IFCTEV+1
      DO 220 I=1,N
      U=YL(I)+B*DY(I)
      YL(I)=YM(I)
      YM(I)=U
220   CONTINUE
C  FINAL STEP
      CALL MUFVAR(NQ,FY,DFY,XN,YM,YM(NQ+1),YM(NQ+2),YM(2*NQ+2),TAU,
     &   DY,DY(NQ+1),DY(NQ+2),DY(2*NQ+2),IVAR,WQ)
      IFCTEV=IFCTEV+1
      DO 2200 I=1,N
      YH=YL(I)+G*DY(I)
      DY(I)=YH-YM(I)
      YM(I)=(YM(I)+YH)*HALF
2200  CONTINUE
      DM=MURNTG(NQ,DY,DY(NQ+1),DY(NQ+2),DY(2*NQ+2),XW1,XW2,TAUW,IVAR)
C
C  STABILITY CHECK
      IF(DM.LT.DMAX) GOTO 2209
C
C  EMERGENCY EXIT
      GOTO 2601
C
C  PREVENTION OF POSSIBLE ORDER INCREASE
2209  IF(J.GT.2.OR.DM.LT.DMAX*HALF) GOTO 2207
      DO 2208 L=JOH,JM
      IF(INCR(L).GT.0) INCR(L)=0
      INCR(L)=INCR(L)-2
2208  CONTINUE
C
C  EXTRAPOLATION
2207  CONTINUE
      DO 234 I=1,N
      V=DT(I,1)
      C=YM(I)
      DT(I,1)=C
      IF(J.EQ.1) GOTO 234
      TA=C
      DO 231 K=2,J
      JK=J-K+1
      B1=D(J,JK)
      W=C-V
      U=W/(B1-ONE)
      C=B1*U
      V=DT(I,K)
      DT(I,K)=U
231   TA=U+TA
      YM(I)=TA
      DY(I)=U
234   CONTINUE
      ERR=0.D0
      IF(J.EQ.1) GOTO 260
      ERR=MURNTG(NQ,DY,DY(NQ+1),DY(NQ+2),DY(2*NQ+2),
     &   XW1,XW2,TAUW,IVAR)
C ERROR (SCALED ROOT MEAN SQUARE)
      ERR=ERR/DSQRT(FN)
      KONV=0
      IF(ERR.LT.EPS) KONV=1
      ERR=ERR/EPH
C
C ORDER CONTROL
      K=J-1
      L=J+K
      FC=ERR**(ONE/DBLE(L))
      IF(FC.LT.FCM) FC=FCM
C  OPTIMAL ORDER DETERMINATION
      OMJ=FC*A(J)
      IF(J.GT.2.AND.OMJ*ONE1.GT.OMJO.OR.K.GT.JOH) GOTO 235
      KO=K
      JO=J
      OMJO=OMJ
      FCO=FC
235   IF(J.LT.KOH.AND.INIT.GT.0) GOTO 260
      IF(KONV.EQ.0) GOTO 236
      IF(KO.LT.K.OR.INCR(J).LT.0) GOTO 20
C  POSSIBLE INCREASE OF ORDER
      IF(NRED(KO).GT.0) NRED(KO)=NRED(KO)-1
      FC=FCO/AL(J,K)
      IF(FC.LT.FCM) FC=FCM
      J1=J+1
      IF(A(J1)*FC*ONE1.GT.OMJO) GOTO 20
      FCO=FC
      KO=JO
      JO=JO+1
      GOTO 20
C
C
C  CONVERGENCE MONITOR
236   RED=ONE/FCO
      JK=KM
      IF(JOH.LT.KM) JK=JOH
      IF(K.GE.JK) GOTO 239
      IF(KO.LT.KOH) RED=AL(KOH,KO)/FCO
237   IF(AL(JK,KO).LT.FCO) GOTO 239
260   CONTINUE
C
C STEPSIZE REDUCTION (DUE TO EXTRAPOLATION TABLE)
239   RED=RED*SAFE
      H=H*RED
2392  IF(NSTEP.EQ.0) GOTO 2390
      NRED(KOH)=NRED(KOH)+1
      DO 2391 L=KOH,KM
2391  INCR(L)=-2-NRED(KOH)
2390  JRED=JRED+1
      IF(JRED.GT.JRMAX) GOTO 32
      GOTO 10
C
C  STEPSIZE REDUCTION (DUE TO STABILITY)
2601  HMAX=G*FJ1*QUART
      RED=HMAX/DABS(H)
      H=HMAX
      IF(JRED.GT.0) GOTO 2390
      GOTO 2392
C
C  PREPARATIONS FOR NEXT BASIC INTEGRATION STEP
20    X=XN
      H1=XEND-X
      DO  2606 I=1,N
2606  Y(I)=YM(I)
      NSTEP=NSTEP+1
      IF(NSTEP.GT.NSTMAX) GO TO 31
C
C STEPSIZE PREDICTION
      IF(FCO.NE.FCM) HR=H
      H=H/FCO
      KOH=KO
      JOH=KOH+1
      IF (NSTEP.GT.1) GOTO 270
      KOS=KOH
      HOS=H
270   CONTINUE
      IF(X+H.NE.X) GO TO 401
      GO TO 33
C
C  SOLUTION EXIT
403   H=HR
      HMAX=HMAXU
      DO 449 I=1,NQ
449   YQ(I)=Y(I)
      UNRM=DSQRT(Y(NQ+1))
      YQ(NQ+1)=TAU
      IF (IVAR.EQ.0) RETURN
      DO 451 I=1,NQ
      DO 452 J=1,NQ
452   GQ(I,J)=Y(2*NQ+1+(J-1)*NQ+I)
451   PQ(I)=Y(NQ+1+I)
      RETURN
C
C  FAIL EXIT
31    CONTINUE
      GOTO 39
32    CONTINUE
      GOTO 39
33    CONTINUE
39    H=ZERO
      HMAX=HMAXU
      RETURN
C
C  END MUDIFX
C
      END
C
C
      SUBROUTINE MUFVAR(N,FY,DFY,T,Y,U,P,G,TAU,DY,DU,DP,DG,IVAR,W)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION Y(N),P(N),G(N,N),DY(N),DP(N),DG(N,N),W(N)
      EXTERNAL FY,DFY
      CALL FY(T,Y,TAU,DY)
      DU=0.D0
      DO 10 I=1,N
10    DU=DU+Y(I)**2
      IF (IVAR.EQ.0) RETURN
      CALL DFY(T,Y,TAU,DG,DP)
      DO 20 I=1,N
      S=DP(I)
      DO 21 J=1,N
21    S=S+DG(I,J)*P(J)
20    DP(I)=S
      DO 30 I=1,N
      DO 31 J=1,N
      S=0.D0
      DO 32 K=1,N
32    S=S+DG(I,K)*G(K,J)
31    W(J)=S
      DO 33 J=1,N
33    DG(I,J)=W(J)
30    CONTINUE
      RETURN
      END
C
C
      FUNCTION MURNTG(N,X,U,P,G,XW1,XW2,TAUW,IVAR)
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)
      DOUBLE PRECISION X(N),P(N),G(N,N),XW1(N),XW2(N)
      DOUBLE PRECISION MURNTG
      MURNTG=0.D0
      DO 10 I=1,N
10    MURNTG=MURNTG+(X(I)/XW1(I))**2
      IF (IVAR.EQ.0) GOTO 90
      DO 20 I=1,N
20    MURNTG=MURNTG+(P(I)*TAUW/XW1(I))**2
      DO 30 I=1,N
      DO 30 J=1,N
30    MURNTG=MURNTG+(G(I,J)*XW1(J)/XW2(I))**2
90    MURNTG=DSQRT(MURNTG)
      RETURN
      END
