C     ALGORITHM 617, COLLECTED ALGORITHMS FROM ACM.
C     ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.10, NO. 3,
C     SEP., 1984, P. 317-324.
      SUBROUTINE VECF (N,X,VFU,NR,IER)
C
C THE SUBROUTINE  VECF  RETURNS IN THE N-VECTOR  VFU
C THE VALUES OF THE FUNCTIONS OF SYSTEM (1) (SEE SUBROUTINE  EVOLV )
C COMPUTED AT THE POINT  X .
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C INDICATES A FAILURE IN COMPUTING A FUNCTION VALUE).
C
C SAMPLE VERSION
C THE PRESENT VERSION OF THE SUBROUTINE  VECF  COMPUTES
C THE FUNCTION VALUES FOR THE SELECTED SAMPLE PROBLEM
C IN THE PRESENT VERSION THE PARAMETER  NR  IS NEVER USED.
C
      DOUBLE PRECISION VFU,X
      DOUBLE PRECISION BIG
      DIMENSION X(N),VFU(N)
      DATA BIG/1.D35/
C
C TEST RANGE OF VARIABLES.
      IF (DMAX1(DABS(X(1)),DABS(X(2))).GT.BIG) GO TO 10
C
C     ROSENBROCK S FUNCTION.
C
      VFU(1) = 1.D0-X(1)
      VFU(2) = 10.D0*(X(2)-X(1)**2)
C
C NORMAL EXIT
      IER = 0
      RETURN
C
C ERROR EXIT
  10  CONTINUE
      IER = -1
      RETURN
C
      END
      SUBROUTINE JCBNA (N,X,C,NR,IDA,IER)
C
C THE SUBROUTINE  JCBNA  COMPUTES THE ANALYTICAL JACOBIAN (IF AVAILABLE)
C THE PARAMETERS ARE DESCRIBED WITHIN THE SUBROUTINE  JCBN
C
C SAMPLE VERSION
C THE PRESENT VERSION COMPUTES THE ANALYTICAL JACOBIAN
C FOR THE SELECTED TEST PROBLEM.
C IN THE PRESENT VERSION THE PARAMETER  NR  IS NEVER USED.
C
      DOUBLE PRECISION C,X
      DIMENSION X(N),C(N,N)
C
      IER = 0
      IDA = 0
C
C ROSENBROCK FUNCTION
      IDA = 1
      C(1,1) = -1.D0
      C(1,2) = 0.D0
      C(2,1) = -20.D0*X(1)
      C(2,2) = 10.D0
      RETURN
      END
      SUBROUTINE DAFNE (N,X,FMIN,PRF,IDER,NU,NO,IDA,
     1    NFEV,NJEV,NSAC,NST,IER,
     2    V,FP,VFU,BUF,XTP,WF,FPS,S,C,CS)
C
C THE SUBROUTINE  DAFNE  ATTEMPTS TO SOLVE A SYSTEM
C OF  N  NONLINEAR EQUATIONS
C     FI(X) = 0,  I = 1,...,N,  X = (X1,...,XN)        (1)
C BY MEANS OF THE DIFFERENTIAL-EQUATION APPROACH DESCRIBED
C WITHIN THE SUBROUTINE  EVOLV .
C THE NUMERICAL INTEGRATION AND RESCALING CAPABILITY OF  EVOLV
C IS EXPLOITED BY  DAFNE  TO OBTAIN A TRAJECTORY WHICH IS AT
C TIMES STOPPED TO RESCALE THE FUNCTION AND THEN RESTARTED FROM
C THE POINT REACHED.
C THIS IS IMPLEMENTED BY PERFORMING
C A SEQUENCE OF CALLS TO THE SUBROUTINE  EVOLV ,
C SETTING SUITABLE VALUES FOR THE CALL PARAMETERS,
C EACH TIME USING AS INITIAL POINT THE FINAL POINT
C RETURNED BY THE PRECEDING CALL, AND COMPUTING
C A NEW VECTOR OF WEIGHTS
C     WF = (WF1,WF2,...,WFN)
C USED FOR FUNCTION RESCALING.
C THE RESCALING WEIGHTS ARE COMPUTED BY CALLING
C THE SUBROUTINE  NEWF .
C
C CALL STATEMENT :
C     CALL DAFNE (N,X,FMIN,PRF,IDER,NU,NO,IDA,
C    1    NFEV,NJEV,NSAC,NST,IER,
C    2    B1,B2,B3,B4,B5,B6,B7,B8,C1,C2)
C
C DESCRIPTION OF PARAMETERS OF THE CALL STATEMENT
C     (FORTRAN IMPLICIT TYPE DEFINITION IS USED.
C     ALL NON-INTEGER PARAMETERS ARE DOUBLE-PRECISION)
C
C N  IS THE (INPUT) NUMBER OF NONLINEAR EQUATIONS.
C
C X  IS A VECTOR OF LENGTH  N  CONTAINING THE X-VARIABLES . (INITIAL
C     VALUES ON INPUT, FINAL ON OUTPUT. INITIAL AND FINAL  X  ARE
C     THE INITIAL AND FINAL ESTIMATES OF THE SOLUTION OF SYSTEM (1)).
C
C FMIN  IS  AN INPUT VECTOR OF DIMENSION  N  CONTAINING THE THRESHOLD
C     VALUES FOR THE FUNCTION VALUES  FI(X) , I = 1,...,N .
C     CONVERGENCE IS CLAIMED IF
C            ABS(FI(X)).LT.FMIN(I) , I = 1,...,N
C
C THE FOLLOWING 5 PARAMETERS (PRF,IDER,NU,NO,IDA)
C ARE USED ONLY IN THE COMPUTATION OF THE JACOBIAN MATRIX
C ( SUBROUTINE  JCBN ).
C
C PRF  IS AN (INPUT) ROUGH ESTIMATE OF THE AVERAGE RELATIVE
C     ERROR AFFECTING THE FUNCTION VALUES SUPPLIED BY SUBROUTINE  VECF .
C     IF THIS ERROR IS ESTIMATED TO BE NOT SEVERAL ORDERS OF MAGNITUDE
C     GREATER THAN MACHINE PRECISION, THE USER MAY SIMPLY PUT  PRF = 0 .
C
C IDER  IS AN INPUT CONTROL VARIABLE TO REQUEST THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDER.EQ.0   ANALYTICAL COMPUTATION BY SUBROUTINE  JCBNA
C     IDER.EQ.1   COMPUTATION BY FORWARD FINITE DIFFERENCES
C     IDER.EQ.2   COMPUTATION BY CENTRAL FINITE DIFFERENCES
C
C NU  AND  NO  ARE NUMBERS DEFINING THE (POSSIBLY) BANDED
C STRUCTURE OF THE JACOBIAN MATRIX, USED FOR FINITE-DIFFERENCES
C JACOBIAN COMPUTATION
C     NU IS THE N. OF NON-ZERO JACOBIAN DIAGONALS UNDER MAIN DIAGONAL
C     NO IS THE N. OF NON-ZERO JACOBIAN DIAGONALS ABOVE MAIN DIAGONAL
C     (FOR A NON-BANDED MATRIX  NU = NO = N-1 )
C
C IDA  IS AN (OUTPUT) INDICATOR OF THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDA.EQ.1   THE DERIVATIVES WERE COMPUTED ANALYTICALLY
C                (AS REQUESTED)
C     IDA.EQ.0   THE DERIVATIVES WERE COMPUTED BY FINITE DIFFERENCES.
C                (CENTRAL FINITE DIFFERENCES IF REQUESTED BY  IDER = 2,
C                FORWARD FINITE DIFFERENCES IF REQUESTED BY  IDER = 1,
C                OR IF REQUESTED ANALYTICAL DERIVATIVES WERE NOT CODED
C                IN SUBROUTINE  JCBNA ).
C
C THE FOLLOWING 4 PARAMETERS (NFEV,NJEV,NSAC,NST) ARE OUTPUT COUNTERS
C     NFEV  OF FUNCTION EVALUATIONS (EXCLUDING THOSE USED
C           FOR JACOBIAN EVALUATION).
C     NJEV  OF JACOBIAN EVALUATIONS.
C     NSAC  OF ACCEPTED TIME INTEGRATION STEPS.
C     NST   OF ATTEMPTED (ACCEPTED + REJECTED) TIME INTEGRATION STEPS.
C
C IER  IS AN OUTPUT ERROR INDICATOR AS FOLLOWS
C     IER = 0    SUCCESS
C     IER = 1    FAILURE (MAX. NUMBER OF ATTEMPTS)
C     IER = 2    FATAL FAILURE IN COMPUTING THE FUNCTION VALUES
C                        AT THE INITIAL POINT OR THE DERIVATIVES.
C
C THE LAST 10 PARAMETERS OF THE CALL STATEMENT
C ARE WORKING STORAGE ARRAYS
C     B1,...,B8  OF DIMENSION  N
C     C1,C2  OF DIMENSION  (N,N)
C NOTE THAT IN THE SUBROUTINE DEFINITION STATEMENT
C THESE ARRAYS ARE NAMED ACCORDING TO THEIR INTERNAL USE.
C THE OUTPUT CONTENT OF THESE ARRAYS SHOULD NOT CONCERN
C THE AVERAGE USER.
C THE INTERESTED USER WILL FIND ON OUTPUT IN  B1, B3, B5, B6
C THE FINAL VALUES OF THE VECTORS  V , VFU , XTP , WF
C ( V  AND  VFU  ARE DESCRIBED IN SUBROUTINE  EVOLV , AND  XTP
C IN SUBROUTINE  JCBN ).
C
C CALL TO  DAFNE
C
C THE PROGRAM CALLING  DAFNE  MUST SET THE INPUT
C CALL PARAMETERS  N , X , FMIN , PRF , IDER , NU , NO .
C
C RETURN FROM DAFNE
C
C THE CALLING PROGRAM RECEIVES FROM  DAFNE  THE OUTPUT
C PARAMETERS   X , IDA , NFEV , NJEV , NSAC , NST , IER .
C
C SUBPROGRAMS CALLED  :  EVOLV , NEWF .
C
C THE USER MUST PROVIDE THE FOLLOWING TWO SUBPROGRAMS
C (FORTRAN IMPLICIT TYPE DEFINITION IS USED.
C ALL NON-INTEGER PARAMETERS ARE DOUBLE-PRECISION)
C
C SUBROUTINE  VECF (N,X,VFU,NR,IER)
C THE SUBROUTINE  VECF  MUST RETURN IN THE N-VECTOR  VFU
C THE (NON-RESCALED) VALUES OF THE FUNCTIONS OF SYSTEM (1)
C COMPUTED AT THE POINT  X .
C NR  IS THE (INPUT) NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     MUST BE USED TO INDICATE THE FACT THAT  X  IS OUTSIDE A SAFETY
C     REGION DEFINED BY THE USER).
C
C SUBROUTINE  JCBNA (N,X,C,NR,IDA,IER)
C THE SUBROUTINE  JCBNA  IS DESIGNED TO GIVE THE USER
C THE OPPORTUNITY OF PROVIDING THE ANALYTICAL JACOBIAN.
C IN THIS CASE THE USER MUST CODE THE SUBROUTINE ACCORDING
C TO THE FOLLOWING SPECIFICATIONS.
C THE SUBROUTINE  JCBNA  MUST RETURN IN THE (N,N) ARRAY  C
C THE JACOBIAN OF THE FUNCTIONS OF SYSTEM (1)
C COMPUTED AT THE POINT  X .
C NR  IS THE (INPUT) NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C IDA  MUST BE SET EQUAL TO 1.
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     MUST BE USED TO INDICATE THE FACT THAT  X  IS OUTSIDE A SAFETY
C     REGION DEFINED BY THE USER).
C IF THE USER DOES NOT WANT TO GIVE THE ANALYTICAL JACOBIAN,
C THE SUBROUTINE  JCBNA  MUST ONLY SET  IDA  EQUAL TO 0.
C
      DOUBLE PRECISION AMU,AMU0,BE,BE0,BUF,C,CS,FMIN,FP,FPS,H,H0
      DOUBLE PRECISION S,PRF,V,VFU,V0,WF,X,XTP
      DIMENSION X(N),FMIN(N),V(N),FP(N),VFU(N),BUF(N),XTP(N),WF(N)
      DIMENSION FPS(N),S(N),C(N,N),CS(N,N)
      DIMENSION NSC(8)
C
C NTRIAL  IS THE DIMENSION OF THE ARRAY  NSC  CONTAINING
C     THE MAX. NUMBER OF STEPS FOR THE SUCCESSIVE CALLS TO  EVOLV
      DATA NTRIAL/8/
      DATA NSC(1),NSC(2),NSC(3),NSC(4),NSC(5),NSC(6),NSC(7),NSC(8)
     1    /  1,     2,     4,     8,    16,     32,   64,    128 /
C
C INITIAL VALUE OF TIME INTEGRATION STEPLENGTH
      H0 = 1.D0
C INITIAL VALUES OF MASS AND FRICTION COEFFICIENTS
      AMU0 = 1.D0
      BE0 = 1.D0
C INITIAL VALUE FOR THE VELOCITIES
      V0 = 0.D0
C
C INITIALIZE FUNCTION WEIGHTS
      DO 10 I = 1,N
         WF(I) = 1.D0
  10     CONTINUE
C
C INITIALIZE TIME INTEGRATION STEPLENGTH
      H = H0
C
C INITIALIZE COUNTERS
      NFEV = 0
      NJEV = 0
      NSAC = 0
      NST = 0
C
C START SEQUENCE OF CALL TO  EVOLV
      DO 30 ITRIAL = 1,NTRIAL
C
C        SET MAX. NUMBER OF STEPS FOR NEXT CALL TO  EVOLV
         NIT = NSC(ITRIAL)
C
C        INITIALIZE VELOCITIES, MASS, AND FRICTION COEFFICIENT
         DO 20 I = 1,N
            V(I) = V0
  20        CONTINUE
         AMU = AMU0
         BE = BE0
C
C        CALL THE BASIC EQUATION-SOLVING SUBROUTINE
         CALL EVOLV (N,X,V,AMU,BE,H,FMIN,NIT,NR,II,IER,
     1       PRF,IDER,NU,NO,IDA,NFEV,NJEV,
     2       FP,C,VFU,BUF,XTP,WF,FPS,CS,S)
C
C        UPDATE COUNTERS
         NSAC = NSAC+II
         NST = NST+NR
C
         IF (IER.EQ.0) GO TO 60
         IF (IER.EQ.4.OR.IER.EQ.5) GO TO 40
C
C        COMPUTE FUNCTION WEIGHTS
         NRP = NR+1
         NFEV = NFEV+1
         CALL NEWF (N,X,WF,PRF,IDER,NU,NO,IDA,
     1       XTP,VFU,C,FP,BUF,NRP,NJEV,IER)
         IF (IER.LT.0) GO TO 50
C
  30     CONTINUE
C
C ERROR EXIT (MAX. NUMBER OF ATTEMPTS REACHED)
      IER = 1
      RETURN
C
C ERROR EXIT (FATAL FAILURE IN COMPUTING THE FUNCTION VALUES
C AT THE INITIAL POINT OR THE DERIVATIVES)
  40  CONTINUE
  50  CONTINUE
      IER = 2
      RETURN
C
C NORMAL EXIT
  60  CONTINUE
      IER = 0
C
      RETURN
      END
      SUBROUTINE CHAH (NSR,VF,GA,IER,II,NR)
C
C THE SUBROUTINE  CHAH  COMPUTES THE VARIATION OF THE TIME INTEGRATION
C STEPLENGTH, AND MAY DECIDE TO REJECT THE LAST STEP PERFORMED.
C
C DESCRIPTION OF PARAMETERS
C
C NSR  IS THE NUMBER OF STEPS FOR WHICH TOTAL ENERGY IS REMEMBERED.
C VF  IS AN ARRAY OF LENGTH  NSR  CONTAINING THE LAST VALUES
C     OF THE TOTAL ENERGY.
C GA  IS THE (OUTPUT) MULTIPLICATIVE FACTOR RETURNED TO UPDATE
C     THE TIME INTEGRATION STEPLENGTH.
C IER  IS AN (INPUT OR OUTPUT) ERROR INDICATOR.
C     A NEGATIVE VALUE INDICATES THAT THE LAST STEP IS TO BE
C     REJECTED.
C II  IS THE NUMBER OF TIME INTEGRATION STEPS.
C NR  IS THE TOTAL NUMBER OF ATTEMPTED TIME INTEGRATION STEPS,
C     INCLUDING THE STEPS ATTEMPTED BUT REJECTED.
C
      DOUBLE PRECISION GA,GAB,GAMAX,GAMIN,GAR,ONE,TOL,UGAMAX,UGAR,VF
      DOUBLE PRECISION GAEPS
      DIMENSION VF(NSR)
C
      DATA GAMIN/2.D-2/
      DATA ILB/-20/
      DATA GAMAX /1.D+2/
      DATA GAB /25.D0/
      DATA ONE /1.D0/
      DATA TOL/1.D-10/
      DATA GAR/1.1D0/
      DATA ISR/0/
      DATA GAEPS/1.D-10/
C
      IF (NR.EQ.1) ILB = -20
      IF (NR.EQ.1) ISR = 0
      IF (IER.LT.0) GO TO 80
      UGAMAX = GAMAX
      IF (II.GT.NSR) UGAMAX = GAB
      IF (ILB+5.GT.II) UGAMAX = ONE
      UGAR = DMIN1(GAR,UGAMAX)
      IP = 0
      IM = 0
      DO 10 I = 2,NSR
         IF (IP.EQ.0.AND.VF(I).GT.VF(1)*(1.D0+TOL)) IP = I-1
         IF (IM.EQ.0.AND.VF(I).LT.VF(1)*(1.D0-TOL)) IM = I-1
         IF (IP.NE.0.AND.IM.NE.0) GO TO 20
  10     CONTINUE
  20  CONTINUE
      IF (IP.EQ.1.AND.IM.EQ.0) GO TO 70
      IF (IM.EQ.1) GO TO 30
      IF (IP.EQ.1) GO TO 40
      IF (IP.GT.1.AND.IM.GT.1) GO TO 50
      IF (IP.GT.1) GO TO 60
      GA = UGAR
      GO TO 90
  30  CONTINUE
      IF (IP.EQ.0) GO TO 80
      GA = ONE
      IF ((IP-1)**2.GT.II) GA = GAMIN
      GO TO 90
  40  CONTINUE
      GA = ONE
      IF (II.GT.NSR) GA = UGAR
      IF ((IM-1)**2.GT.II*4) GA = UGAMAX
      GO TO 90
  50  CONTINUE
      GA = UGAMAX
      GO TO 90
  60  CONTINUE
      GA = UGAMAX
      GO TO 90
  70  CONTINUE
      GA = UGAMAX
      GO TO 90
  80  CONTINUE
      ISR = ISR+1
      IER = -3
      GA = GAMIN
      IF (ISR.GE.3) GA = GAEPS
      ILB = II
      RETURN
  90  CONTINUE
      IER = 1
      ISR = 0
      RETURN
      END
      SUBROUTINE CREXTP (N,X,SPRR,XTP)
C
C THE SUBROUTINE  CREXTP  INITIALIZES THE ARRAY  XTP  (SEE BELOW)
C
C DESCRIPTION OF PARAMETERS
C
C N  AND  X  ARE DEFINED IN SUBROUTINE  EVOLV .
C SPRR  IS A VARIABLE WHOSE VALUE IS COMPUTED BY SUBROUTINE
C     JCBN , DEPENDING ON THE ESTIMATED NUMERICAL PRECISION.
C XTP  IS AN ARRAY OF LENGTH  N  USED BY
C     SUBROUTINE  JCBN  TO COMPUTE THE STEPLENGTHS USED
C     FOR COMPUTING THE FINITE-DIFFERENCES JACOBIAN.
C
      DOUBLE PRECISION XTP,SPRR,X
      DIMENSION X(N),XTP(N)
C
      DO 10 I = 1,N
         XTP(I) = DMAX1(SPRR,DABS(X(I)))
  10     CONTINUE
      RETURN
      END
      SUBROUTINE DERIV (N,X,VFU,FP,C,WF,NR,NJEV,IER,
     1    PRF,IDER,NU,NO,IDA,BUF,XTP)
C
C THE SUBROUTINE  DERIV RETURNS IN THE N-VECTOR  FP
C THE GRADIENT OF THE WEIGHTED SUM
C     FW = W1*F1(X)**2 + ... + WN*FN(X)**2
C OF THE SQUARES OF THE LEFT-HAND SIDE FUNCTIONS OF SYSTEM (1)
C (SEE  EVOLV ), AND IN THE (N BY N) MATRIX  C  THE (APPROXIMATE)
C HESSIAN MATRIX
C
C          T
C     2 * J  * W * J
C
C OF  FW , COMPUTED AT THE POINT  X , WHERE  W  IS
C THE DIAGONAL MATRIX OF THE RESCALING
C WEIGHTS, AND  J  IS THE JACOBIAN MATRIX OF THE (NON-RESCALED)
C FUNCTIONS IN (1), OBTAINED BY CALLING THE SUBROUTINE  JCBN .
C VFU  IS AN INPUT ARRAY OF LENGTH  N  CONTAINING THE VALUES OF
C     THE FUNCTIONS IN SYSTEM (1).
C WF  IS THE N-VECTOR CONTAINING THE RESCALING WEIGHTS.
C NR  IS THE (INPUT) NUMBER OF THE TIME INTEGRATION STEP IN
C     PROGRESS (COUNTING ALSO THE REJECTED STEPS).
C NJEV  IS THE COUNTER OF THE JACOBIAN EVALUATIONS
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     INDICATES A FAILURE IN COMPUTING THE JACOBIAN).
C PRF , IDER , NU , NO , IDA  ARE DESCRIBED WITHIN THE
C     SUBROUTINE  DAFNE .
C BUF  IS A WORKING STORAGE ARRAY OF LENGTH  N
C XTP  IS DESCRIBED IN THE SUBROUTINE  JCBN .
C
C SUBPROGRAMS CALLED  :  JCBN
C
      DOUBLE PRECISION BUF,C,FP,PRF,VFU,WF,X,XTP
      DIMENSION X(N),C(N,N),FP(N),VFU(N),BUF(N),XTP(N),WF(N)
C
C COMPUTE THE JACOBIAN MATRIX
      CALL JCBN (N,X,VFU,C,FP,BUF,XTP,PRF,IDER,NU,NO,IDA,
     1    NR,NJEV,IER)
      IF (IER.LT.0) RETURN
C
C COMPUTE THE GRADIENT VECTOR  FP
      DO 20 I = 1,N
         FP(I) = 0.D0
         DO 10 J = 1,N
            FP(I) = FP(I)+VFU(J)*C(J,I)*WF(J)
  10        CONTINUE
         FP(I) = FP(I)*2.D0
  20     CONTINUE
C
      DO 70 K = 1,N
C
C     COMPUTE COLUMN  K  OF THE APPROX. HESSIAN AND
C     OVERWRITE IT ON COLUMN  K  OF THE JACOBIAN.
C
C        COMPUTE THE BELOW-DIAGONAL PART OF COLUMN  K ...
         DO 40 J = K,N
            BUF(J) = 0.D0
            DO 30 I = 1,N
               BUF(J) = BUF(J)+C(I,J)*C(I,K)*WF(I)
  30           CONTINUE
  40        CONTINUE
C
C        ... AND OVERWRITE IT ON COLUMN  K  OF THE JACOBIAN
         DO 50 J = K,N
            C(J,K) = BUF(J)*2.D0
  50        CONTINUE
C
C        FILL BY SYMMETRY THE ABOVE-DIAGONAL PART OF COLUMN  K
         DO 60 J = 1,K
            C(J,K) = C(K,J)
  60        CONTINUE
C
  70     CONTINUE
      RETURN
      END
      SUBROUTINE DLIN (N,FP,C,FPS,CS,S,IER)
C
C THE SUBROUTINE DLIN COMPUTES THE MINIMUM-NORM
C LEAST-SQUARES SOLUTION OF THE LINEAR SYSTEM
C WITH POSSIBLY SINGULAR  N BY N  COEFFICIENT MATRIX  C
C AND RIGHT-HAND SIDE N-VECTOR  FP ,  AND RETURNS
C IN  FP  THE SOLUTION VECTOR (THE CONTENT OF  C
C MAY BE ALTERED).
C FPS , CS , S  ARE WORKING STORAGE ARRAYS OF LENGTH N , (N,N) , N .
C IER IS AN OUTPUT INDICATOR AS FOLLOWS
C     IER = -1       SOLUTION NOT FOUND
C     IER = 0        SOLUTION FOUND WITHOUT REGULARIZATION
C     IER = 1        SOLUTION FOUND BY REGULARIZATION
C
C SUBPROGRAMS CALLED  :  QMS , MOVM , DLINX .
C
      DOUBLE PRECISION C,FP
      DOUBLE PRECISION CS,D,DEPS1,DEPS2
      DOUBLE PRECISION DMEM,DMIN,DMAX,DSTEP,FPS,S,SP,SPV,SR
      DIMENSION FP(N),C(N,N)
      DIMENSION FPS(N),CS(N,N),S(N)
C
C THE REGULARIZATION PARAMETER  D  STARTS FROM THE PRECEDING
C VALUE  DMEM  AND IS REPEATEDLY UPDATED BY A FACTOR  DSTEP
C NOT BEYOND THE BOUNDARY VALUES  DMIN , DMAX .
C
C USE THE FIRST VALUES (-60,60) OR THE SECOND ONES (-30,30) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1.D120 (E.G. UNIVAC 1100), OR BELOW
C 1.D120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA DMIN,DMAX/1.D-60,1.D60/
C     DATA DMIN,DMAX/1.D-30,1.D30/
C
      DATA DSTEP/10.D0/
      DATA DMEM/1.D0/
C
C DEPS1  AND  DEPS2  ARE TOLERANCE FACTORS USED BY SUBROUTINE  DLINX .
      DATA DEPS1,DEPS2/1.D-16,1.D-16/
C
      DATA NMAX/1000/
C
C RESCALE BY NORMALIZING MAIN DIAGONAL ELEMENTS TO ONE .
C
C     COMPUTE SCALE FACTORS
      DO 10 I = 1,N
         S(I) = 1.D0/DSQRT(C(I,I))
  10     CONTINUE
C
      DO 30 I = 1,N
C
C        RESCALE COEFFICIENT MATRIX
         DO 20 J = 1,N
            C(I,J) = C(I,J)*S(I)*S(J)
  20        CONTINUE
C
C        RESCALE RIGHT-HAND-SIDE VECTOR
         FP(I) = FP(I)*S(I)
C
  30     CONTINUE
C
C     SAVE  FP  AND  C  IN  FPS  AND  CS
      CALL MOVM (N,FP,C,FPS,CS)
C
C     TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C     SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C     SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C     WITH TOLERANCE PARAMETER  DEPS1 .
      CALL DLINX (N,FPS,CS,DEPS1,IER)
C
      IF (IER.GT.0) IER = 0
      IF (IER.EQ.0) GO TO 140
C
C     THE LINEAR SYSTEM IS SINGULAR
C     START SEARCH FOR REGULARIZATION FACTOR  D
C     STARTING FROM PRECEDING VALUE  DMEM .
      D = DMEM
C
  40  CONTINUE
C
C DECREASE REGULARIZATION PARAMETER
      D = D/DSTEP
C
      DO 50 I = 1,NMAX
C
C        REGULARIZE THE LINEAR SYSTEM  C ,  FP , WITH REGULARIZATION
C        PARAMETER  D , OBTAINING THE LINEAR SYSTEM  FPS ,  CS .
         CALL QMS (N,FP,C,FPS,CS,D)
C
C        TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C        SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C        SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C        WITH TOLERANCE PARAMETER  DEPS2 .
         CALL DLINX (N,FPS,CS,DEPS2,IER)
C
         IF (IER.GE.0.AND.I.EQ.1.AND.D.GT.DMIN) GO TO 40
         IF (IER.GE.0) GO TO 80
         IF (D.GT.DMAX) GO TO 60
C
C        INCREASE THE REGULARIZATION PARAMETER
         D = D*DSTEP
C
  50     CONTINUE
C
  60  CONTINUE
  70  CONTINUE
C
C ERROR EXIT
      IER = -1
      RETURN
C
  80  CONTINUE
C INCREASE THE REGULARIZATION PARAMETER
      DMEM = D*DSTEP
C
C     THE LINEAR SYSTEM IS ESTIMATED NONSINGULAR BY SUBROUTINE  DLINX
C     LOOK FOR LEAST-SQUARES SOLUTION
      SP = 0.D0
C
      DO 110 I = 1,NMAX
         SPV = SP
         SP = 0.D0
         DO 100 J = 1,N
            SR = 0.D0
            DO 90 K = 1,N
               SR = SR+C(J,K)*FPS(K)
  90           CONTINUE
            SR = SR-FP(J)
            SP = SP+SR*SR
 100        CONTINUE
         IF (SP.GT.SPV.AND.I.GT.1) GO TO 120
C
C        INCREASE THE REGULARIZATION PARAMETER
         D = D*DSTEP
         IF (D.GT.DMAX) GO TO 130
C
C        REGULARIZE THE LINEAR SYSTEM  C ,  FP , WITH REGULARIZATION
C        PARAMETER  D , OBTAINING THE LINEAR SYSTEM  FPS ,  CS .
         CALL QMS (N,FP,C,FPS,CS,D)
C
C        TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C        SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C        SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C        WITH TOLERANCE PARAMETER  DEPS2 .
         CALL DLINX (N,FPS,CS,DEPS2,IER)
C
         IF (IER.LT.0) GO TO 70
C
 110     CONTINUE
C
 120  CONTINUE
      D = D/DSTEP
C
C     REGULARIZE THE LINEAR SYSTEM  C ,  FP , WITH REGULARIZATION
C     PARAMETER  D , OBTAINING THE LINEAR SYSTEM  FPS ,  CS .
      CALL QMS (N,FP,C,FPS,CS,D)
C
C     TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C     SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C     SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C     WITH TOLERANCE PARAMETER  DEPS2 .
      CALL DLINX (N,FPS,CS,DEPS2,IER)
C
 130  CONTINUE
      IER = 1
 140  CONTINUE
C
C REVERSE SCALING FOR VECTOR  FP
      DO 150 I = 1,N
         FP(I) = FPS(I)*S(I)
 150     CONTINUE
C
      RETURN
C
      END
      SUBROUTINE DLINX (N,B,A,DEPS,IER)
C
C THE SUBROUTINE  DLINX  ATTEMPTS TO COMPUTE THE SOLUTION
C OF THE LINEAR SYSTEM WITH  N BY N  COEFFICIENT MATRIX  A
C AND RIGHT-HAND SIDE N-VECTOR  B ,  AND RETURNS
C IN  B  THE SOLUTION VECTOR (THE CONTENT OF  A
C MAY BE ALTERED).
C DEPS  IS AN INPUT TOLERANCE FACTOR USED TO DETECT FAILURE.
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE INDICATES FAILURE)
C
      DOUBLE PRECISION A,AM,DEPS,P,PI,B,Z,ERR
C
      DIMENSION A(N,N),B(N)
      IF(N.LE.0)GO TO 230
C
      IER = 0
      IRM = 0
      ICM = 0
      P = 0.D0
      DO 30 IR = 1,N
         DO 10 IC = 1,N
            AM = DABS(A(IR,IC))
            IF(AM.LE.P)GO TO 10
            P = AM
            IRM = IR
            ICM = IC
  10        CONTINUE
  30     CONTINUE
      IF(IRM.EQ.0.OR.ICM.EQ.0)GO TO 240
      ERR = DEPS*P
C
      DO 170 K = 1,N
C
C
         IF(P.LE.0.D0)GO TO 230
         IF(IER.NE.0)GO TO 70
         IF(P.GT.ERR)GO TO 70
         IER = K-1
  70     PI = 1.D0/A(IRM,ICM)
         I = IRM - K
         J = ICM - K
C
         Z = PI*B(IRM)
         B(IRM) = B(K)
         B(K) = Z
C
         IF(K.GE.N)GO TO 180
C
         IF(J.LE.0)GO TO 120
         DO 110 IR = K,N
            Z = A(IR,K)
            A(IR,K) = A(IR,ICM)
 110        A(IR,ICM) = Z
C
 120     DO 130 IC = K,N
            Z = PI*A(IRM,IC)
            A(IRM,IC) = A(K,IC)
 130        A(K,IC) = Z
C
         A(K,K) = J
C
         P = 0.D0
         IRS = K + 1
         DO 160 IR = IRS,N
            PI = -A(IR,K)
            IST = K+1
            DO 150 IC = IST,N
               A(IR,IC) = A(IR,IC)+PI*A(K,IC)
               Z = DABS(A(IR,IC))
               IF(Z.LE.P)GO TO 150
               P = Z
               IRM = IR
               ICM = IC
 150           CONTINUE
 160        B(IR) = B(IR)+PI*B(K)
 170     CONTINUE
C
 180  IF (N-1) 230,220,190
 190  IR = N
      IC = N+1
      DO 210 I = 2,N
         IR = IR - 1
         IC = IC - 1
         L = A(IR,IR)+.5D0
         Z = B(IR)
         LL = IR
         DO 200 K = IC,N
            LL = LL+1
 200        Z = Z-A(IR,K)*B(LL)
         K = IR+L
         B(IR) = B(K)
 210     B(K) = Z
 220  RETURN
C
C     ERROR RETURN
 230  IER = -1
      RETURN
 240  CONTINUE
      IER = -2
      RETURN
      END
      SUBROUTINE EVOLV (N,X,V,AMU,BE,H,FMIN,NIT,NR,II,IER,
     1    PRF,IDER,NU,NO,IDA,NFEV,NJEV,
     2    FP,C,VFU,BUF,XTP,WF,FPS,CS,S)
C
C THE SUBROUTINE  EVOLV  ATTEMPTS TO SOLVE A SYSTEM
C OF  N  NONLINEAR EQUATIONS
C     FI(X) = 0,  I = 1,...,N,  X = (X1,...,XN)        (1)
C THE SOLUTION IS SOUGHT AT A ZERO-VALUED MINIMUM OF THE WEIGHTED SUM
C     FW = W1*F1(X)**2 + ... + WN*FN(X)**2             (2)
C OF THE SQUARES OF THE LEFT-HAND SIDE FUNCTIONS OF SYSTEM (1).
C THE MINIMIZER IS SOUGHT AS THE ASYMPTOTIC VALUE (AS TIME  T  GOES TO
C INFINITY) OF THE SOLUTION TRAJECTORY  X(T) = (X1(T),...,XN(T))  OF
C THE ASSOCIATED DYNAMICAL SYSTEM DEFINED BY THE DIFFERENTIAL EQUATIONS
C     AMU*(DV/DT) = -BE*V -(DFW/DX)                    (3)
C WHERE  V = (V1,...,VN)  AND  (DV/DT)  ARE THE VECTORS OF THE
C VELOCITIES  VI(T) = DXI/DT  AND OF THE ACCELERATIONS  DVI/DT,
C I = 1,...,N,  (DFW/DX)  IS THE GRADIENT OF  FW  W.R.T.  X ,  AND
C AMU  AND  BE  ARE THE (SCALAR) MASS AND FRICTION COEFFICIENT.
C THE NUMERICAL INTEGRATION OF (3) IS PERFORMED WITH THE
C LAMBERT-SIGURDSSON LINEARLY IMPLICIT SCHEME
C WITH VARIABLE TIME INTEGRATION STEPLENGTH AND WITH
C THE POSSIBILITY OF REJECTING A STEP WHICH IS ESTIMATED
C NON SATISFACTORY.
C NOTE THAT THE WEIGHTING OF THE SQUARES OF THE FUNCTIONS IN  (2)
C IS EQUIVALENT TO RESCALING THE FUNCTION VALUES WITH FACTORS SQRT(W1),
C ..., SQRT(WN).
C
C DESCRIPTION OF PARAMETERS
C
C N  IS THE NUMBER OF NONLINEAR EQUATIONS.
C X,V  ARE VECTORS OF LENGTH  N  CONTAINING X-VARIABLES AND VELOCITIES.
C     (INITIAL VALUES ON INPUT, FINAL ON OUTPUT. INITIAL AND FINAL  X
C     ARE THE INITIAL AND FINAL ESTIMATES OF THE SOLUTION OF SYSTEM (1))
C AMU  IS THE VALUE OF THE  MASS COEFFICIENT.
C BE  IS THE FRICTION COEFFICIENT
C     (INITIAL VALUE ON INPUT, FINAL VALUE ON OUTPUT).
C H  IS THE TIME STEPLENGTH FOR THE NUMERICAL INTEGRATION OF SYSTEM (3)
C     (INITIAL VALUE ON INPUT, FINAL VALUE ON OUTPUT).
C FMIN  IS  A VECTOR OF DIMENSION  N  CONTAINING THE THRESHOLD
C     VALUES FOR THE FUNCTION VALUES  FI(X) , I = 1,...,N  (SEE BELOW).
C NIT  IS THE MAXIMUM ALLOWED NUMBER OF TIME INTEGRATION STEPS.
C NR  IS THE TOTAL NUMBER OF ATTEMPTED TIME INTEGRATION STEPS
C     (I.E. INCLUDING STEPS ATTEMPTED BUT REJECTED)
C II  IS THE NUMBER OF ACCEPTED STEPS.
C IER  IS AN OUTPUT INDICATOR OF THE STOPPING CONDITION, AS FOLLOWS
C     IER = 0   CONVERGENCE REACHED  AT STEP  II ,  I.E.
C               ABS(FI(X)).LT.FMIN(I) , I = 1,...,N
C     IER = 1   NUMBER OF ACCEPTED INTEGRATION STEPS  .EQ.NIT
C     IER = 2   ESTIMATED RATE OF CONVERGENCE TOO SLOW
C     IER = 3   ATTEMPT TO REDUCE THE TIME INTEGRATION STEPLENGTH
C               BELOW THE ALLOWED MINIMUM.
C     IER = 4   FATAL FAILURE IN COMPUTING THE DERIVATIVES.
C     IER = 5   FATAL FAILURE IN COMPUTING THE FUNCTION VALUES AT THE
C               INITIAL POINT.
C THE PARAMETERS  PRF , IDER , NU , NO , IDA
C ARE USED ONLY IN THE COMPUTATION OF THE NUMERICAL JACOBIAN
C ( SUBROUTINE  JCBN ), AND ARE DESCRIBED WITHIN THE SUBROUTINE  DAFNE .
C THE FOLLOWING 2 PARAMETERS ( NFEV , NJEV ) ARE OUTPUT COUNTERS
C     NFEV  OF FUNCTION EVALUATIONS (EXCLUDING THOSE USED
C           FOR JACOBIAN EVALUATION).
C     NJEV  OF JACOBIAN EVALUATIONS.
C FP,C  ARE WORKING STORAGE ARRAYS OF DIMENSION  N  AND  (N,N) .
C VFU  IS A VECTOR OF DIMENSION  N  USED FOR INTERNAL STORAGE
C     OF THE VALUES OF THE FUNCTIONS IN (1).
C BUF  AND  XTP  ARE WORKING STORAGE ARRAY OF LENGTH  N .
C WF  IS THE VECTOR OF THE WEIGHTS DESCRIBED IN SUBROUTINE  DAFNE .
C FPS , CS , S  ARE WORKING STORAGE ARRAYS OF DIMENSION  N , (N,N) , N
C     NEEDED BY SUBROUTINE  DLIN .
C
C SUBPROGRAMS CALLED  :  CHAH , FUNZ , DERIV , DLIN .
C
      DOUBLE PRECISION AMU,B,BE,BEMAX,BUF,C,CS,EMAX,EMIN
      DOUBLE PRECISION ET,F,FMIN,FP,FPS,GA,H,HMAX,HMIN
      DOUBLE PRECISION PRF,RAPMIN,S,V,VF,VFU,WF,X,XTP
      DOUBLE PRECISION FUNZ
      DIMENSION X(N),V(N),FMIN(N),FP(N),C(N,N),VFU(N),BUF(N)
      DIMENSION XTP(N),WF(N),FPS(N),CS(N,N),S(N)
      DIMENSION VF(30)
C
      DATA RAPMIN/1.0001D0/
C
C EXTREME VALUES FOR  H  AND  BE
C USE THE FIRST VALUES (-50,50,30) OR THE SECOND ONES (-25,25,17) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1.D120 (E.G. UNIVAC 1100), OR BELOW
C 1.D120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA HMIN,HMAX/1.D-50,1.D+50/,BEMAX/1.D30/
C     DATA HMIN,HMAX/1.D-25,1.D+25/,BEMAX/1.D17/
C
C NSR IS THE NUMBER OF STEPS FOR WHICH TOTAL ENERGY IS REMEMBERED
C NSR MUST BE .LE. DIMENSION OF VF
      DATA NSR/30/
C
      DATA ET/0.D0/
C
C INITIALIZE
      II = 0
      NR = 0
      F = FUNZ(N,X,VFU,NR,NFEV,WF,IER)
      IF (IER.LT.0) GO TO 190
      DO 10 ID = 1,NSR
         VF(ID) = F*DBLE(FLOAT(ID))**4
  10     CONTINUE
      H = DMAX1(H,HMIN)
C
C START NUMERICAL INTEGRATION LOOP  ( NIT  STEPS)
C
      DO 150 IT = 1,NIT
         II = IT
         DO 20 ID = 2,NSR
            IID = NSR+1-ID
            VF(IID+1) = VF(IID)
  20        CONTINUE
C
C        START OF TIME INTEGRATION STEP
C
  30     CONTINUE
         NR = NR+1
C
C        COMPUTE GRADIENT AND HESSIAN
         CALL DERIV (N,X,VFU,FP,C,WF,NR,NJEV,IER,
     1       PRF,IDER,NU,NO,IDA,BUF,XTP)
         IF (IER.LT.0) GO TO 200
C
         DO 50 L = 1,N
            DO 40 J = 1,N
               C(L,J) = C(L,J)*H
  40           CONTINUE
            C(L,L) = C(L,L)+BE+AMU/H
            FP(L) = -FP(L)*H+V(L)*AMU
  50        CONTINUE
C
C        SOLUTION OF LINEAR SYSTEM
         CALL DLIN (N,FP,C,FPS,CS,S,IER)
C
C        SAVE AND UPDATE X-VARIABLES AND VELOCITIES
         DO 60 L = 1,N
            C(1,L) = X(L)
            BUF(L) = VFU(L)
            B = V(L)
            V(L) = FP(L)/H
            X(L) = X(L)+FP(L)
            FP(L) = B
  60        CONTINUE
         IF (IER.LT.0) GO TO 100
C
C        COMPUTE AND SAVE TOTAL ENERGY
         F = FUNZ(N,X,VFU,NR,NFEV,WF,IER)
         IF (IER.LT.0) GO TO 110
         ET = F
         DO 70 J = 1,N
            ET = ET+AMU*0.5D0*V(J)**2
  70        CONTINUE
         VF(1) = ET
C
C
C        END OF TIME INTEGRATION STEP
C
C        TEST CONVERGENCE
         DO 80 J = 1,N
            IF (DABS(VFU(J)).GT.FMIN(J)) GO TO 90
  80        CONTINUE
         GO TO 160
  90     CONTINUE
C
 100     CONTINUE
 110     CONTINUE
C
C        COMPUTE TIME INTEGRATION STEPLENGTH
         CALL CHAH (NSR,VF,GA,IER,II,NR)
         H = H*GA
         IF (H.LT.HMIN) GO TO 170
         H = DMIN1(H,HMAX)
C
C        IF STEP REJECTED RESTORE PRECEDING CONDITIONS
         IF (IER.GE.0) GO TO 130
         DO 120 L = 1,N
            X(L) = C(1,L)
            V(L) = FP(L)
            VFU(L) = BUF(L)
 120        CONTINUE
         GO TO 30
C
 130     CONTINUE
C
C        TEST IF CONVERGENCE TOO SLOW
         EMAX = ET
         EMIN = ET
         DO 140 ID = 2,NSR
            EMAX = DMAX1(VF(ID),EMAX)
            EMIN = DMIN1(VF(ID),EMIN)
 140        CONTINUE
         IF (EMAX.LT.EMIN*RAPMIN) GO TO 180
C
C        COMPUTE FRICTION COEFFICIENT
         IF (II.GT.10) BE = DMIN1(BEMAX,BE+BE)
C
 150     CONTINUE
C
C END OF NUMERICAL INTEGRATION LOOP  ( NIT  STEPS)
C
C EXIT CONDITIONS
C
C MAX  NUMBER OF TIME INTEGRATION STEPS
      IER = 1
      RETURN
C
C CONVERGENCE
 160  CONTINUE
      IER = 0
      RETURN
C
C ATTEMPT TO REDUCE THE TIME INTEGRATION STEP BELOW THE ALLOWED MINIMUM
 170  CONTINUE
      IER = 3
      RETURN
C
C ESTIMATED RATE OF CONVERGENCE TOO SLOW
 180  CONTINUE
      IER = 2
      RETURN
C
C FATAL FAILURE IN COMPUTING FUNCTION VALUES AT INITIAL POINT.
 190  CONTINUE
      IER = 5
      RETURN
C
C FATAL FAILURE IN COMPUTING DERIVATIVES.
 200  CONTINUE
      IER = 4
      RETURN
C
      END
      DOUBLE PRECISION FUNCTION FUNZ (N,X,VFU,NR,NFEV,WF,IER)
C
C THE FUNCTION  FUNZ  RETURNS AS ITS
C FUNCTION VALUE THE WEIGHTED SUM
C     FW = W1*F1(X)**2 + ... + WN*FN(X)**2
C OF THE SQUARES OF THE LEFT-HAND SIDE FUNCTIONS OF SYSTEM (1)
C (SEE SUBROUTINE  EVOLV ), WHICH AMOUNTS TO RESCALE THE
C FUNCTION VALUES WITH FACTORS  SQRT(W1),..., SQRT(WN),
C AND RETURNS IN THE VECTOR  VFU  OF DIMENSION  N  THE VALUES
C OF THE (NON-RESCALED) FUNCTIONS IN (1), COMPUTED AT  X
C BY CALLING THE SUBROUTINE  VECF .
C NR  IS THE NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C NFEV  IS THE COUNTER OF FUNCTION EVALUATIONS (EXCLUDING
C     THOSE USED FOR JACOBIAN EVALUATION.
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     INDICATES A FAILURE IN COMPUTING A FUNCTION VALUE).
C THE WEIGHTED SUM OF THE SQUARES IS COMPUTED BY MEANS OF
C THE WEIGHTS  W1,...,WN  CONTAINED IN THE N-VECTOR  WF .
C
C SUBPROGRAMS CALLED  :  VECF
C   ( VECF  MUST BE SUPPLIED BY THE USER, SEE SUBROUTINE  DAFNE ).
C
      DOUBLE PRECISION F,VFU,WF,X
      DIMENSION X(N),VFU(N),WF(N)
C
C UPDATE FUNCTION EVALUATIONS COUNTER
      NFEV = NFEV+1
C
C COMPUTE FUNCTION VALUES
      CALL VECF (N,X,VFU,NR,IER)
      IF (IER.LT.0) RETURN
C
C COMPUTE THE WEIGHTED SUM OF THE SQUARES
      F = 0.D0
      DO 10 I = 1,N
         F = F+WF(I)*VFU(I)**2
  10     CONTINUE
C
      FUNZ = F
      RETURN
      END
      SUBROUTINE JCBN (N,X,VFU,C,FP,BUF,XTP,PRF,IDER,NU,NO,IDA,
     1    NR,NJEV,IER)
C
C THE SUBROUTINE  JCBN  COMPUTES THE JACOBIAN MATRIX
C OF THE SYSTEM (1) (SEE SUBROUTINE  EVOLV ).
C
C DESCRIPTION OF PARAMETERS
C
C THE PARAMETERS  N , X , VFU , FP , BUF  ARE DESCRIBED WITHIN THE
C SUBROUTINE  EVOLV .
C
C C  IS AN OUTPUT N BY N ARRAY CONTAINING THE COMPUTED JACOBIAN
C
C XTP  IS AN N-VECTOR INITIALIZED BY SUBROUTINE  CREXTP ,
C     UPDATED BY SUBROUTINE  NEWXTP , AND USED HERE FOR
C     COMPUTING THE FINITE-DIFFERENCES JACOBIAN.
C
C PRF  IS AN (INPUT) ROUGH ESTIMATE OF THE AVERAGE RELATIVE
C     ERROR AFFECTING THE FUNCTION VALUES SUPPLIED BY SUBROUTINE  VECF .
C     IF THIS ERROR IS ESTIMATED TO BE NOT SEVERAL ORDERS OF MAGNITUDE
C     GREATER THAN MACHINE PRECISION, THE USER MAY SIMPLY PUT  PRF = 0 .
C
C IDER  IS AN INPUT CONTROL VARIABLE TO REQUEST THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDER.EQ.0   ANALYTICAL COMPUTATION BY SUBROUTINE  JCBNA
C     IDER.EQ.1   COMPUTATION BY FORWARD FINITE DIFFERENCES
C     IDER.EQ.2   COMPUTATION BY CENTRAL FINITE DIFFERENCES
C
C NU  AND  NO  ARE NUMBERS DEFINING THE (POSSIBLY) BANDED
C STRUCTURE OF THE JACOBIAN MATRIX, USED FOR FINITE-DIFFERENCES
C JACOBIAN COMPUTATION
C     NU IS THE N. OF NON-ZERO JACOBIAN DIAGONALS UNDER MAIN DIAGONAL
C     NO IS THE N. OF NON-ZERO JACOBIAN DIAGONALS ABOVE MAIN DIAGONAL
C     (FOR A NON-BANDED MATRIX  NU = NO = N-1 )
C
C IDA  IS AN (OUTPUT) INDICATOR OF THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDA.EQ.1   THE DERIVATIVES WERE COMPUTED ANALYTICALLY
C                (AS REQUESTED)
C     IDA.EQ.0   THE DERIVATIVES WERE COMPUTED BY FINITE DIFFERENCES.
C                (CENTRAL FINITE DIFFERENCES IF REQUESTED BY  IDER = 2,
C                FORWARD FINITE DIFFERENCES IF REQUESTED BY  IDER = 1,
C                OR IF REQUESTED ANALYTICAL DERIVATIVES WERE NOT CODED
C                IN SUBROUTINE  JCBNA ).
C
C NR  IS THE NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C NJEV  IS THE COUNTER OF THE JACOBIAN EVALUATIONS
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     INDICATES A FAILURE IN COMPUTING A FUNCTION VALUE).
C
C SUBPROGRAMS CALLED  :  PRMACH , CREXTP , NEWXTP , JCBNA
C ( JCBNA  MUST BE SUPPLIED BY THE USER, SEE SUBROUTINE  DAFNE ).
C
      DOUBLE PRECISION BUF,C,E,EE,XTP,FP
      DOUBLE PRECISION PR,PRF,PRM,PRR,SPRR
      DOUBLE PRECISION VFU,X,XX
      DIMENSION X(N),FP(N),C(N,N),BUF(N),VFU(N),XTP(N)
      DATA PRM/0.D0/
C
C UPDATE JACOBIAN EVALUATIONS COUNTER
      NJEV = NJEV+1
C
      IDA = 0
      IER = 0
C
C TRY TO COMPUTE ANALYTICAL JACOBIAN, IF REQUIRED
      IF (IDER.EQ.0) CALL JCBNA (N,X,C,NR,IDA,IER)
C
C ERROR EXIT
      IF (IER.LT.0) RETURN
C
C EXIT DUE TO SUCCESSFUL COMPUTATION OF THE
C ANALYTICAL JACOBIAN
      IF (IDA.EQ.1) RETURN
C
C COMPUTE FINITE-DIFFERENCE JACOBIAN (IF REQUIRED,
C OR DUE TO AN UNSUCCESSFUL ATTEMPT TO COMPUTE ANALYTICAL
C JACOBIAN).
C
      IDERM = IDER
      IF (IDER.GT.0) GO TO 10
C
C COMPUTE FORWARD FINITE-DIFFERENCE JACOBIAN IF
C ANALYTICAL JACOBIAN WAS NOT AVAILABLE.
      IDER = 1
C
  10  CONTINUE
C
C INITIALIZATION (PERFORMED ONLY IN THE FIRST
C TIME INTEGRATION STEP).
C
      IF (NR.GT.1) GO TO 20
C
C ESTIMATE MACHINE PRECISION  PRM
      IF (PRM.LE.0.D0) CALL PRMACH (PRM)
C
C ESTIMATE ACTUAL WORKING PRECISION  PR
      PR = DMAX1(PRF,PRM)
      PRR = PR
      IF (IDER.EQ.2) PRR = PR**(2.D0/3.D0)
      SPRR = DSQRT(PRR)
C
C INITIALIZE VECTOR  XTP
      CALL CREXTP (N,X,SPRR,XTP)
C
  20  CONTINUE
C END OF INITIALIZATION
C
      NRP = NR+1
C
C TEST APPLICABILITY OF BANDED JACOBIAN ALGORITHM
      ND = NO+NU+1
      IF (ND.LT.N) GO TO 70
C
C NON-BANDED JACOBIAN ALGORITHM
      DO 60 I = 1,N
C
C        COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
         E = SPRR*DMAX1(DABS(X(I)),XTP(I))
         E = DSIGN(E,X(I))
C
         XX = X(I)
         X(I) = X(I)+E
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,X,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         EE = 1.D0/E
C
C        COMPUTE FORWARD FINITE-DIFFERENCE DERIVATIVES
         DO 30 J = 1,N
            C(J,I) = (BUF(J)-VFU(J))*EE
  30        CONTINUE
C
         IF (IDER.EQ.1) GO TO 50
C
C        COMPUTATION FOR CENTRAL DIFFERENCE DERIVATIVES
         X(I) = XX-E
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,X,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         DO 40 J = 1,N
            C(J,I) = .5D0*(C(J,I)-(BUF(J)-VFU(J))*EE)
  40        CONTINUE
C
  50     CONTINUE
         X(I) = XX
C
  60     CONTINUE
C END OF COMPUTATIONS FOR CENTRAL DIFFERENCE DERIVATIVES
C
      GO TO 160
C
C ALGORITHM FOR BANDED JACOBIAN
  70  CONTINUE
C
      DO 80 I = 1,N
         FP(I) = X(I)
  80     CONTINUE
C
      DO 150 I = 1,ND
C
         DO 90 J = I,N,ND
C
C           COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
            E = SPRR*DMAX1(DABS(X(J)),XTP(J))
            E = DSIGN(E,X(J))
C
            FP(J) = X(J)+E
  90        CONTINUE
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,FP,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         DO 110 J = I,N,ND
C
C           COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
            E = SPRR*DMAX1(DABS(X(J)),XTP(J))
            E = DSIGN(E,X(J))
C
            FP(J) = X(J)-E*DBLE(FLOAT(IDER-1))
            EE = 1.D0/E
            K1 = MAX0(1,J-NO)
            K2 = MIN0(N,J+NU)
C
C           FORWARD FINITE DIFFERENCES
            DO 100 K = K1,K2
               C(K,J) = (BUF(K)-VFU(K))*EE
 100           CONTINUE
C
 110        CONTINUE
C
         IF (IDER.EQ.1) GO TO 140
C
C        CENTRAL FINITE DIFFERENCES
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,FP,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         DO 130 J = 1,N,ND
C
C           COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
            E = SPRR*DMAX1(DABS(X(J)),XTP(J))
            E = DSIGN(E,X(J))
C
            FP(J) = X(J)
            EE = 1.D0/E
            K1 = MAX0(1,J-NO)
            K2 = MIN0(N,J+NU)
C
            DO 120 K = K1,K2
               C(K,J) = .5D0*(C(K,J)-(BUF(K)-VFU(K))*EE)
 120           CONTINUE
 130        CONTINUE
 140     CONTINUE
 150     CONTINUE
 160  CONTINUE
C
C UPDATE VECTOR  XTP
      CALL NEWXTP (N,X,VFU,C,XTP)
C
      IDER = IDERM
      RETURN
      END
      SUBROUTINE MOVM (N,FP,C,FPS,CS)
C
C THE SUBROUTINE  MOVM  SAVES THE ARRAYS  FP  AND  C
C (OF DIMENSION  N  AND  (N,N) ) IN THE ARRAYS  FPS  AND  CS .
C
      DOUBLE PRECISION FP,C,FPS,CS
      DIMENSION FP(N),C(N,N),FPS(N),CS(N,N)
C
      DO 20 I = 1,N
         FPS(I) = FP(I)
         DO 10 J = 1,N
            CS(I,J) = C(I,J)
  10        CONTINUE
  20     CONTINUE
      RETURN
      END
      SUBROUTINE NEWF (N,X,WF,PRF,IDER,NU,NO,IDA,
     1    XTP,VFU,C,FP,BUF,NRP,NJEV,IER)
C
C THE SUBROUTINE NEWF COMPUTES THE NEW WEIGHTS FOR FUNCTION RESCALING
C THE PARAMETERS ARE DESCRIBED WITHIN THE SUBROUTINE  DAFNE .
C
      DOUBLE PRECISION X,VFU,FP,BUF,C,WF,PRF,XTP
      DOUBLE PRECISION RWMAX,RWMIN
      DOUBLE PRECISION BMAX,BMIN
      DOUBLE PRECISION SB,B,RW
      DIMENSION X(N),VFU(N),FP(N),BUF(N),C(N,N),WF(N),XTP(N)
C
C RWMAX  AND  RWMIN  ARE THE MAX AND MIN VALUES OF THE
C MULTIPLICATIVE FACTOR USED TO UPDATE THE WEIGHTS
      DATA RWMAX,RWMIN/1.D+3,1.D-3/
C
C BMAX  AND  BMIN  DEFINE THE ALLOWED RANGE FOR THE SQUARED NORM  B .
C USE THE FIRST VALUES (40,-40) OR THE SECOND ONES (30,-30) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1.D120 (E.G. UNIVAC 1100), OR BELOW
C 1.D120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA BMAX,BMIN/1.D+40,1.D-40/
C     DATA BMAX,BMIN/1.D+30,1.D-30/
C
C COMPUTE FUNCTION VALUES
      CALL VECF (N,X,VFU,NRP,IER)
      IF (IER.LT.0) RETURN
C
C COMPUTE JACOBIAN
      CALL JCBN (N,X,VFU,C,FP,BUF,XTP,PRF,IDER,NU,NO,IDA,
     1    NRP,NJEV,IER)
C
      SB = 1.D0
      DO 20 I = 1,N
C
C        COMPUTE SQUARED NORM OF GRADIENT OF FUNCTION  I
         B = 0.D0
         DO 10 J = 1,N
            B = B+C(I,J)**2
  10        CONTINUE
C
         B = DMAX1(B,BMIN)
         B = DMIN1(B,BMAX)
         BUF(I) = WF(I)*B
         SB = SB*BUF(I)**(1.D0/DBLE(FLOAT(N)))
C
  20     CONTINUE
C
C COMPUTE WEIGHTS
      DO 30 I = 1,N
         RW = SB/BUF(I)
         RW = DMIN1(RW,RWMAX)
         RW = DMAX1(RW,RWMIN)
         WF(I) = WF(I)*RW
  30     CONTINUE
C
      RETURN
      END
      SUBROUTINE NEWXTP (N,X,VFU,C,XTP)
C
C THE SUBROUTINE  NEWXTP  UPDATES THE ARRAY  XTP
C USED BY SUBROUTINE  JCBN  TO COMPUTE THE STEPLENGTHS USED
C FOR COMPUTING THE FINITE-DIFFERENCES JACOBIAN.
C
C THE PARAMETERS ARE DEFINED IN THE SUBROUTINE  EVOLV .
C
      DOUBLE PRECISION C,ECIJ,XTP,RM,RP,VFI,VFU,X
      DIMENSION X(N),VFU(N),C(N,N),XTP(N)
C
C RP  AND  RM  ARE THE VALUES OF THE MULTIPLICATIVE
C     FACTOR USED TO UPDATE THE VECTOR  XTP .
      DATA RP,RM/10.D0,.1D0/
C
      DO 20 I = 1,N
         IM = 0
         VFI = DABS(VFU(I))
         DO 10 J = 1,N
            ECIJ = DABS(C(I,J)*XTP(J))
            IC = 1
            IF (VFI.LT.ECIJ*RP) IC = 2
            IF (VFI.LT.ECIJ*RM) IC = 3
            IM = MAX0(IM,IC)
  10        CONTINUE
         IF (IM.EQ.1) GO TO 50
         IF (IM.EQ.2) XTP(I) = -DABS(XTP(I))
         IF (IM.EQ.3) XTP(I) = DABS(XTP(I))
  20     CONTINUE
      DO 40 J = 1,N
         IM = 0
         DO 30 I = 1,N
            VFI = DABS(VFU(I))
            ECIJ = DABS(C(I,J)*XTP(J))
            IC = 1
            IF (VFI.LT.ECIJ*RP) IC = 2
            IF (VFI.LT.ECIJ*RM) IC = 3
            IF (IC.EQ.2.AND.XTP(I).LT.0.D0) IC = 4
            IM = MAX0(IM,IC)
  30        CONTINUE
         IF (IM.EQ.1) XTP(J) = XTP(J)*RP
         IF (IM.EQ.3) XTP(J) = XTP(J)*RM
  40     CONTINUE
      GO TO 70
  50  CONTINUE
      DO 60 J = 1,N
         XTP(J) = XTP(J)*RP
  60     CONTINUE
  70  CONTINUE
      DO 80 J = 1,N
         XTP(J) = DMAX1(DABS(XTP(J)),DABS(X(J)))
  80     CONTINUE
      RETURN
      END
      SUBROUTINE PRMACH (PRM)
C
C THE SUBROUTINE  PRMACH  ESTIMATES THE MACHINE PRECISION  PRM
C (SOFTWARE PRECISION, ESTIMATED ON SOME INTRINSIC
C FORTRAN FUNCTIONS).
C
      DOUBLE PRECISION A,E,P,PRM
C
      P = 0.D0
      DO 10 I = 1,64
         A = DBLE(FLOAT(I))
         E = 0.D0
         E = E+(DSQRT(A*A)-A)**2
         E = E+(DSQRT(A)**2-A)**2
         E = E+(DEXP(DLOG(A))-A)**2
         E = E+(DLOG(DEXP(A))-A)**2
         E = E+(A*DSIN(A)**2+A*DCOS(A)**2-A)**2
         P = P+E/(5.D0*A**2)
  10     CONTINUE
      P = DSQRT(P/64.D0)
      PRM = P
C
      RETURN
      END
      SUBROUTINE QMS (N,FP,C,FPS,CS,D)
C
C THE SUBROUTINE  QMS  REGULARIZES THE LINEAR SYSTEM
C WITH  N BY N  COEFFICIENT MATRIX  C  AND RIGHT-HAND-SIDE
C N-VECTOR  FP  WITH REGULARIZATION PARAMETER  D
C OBTAINING THE LINEAR SYSTEM  CS , FPS .
C
      DOUBLE PRECISION FP,C,FPS,CS,D
      DIMENSION FP(N),C(N,N),FPS(N),CS(N,N)
C
      DO 30 I = 1,N
         FPS(I) = 0.D0
         DO 20 J = 1,N
            FPS(I) = FPS(I)+C(J,I)*FP(J)
            CS(I,J) = 0.D0
            DO 10 K = 1,N
               CS(I,J) = CS(I,J)+C(K,I)*C(K,J)
  10           CONTINUE
  20        CONTINUE
         CS(I,I) = CS(I,I)+D
  30     CONTINUE
C
      RETURN
      END
C MAIN (TEST VERSION).
C
C THE PRESENT VERSION OF THE MAIN PROGRAM IS DESIGNED TO TEST
C THE  DAFNE  ALGORITHM, AS IMPLEMENTED IN THE PROPOSED SUBROUTINES,
C ON THE TEST PROBLEMS DEFINED IN THE SUBROUTINE  VECFCN .
C BY MEANS OF A SUITABLE DECK OF DATA CARDS THE USER CAN
C SETUP A RUN CONTAINING AN ARBITRARILY LONG SEQUENCE
C OF TESTS ON THE ABOVE PROBLEMS, WITH THE POSSIBILITY
C OF RESCALING THE VARIABLES, THE FUNCTIONS, AND THE INITIAL POINT,
C AS SUGGESTED BY MORE, GARBOW, AND HILLSTROM,
C AND OF CHANGING THE VALUES OF A NUMBER OF
C NUMERICAL OR CONTROL VARIABLES.
C FOR EACH TEST THE DATA CARDS MUST BE
C   A CARD WITH PROBLEM NUMBER AND TITLE
C   AN ARBITRARY NUMBER OF OPTIONAL DATA CARDS FOR CHANGING
C     THE VALUES OF A NUMBER OF VARIABLES, AS EXPLAINED BELOW
C   A LAST (BLANK) DATA CARD, TO START TEST EXECUTION.
C ANOTHER BLANK CARD IS REQUIRED TO INDICATE THE END OF THE TESTS.
C
      DOUBLE PRECISION FAC,FAC0,FMIN
      DOUBLE PRECISION PRF,PRM,RFACTR,RVAL
      DOUBLE PRECISION X,XX
      DOUBLE PRECISION B1,B2,B3,B4,B5,B6,B7,B8,C1,C2
C
C COMMON/MORE/    USED BY SUBROUTINES  JCBNA, VECF .
      COMMON/MORE/NPROB
C COMMON /CHEM/   USED BY SUBROUTINE VECFCN .
      COMMON /CHEM/RFACTR
C COMMON/RISC/    USED BY SUBROUTINES JCBNA, VECF .
      COMMON/RISC/ISCL
C
      DIMENSION IA(70)
      DIMENSION FAC0(10)
C
      DIMENSION X(100),XX(100)
      DIMENSION FMIN(100)
      DIMENSION B1(100),B2(100),B3(100),B4(100),B5(100),B6(100)
      DIMENSION B7(100),B8(100),C1(100,100),C2(100,100)
C
C NMAX IS THE DIMENSION OF PRECEDING ARRAYS.
      DATA NMAX/100/
C
      DATA IB/1H /
C
C ESTIMATE MACHINE PRECISION  PRM
      CALL PRMACH (PRM)
      WRITE (6,370) PRM
C
  10  CONTINUE
C
C SET DEFAULT VALUES
C
C   NUMBER OF EQUATIONS
      N = 1
C   INITIAL POINT
      XX(1) = 0.D0
C   THRESHOLD VALUE FOR THE FUNCTION VALUES
      DO 20 I = 1,NMAX
         FMIN(I) = 1.D-10
  20     CONTINUE
C   CONTROL VARIABLE FOR PRINTING (IF ICP.NE.0) ANY MODIFICATION
C     TO THE DEFAULT VALUES
      ICP = 0
C   CONTROL VARIABLE USED BY SUBROUTINE  JCBN
C     IDER = 0    ANALYTICAL DERIVATIVES
C     IDER = 1    NUMERICAL FORWARD DIFFERENCES DERIVATIVES
C     IDER = 2    NUMERICAL CENTRAL DIFFERENCES DERIVATIVES
      IDER = 1
C   CONTROL VARIABLE FOR PRINTING DATA AND RESULTS OF EACH TEST
C   IN EXTENDED FORM IF  ITAB = 0 , AND IN TABLE FORM IF  ITAB = 1 .
      ITAB = 0
C   NUMBER OF RUNS STARTING FROM INITIAL POINTS DIFFERING BY A
C     MULTIPLICATION FACTOR  FAC0
      NRIP = 3
C   CONTROL VARIABLE FOR THE CHOICE OF THE INITIAL POINT
C     INIT = 1    CHOICE AND RESCALING OF THE INITIAL POINT
C                 PERFORMED BY THE SUBROUTINE  INITPT .
C     INIT.NE.1   INITIAL POINT READ FROM DATA CARD, AFTER READING  N ,
C                 RESCALING PERFORMED BY THIS MAIN PROGRAM.
      INIT = 1
C   MULTIPLICATION FACTORS FOR CHANGING THE INITIAL POINT
      FAC0(1) = 1.D0
      FAC0(2) = 10.D0
      FAC0(3) = 100.D0
C   ESTIMATED RELATIVE PRECISION OF FUNCTION VALUES
      PRF = 1.D-14
C
C   NUMBERS DEFINING THE (POSSIBLY) BANDED STRUCTURE OF
C   THE JACOBIAN MATRIX
C   NO IS THE N. OF NON-ZERO JACOBIAN DIAGONALS ABOVE MAIN DIAGONAL
C   NU IS THE N. OF NON-ZERO JACOBIAN DIAGONALS UNDER MAIN DIAGONAL
C   (FOR A NON-BANDED MATRIX  NU = NO = N-1 )
      NO = N-1
      NU = N-1
C
C   NRS1  AND  NRS2  ARE THE EXTREME VALUES FOR  ISCL , I.E. THE
C   PROBLEMS ARE RUN WITH  ISCL = NRS1,...,NRS2  (STEP 1), WHERE
C   ISCL  IS A CONTROL VARIABLE USED BY SUBROUTINE  VECF
C     IF ISCL.EQ.1 THE PROBLEMS ARE RUN WITHOUT RESCALING
C     IF ISCL.EQ.2 THE PROBLEMS ARE RUN WITH RESCALED VARIABLES
C     IF ISCL.EQ.3 THE PROBLEMS ARE RUN WITH RESCALED FUNCTIONS
C
      NRS1 = 1
      NRS2 = 1
C
C   FACTOR USED BY VECFCN (CHEMICAL EQUILIBRIUM PROBLEM 3)
      RFACTR = 1.D0
C
  30  CONTINUE
C
C READ TEST PROBLEM DATA
C
C READ PROBLEM NUMBER AND NAME
      READ (5,380) NPROB,IA
C
      IF (ICP.NE.0) WRITE (6,390)
C
C TEST FOR END-OF-TESTS CARD
      IF (NPROB.LE.0) GO TO 360
C
  40  CONTINUE
C
C POSSIBLE CHANGES IN NUMERICAL VALUES
C   IVAR.LT.0         RESET DEFAULT VALUES
C   IVAR.EQ.0         END OF CHANGES
C   1.LE.IVAR.LE.18   THE VARIABLE NUMBERED  IVAR
C                     IS TO BE GIVEN A NEW VALUE
C   IVAR.GT.18        DATA CARD IGNORED
C
C VARIABLES THAT CAN BE GIVEN NEW VALUES, AND CORRESPONDING
C IDENTIFICATION NUMBER  IVAR
C
C N       FMIN IDER    NRIP NRS1 NRS2 FAC0 ITAB    INIT PRF   NU NO ICP
C 1 2  3   4    5    6   7    8    9   10   11  12  13  14 15 16 17 18
C
C IVAL       NEW VALUE FOR VARIABLE NUMBERED  IVAR , IF INTEGER
C              (ELSE  IVAL=0 )
C RVAL       NEW VALUE FOR VARIABLE NUMBERED  IVAR , IF REAL
C              (ELSE  RVAL=0. )
C
C READ (AND PRINT) DATA MODIFYING THE DEFAULT VALUES
      READ (5,400) IVAR,IVAL,RVAL
      IF (ICP.NE.0.AND.IVAR.NE.0) WRITE (6,410) IVAR,IVAL,RVAL
C
      IF (IVAR) 10,240,50
  50  CONTINUE
      IF (IVAR.GE.19) GO TO 40
      GO TO (60,70,80,90,100,110,120,130,140,150,160,170,180,190,200,210
     1    ,220,230), IVAR
  60  CONTINUE
      N = IVAL
      NU = N-1
      NO = NU
      IF (INIT.NE.1) READ (5,420) (XX(I),I = 1,N)
      GO TO 40
  70  CONTINUE
      GO TO 40
  80  CONTINUE
      GO TO 40
  90  CONTINUE
      READ (5,420) (FMIN(I),I = 1,N)
      GO TO 40
 100  IDER = IVAL
      GO TO 40
 110  CONTINUE
      GO TO 40
 120  NRIP = IVAL
      GO TO 40
 130  CONTINUE
      NRS1 = IVAL
      GO TO 40
 140  CONTINUE
      NRS2 = IVAL
      GO TO 40
 150  READ (5,420) (FAC0(I),I = 1,NRIP)
      GO TO 40
 160  CONTINUE
      ITAB = IVAL
      IF (ITAB.EQ.1) WRITE (6,390)
      IF (ITAB.EQ.1) WRITE (6,500)
      GO TO 40
 170  CONTINUE
      GO TO 40
 180  INIT = IVAL
      GO TO 40
 190  CONTINUE
      PRF = RVAL
      GO TO 40
 200  CONTINUE
      GO TO 40
 210  NU = IVAL
      GO TO 40
 220  NO = IVAL
      GO TO 40
 230  ICP = IVAL
      GO TO 40
 240  CONTINUE
C
C START RESCALING LOOP
      DO 350 JRIS = NRS1,NRS2
         ISCL = JRIS
C
C        START INITIAL-POINT LOOP
         DO 340 INPT = 1,NRIP
            FAC = FAC0(INPT)
C
C           CHANGE INITIAL POINT
            IF (INIT.EQ.1) GO TO 280
            IX0 = 1
            DO 250 K = 1,N
               IF (XX(K).EQ.0.D0) GO TO 250
               IX0 = 0
               GO TO 260
 250           CONTINUE
 260        CONTINUE
            DO 270 K = 1,N
               IF (IX0.EQ.0) X(K) = XX(K)*FAC
               IF (IX0.EQ.1.AND.FAC.EQ.1.D0) X(K) = XX(K)
               IF (IX0.EQ.1.AND.FAC.NE.1.D0) X(K) = FAC
 270           CONTINUE
            GO TO 290
 280        CONTINUE
            CALL INITPT (N,X,NPROB,FAC)
 290        CONTINUE
C
C           REVERSE SCALING TO PRESERVE TRUE INITIAL POINT
            IF (ISCL.NE.2.OR.N.EQ.1) GO TO 310
            DO 300 K = 1,N
               X(K) = X(K)*.1D0**(DBLE(FLOAT(5*(K+K-N-1)))/
     1             DBLE(FLOAT(N-1)))
 300           CONTINUE
 310        CONTINUE
C
C           PRINT TEST PROBLEM DATA
            IF (ITAB.EQ.1) GO TO 320
            WRITE (6,430) IA,NPROB,N
            WRITE (6,520) ISCL
            WRITE (6,530) FAC
            WRITE (6,540)
            WRITE (6,550) (IB,I,X(I),I = 1,N)
            WRITE (6,440)
            WRITE (6,450) (IB,I,FMIN(I),I = 1,N)
            WRITE (6,460)
            WRITE (6,470) IDER,PRF,NU,NO
 320        CONTINUE
C
            CALL DAFNE (N,X,FMIN,PRF,IDER,NU,NO,IDA,
     1          NFEV,NJEV,NSAC,NST,IER,
     2          B1,B2,B3,B4,B5,B6,B7,B8,C1,C2)
C
C           RUN TERMINATION MESSAGES
C
C           INDICATOR OF DERIVATIVES MODE
C           IDERM = 0  ANALYTICAL DERIVATIVES           (IDER=0, IDA=0)
C           IDERM = 1  FORWARD-DIFFERENCES DERIVATIVES  (IDER=1, IDA=1)
C           IDERM = 2  CENTRAL-DIFFERENCES DERIVATIVES  (IDER=2, IDA=1)
C           IDERM =-1  REQUESTED ANALYTICAL DERIVATIVES
C                      NOT AVAILABLE, FORWARD-DIFFERENCES
C                      DERIVATIVES USED                 (IDER=0, IDA=1)
C
            IDERM = IDER
            IF (IDER.EQ.0.AND.IDA.EQ.0) IDERM = -1
            IF (ITAB.EQ.1) GO TO 330
            WRITE (6,480)
            WRITE (6,490) IER
            WRITE (6,560)
            WRITE (6,550) (IB,I,X(I),I = 1,N)
            WRITE (6,570)
            WRITE (6,580) (IB,I,B3(I),I = 1,N)
            WRITE (6,590) IDERM
            WRITE (6,600) NFEV,NJEV,NSAC,NST
 330        CONTINUE
            IF (ITAB.EQ.1) WRITE (6,510) NPROB,N,ISCL,INPT,
     1          IDERM,IER,NSAC,NST,NFEV,NJEV
C
 340        CONTINUE
C        END INITIAL-POINT LOOP
C
 350     CONTINUE
C END RESCALING LOOP
C
C TAKE NEW TEST PROBLEM
      GO TO 30
C
 360  CONTINUE
C END OF TEST PROBLEMS
C
      WRITE (6,390)
C
      STOP
C
 370  FORMAT (1H0,36H ESTIMATED MACHINE PRECISION    PRM=,G20.10/)
 380  FORMAT (I10,70A1)
 390  FORMAT (1H1)
 400  FORMAT (2I5,G20.10)
 410  FORMAT (10X,18H VARIABLE NUMBERED,I5/10X,10H NEW VALUE,I10,
     1    11H IF INTEGER/20X,G20.10,8H IF REAL/)
 420  FORMAT (4G15.10)
 430  FORMAT (1H1,10X,70A1//10X,36H PROBLEM NUMBER               NPROB=,
     1    I4/10X,36H NUMBER OF EQUATIONS              N=,I4/)
 440  FORMAT (1H0,9X,41H THRESHOLD VALUES FOR THE FUNCTION VALUES)
 450  FORMAT (5(1A1,5HFMIN(,I3,2H)=,G14.7))
 460  FORMAT (1H0,9X,42H VARIABLES USED FOR COMPUTING THE JACOBIAN)
 470  FORMAT (10X,6H IDER=,I2/10X,6H PRF =,G14.7/10X,5H NU =,I4/10X,
     1    5H NO =,I4/)
 480  FORMAT (1H0,19X,18H NUMERICAL RESULTS)
 490  FORMAT (1H0,9X,31H ERROR INDICATOR          IER =,I2)
 500  FORMAT (1H0,1H*,7H  NPROB,7H   N   ,7H  ISCL ,7H  INPT ,7H  IDERM,
     1    7H  IER  ,7H  NSAC ,7H  NST  ,7H  NFEV ,7H  NJEV ,1H*)
 510  FORMAT (1H0,1H*,10(I5,2X),1H*)
 520  FORMAT (1H0,9X,40H CONTROL VARIABLE FOR PROBLEM RESCALING ,
     1    10H    ISCL =,I2)
 530  FORMAT (1H0,9X,40H SCALING FACTOR FOR THE INITIAL POINT   ,
     1    10H     FAC =,G14.7)
 540  FORMAT (1H0,9X,14H INITIAL POINT)
 550  FORMAT (5(1A1,5H   X(,I3,2H)=,G14.7))
 560  FORMAT (1H0,9X,12H FINAL POINT)
 570  FORMAT (1H0,9X,22H FINAL FUNCTION VALUES)
 580  FORMAT (5(1A1,5H B3 (,I3,2H)=,G14.7))
 590  FORMAT (1H0,9X,42H INDICATOR OF DERIVATIVES MODE     IDERM =,I2)
 600  FORMAT (1H0,9X,12H COUNTERS OF
     1    /12X,45H FUNCTION EVALUATIONS                  NFEV =,I4
     2    /12X,45H JACOBIAN EVALUATIONS                  NJEV =,I4
     3    /12X,45H ACCEPTED TIME INTEGRATION STEPS       NSAC =,I4
     4    /12X,45H ATTEMPTED TIME INTEGRATION STEPS       NST =,I4)
      END
      SUBROUTINE VECF (N,X,VFU,NR,IER)
C
C SUBROUTINE  VECF (N,X,VFU,IER)
C THE SUBROUTINE  VECF  RETURNS IN THE N-VECTOR  VFU
C THE VALUES OF THE FUNCTIONS OF SYSTEM (1) (SEE SUBROUTINE  EVOLV )
C COMPUTED AT THE POINT  X .
C NR  IS THE NUMBER OF TIME INTEGRATION STEPS IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE INDICATES
C     THAT  X  IS OUTSIDE A SAFETY REGION DEFINED BY THE USER).
C
C TEST VERSION
C THE PRESENT VERSION OF THE SUBROUTINE  VECF  COMPUTES
C THE FUNCTION VALUES FOR THE PROBLEM NUMBERED
C NPROB  IN THE LIST OF THE SELECTED TEST PROBLEMS
C BY CALLING THE SUBROUTINE  VECFCN .
C IN ORDER TO AVOID FAILURES (E.G. OVERFLOWS) INSIDE THE SUBROUTINE
C VECFCN  THE PRESENT VERSION OF  VECF  CALLS THE SUBROUTINE  TRANGE
C WHICH TESTS IF THE FUNCTION ARGUMENTS ARE WITHIN A SAFETY
C REGION .
C
      DOUBLE PRECISION VFU,X
      DOUBLE PRECISION XM
      DIMENSION X(N),VFU(N)
      DIMENSION XM(100)
      COMMON/MORE/NPROB
      COMMON/RISC/ISCL
C
C SCALE VARIABLES
      IF (ISCL.NE.2.OR.N.EQ.1) GO TO 20
      DO 10 K = 1,N
         XM(K) = X(K)
         X(K) = X(K)*10.D0**(DBLE(FLOAT(5*(K+K-N-1)))/DBLE(FLOAT(N-1)))
  10     CONTINUE
  20  CONTINUE
C
C CALL PROTECTION SUBROUTINE
      CALL TRANGE (N,X,NPROB,IER)
C
      IF (IER.LT.0) GO TO 100
C
C CALL FUNCTION COMPUTING SUBROUTINE
      CALL VECFCN (N,X,VFU,NPROB)
C
      IF (N.EQ.1) GO TO 90
      GO TO (80,30,50), ISCL
C
C RESTORE NON-SCALED VARIABLES
  30  CONTINUE
      DO 40 K = 1,N
         X(K) = XM(K)
  40     CONTINUE
      GO TO 70
C
C SCALE FUNCTION VALUES
  50  CONTINUE
      DO 60 K = 1,N
         VFU(K) = VFU(K)*10.D0**(DBLE(FLOAT(5*(K+K-N-1)))/
     1       DBLE(FLOAT(N-1)))
  60     CONTINUE
  70  CONTINUE
  80  CONTINUE
  90  CONTINUE
      RETURN
 100  CONTINUE
      RETURN
C
      END
      SUBROUTINE JCBNA (N,X,C,NR,IDA,IER)
C
C THE SUBROUTINE  JCBNA  COMPUTES THE ANALYTICAL JACOBIAN (IF AVAILABLE)
C THE PARAMETERS ARE DESCRIBED WITHIN THE SUBROUTINE  JCBN
C
C TEST VERSION
C THE PRESENT VERSION PROVIDES THE ANALYTICAL JACOBIAN
C FOR SOME OF THE SELECTED TEST PROBLEMS DEFINED IN
C THE SUBROUTINE  VECFCN  AND COMPUTES THE ANALYTICAL
C JACOBIAN FOR THE PROBLEM NUMBERED  NPROB  IN THE
C TEST PROBLEM LIST.
C IF THE ANALYTICAL JACOBIAN OF PROBLEM  NPROB  IS NOT CODED,  JCBNA
C RETURNS  IDA = 0, OTHERWISE  IDA = 1 .
C IN ORDER TO AVOID FAILURES (E.G. OVERFLOWS) INSIDE THE SUBROUTINE
C VECFCN  THE PRESENT VERSION OF  JCBNA  CALLS THE SUBROUTINE  TRANGE
C WHICH TESTS IF THE FUNCTION ARGUMENTS ARE WITHIN A SAFETY
C REGION .
C
      DOUBLE PRECISION C,DSQRT5,DSQR10,SJ,X,XM
      DIMENSION X(N),C(N,N)
      DIMENSION XM(100)
      COMMON/RISC/ISCL
      COMMON/MORE/NPROB
      IF (ISCL.NE.2.OR.N.EQ.1) GO TO 20
      DO 10 K = 1,N
         XM(K) = X(K)
         X(K) = X(K)*10.D0**(DBLE(FLOAT(5*(K+K-N-1)))/DBLE(FLOAT(N-1)))
  10     CONTINUE
  20  CONTINUE
      IDA = 0
C
      CALL TRANGE (N,X,NPROB,IER)
C
      IF (IER.LT.0) RETURN
      GO TO (30,40,60,70,80,90,100,110,160,170,180,220,230,240,250,260,
     1    270), NPROB
C
C ROSENBROCK FUNCTION
  30  CONTINUE
      IDA = 1
      C(1,1) = -1.D0
      C(1,2) = 0.D0
      C(2,1) = -20.D0*X(1)
      C(2,2) = 10.D0
      GO TO 280
C
C POWELL SINGULAR FUNCTION
  40  CONTINUE
      IDA = 1
      IF (NR.GT.1) GO TO 50
      DSQRT5 = DSQRT(5.D0)
      DSQR10 = DSQRT(10.D0)
  50  CONTINUE
      C(1,1) = 1.D0
      C(1,2) = 10.D0
      C(1,3) = 0.D0
      C(1,4) = 0.D0
      C(2,1) = 0.D0
      C(2,2) = 0.D0
      C(2,3) = DSQRT5
      C(2,4) = -C(2,3)
      C(3,1) = 0.D0
      C(3,2) = 2.D0*X(2)-4.D0*X(3)
      C(3,3) = 8.D0*X(3)-4.D0*X(2)
      C(3,4) = 0.D0
      C(4,1) = DSQR10*(X(1)-X(4))
      C(4,2) = 0.D0
      C(4,3) = 0.D0
      C(4,4) = -C(4,1)
      GO TO 280
C
C POWELL BADLY SCALED FUNCTION
  60  CONTINUE
      IDA = 1
      C(1,1) = 1.D4*X(2)
      C(1,2) = 1.D4*X(1)
      C(2,1) = -DEXP(-X(1))
      C(2,2) = -DEXP(-X(2))
      GO TO 280
C
  70  CONTINUE
      GO TO 280
C
  80  CONTINUE
      GO TO 280
C
  90  CONTINUE
      GO TO 280
C
 100  CONTINUE
      GO TO 280
C
C BROWN ALMOST LINEAR FUNCTION
 110  CONTINUE
      IDA = 1
      DO 130 I = 1,N
         DO 120 J = 1,N
            C(I,J) = 1.D0
 120        CONTINUE
         C(I,I) = 2.D0
 130     CONTINUE
      DO 150 J = 1,N
         C(N,J) = 1.D0
         DO 140 I = 1,N
            IF (I.NE.J) C(N,J) = C(N,J)*X(I)
 140        CONTINUE
 150     CONTINUE
      GO TO 280
C
 160  CONTINUE
      GO TO 280
C
 170  CONTINUE
      GO TO 280
C
C TRIGONOMETRIC FUNCTION
 180  CONTINUE
      IDA = 1
      DO 200 J = 1,N
         SJ = DSIN(X(J))
         DO 190 I = 1,N
            C(I,J) = SJ
 190        CONTINUE
 200     CONTINUE
      DO 210 I = 1,N
         C(I,I) = -DCOS(X(I))+DBLE(FLOAT(I))*DSIN(X(I))+C(I,I)
 210     CONTINUE
      GO TO 280
C
 220  CONTINUE
      GO TO 280
C
 230  CONTINUE
      GO TO 280
C
 240  CONTINUE
      GO TO 280
C
C CHEMICAL EQUILIBRIUM PROBLEM 1
 250  CONTINUE
      IDA = 1
      C(1,1) = 0.D0
      C(1,2) = 1.D0
      C(2,1) = 1.D0
      C(2,2) = -1.D0
      GO TO 280
C
C CHEMICAL EQUILIBRIUM PROBLEM 2
 260  CONTINUE
      IDA = 1
      C(1,1) = 1.D0
      C(1,2) = 1.D0
      C(1,3) = 0.D0
      C(1,4) = 1.D0
      C(1,5) = 0.D0
      C(1,6) = 0.D0
      C(2,1) = 0.D0
      C(2,2) = 0.D0
      C(2,3) = 0.D0
      C(2,4) = 0.D0
      C(2,5) = 1.D0
      C(2,6) = 1.D0
      C(3,1) = 1.D0
      C(3,2) = 1.D0
      C(3,3) = 1.D0
      C(3,4) = 0.D0
      C(3,5) = 2.D0
      C(3,6) = 1.D0
      C(4,1) = 1.D0
      C(4,2) = -.1D0
      C(4,3) = 0.D0
      C(4,4) = 0.D0
      C(4,5) = 0.D0
      C(4,6) = 0.D0
      C(5,1) = 1.D0
      C(5,2) = 0.D0
      C(5,3) = -1.D4*X(4)
      C(5,4) = -1.D4*X(3)
      C(5,5) = 0.D0
      C(5,6) = 0.D0
      C(6,1) = 0.D0
      C(6,2) = 0.D0
      C(6,3) = -55.D14*X(6)
      C(6,4) = 0.D0
      C(6,5) = 1.D0
      C(6,6) = -55.D14*X(3)
      GO TO 280
C
 270  CONTINUE
      GO TO 280
C
C RESCALE
 280  CONTINUE
      IF (IDA.EQ.0.OR.N.EQ.1) RETURN
      GO TO (360,290,320), ISCL
C
C RESCALE COLUMNS OF JACOBIAN
 290  CONTINUE
      DO 310 K = 1,N
         X(K) = XM(K)
         DO 300 J = 1,N
            C(J,K) = C(J,K)*10.D0**(DBLE(FLOAT(5*(K+K-N-1)))/
     1          DBLE(FLOAT(N-1)))
 300        CONTINUE
 310     CONTINUE
      GO TO 350
C
C RESCALE ROWS OF JACOBIAN
 320  CONTINUE
      DO 340 K = 1,N
         DO 330 J = 1,N
            C(K,J) = C(K,J)*10.D0**(DBLE(FLOAT(5*(K+K-N-1)))/
     1          DBLE(FLOAT(N-1)))
 330        CONTINUE
 340     CONTINUE
C
 350  CONTINUE
 360  CONTINUE
      RETURN
      END
      SUBROUTINE INITPT(N,X,NPROB,FACTOR)
C
C THE SUBROUTINE  INITPT  SETS THE STANDARD INITIAL
C POINTS FOR THE TEST FUNCTIONS DEFINED IN THE SUBROUTINE  VECFN .
C FROM MINPACK OF J. MORE
C
      DOUBLE PRECISION C1,FACTOR,H,HALF
      DOUBLE PRECISION ONE,THREE,TJ,X,ZERO
      DIMENSION X(N)
      DATA ZERO,HALF,ONE,THREE,C1 /0.0D0,.5D0,1.D0,3.D0,1.2D0/
      GO TO (10,20,30,40,50,60,80,100,120,120,140,160,180,180,
     1    300,310,320,320,320,320),NPROB
C
C ROSENBROCK
C
   10 CONTINUE
      X(1) = -C1
      X(2) = ONE
      GO TO 200
C
C POWELL SINGULAR
C
   20 CONTINUE
      X(1) = THREE
      X(2) = -ONE
      X(3) = ZERO
      X(4) = ONE
      GO TO 200
C
C POWELL BADLY SCALED
C
   30 CONTINUE
      X(1) = ZERO
      X(2) = ONE
      GO TO 200
C
C WOOD
C
   40 CONTINUE
      X(1) = -THREE
      X(2) = -ONE
      X(3) = -THREE
      X(4) = -ONE
      GO TO 200
C
C HELICAL VALLEY
C
   50 CONTINUE
      X(1) = -ONE
      X(2) = ZERO
      X(3) = ZERO
      GO TO 200
C
C WATSON
C
   60 CONTINUE
      DO 70 J = 1,N
         X(J) = ZERO
   70    CONTINUE
      GO TO 200
C
C CHEBYQUAD
C
   80 CONTINUE
      H = ONE/DBLE(FLOAT(N+1))
      DO 90 J = 1,N
         X(J) = DBLE(FLOAT(J))*H
   90    CONTINUE
      GO TO 200
C
C BROWN
C
  100 CONTINUE
      DO 110 J = 1,N
         X(J) = HALF
  110    CONTINUE
      GO TO 200
C
C DISCRETE B.V.P.   AND INTEGRAL EQUAT.
C
  120 CONTINUE
      H = ONE/DBLE(FLOAT(N+1))
      DO 130 J = 1,N
         TJ = DBLE(FLOAT(J))*H
         X(J) = TJ*(TJ-ONE)
  130    CONTINUE
      GO TO 200
C
C TRIGONOMETRIC
C
  140 CONTINUE
      H = ONE/DBLE(FLOAT(N))
      DO 150 J = 1,N
         X(J) = H
  150    CONTINUE
      GO TO 200
C
C VARIABLY DIMENSIONED FUNCTIONS
C
  160 CONTINUE
      H = ONE/DBLE(FLOAT(N))
      DO 170 J = 1,N
         X(J) = ONE-DBLE(FLOAT(J))*H
  170    CONTINUE
      GO TO 200
C
C BROYDEN TRIDIAGONAL AND BANDED
C
  180 CONTINUE
      DO 190 J = 1,N
         X(J) = -ONE
  190    CONTINUE
      GO TO 200
C
C CHEMICAL 1
C
  300 CONTINUE
      X(1) = ZERO
      X(2) = ZERO
      GO TO 200
C
C CHEMICAL 2
C
  310 CONTINUE
      DO 315 J = 1,6
         X(J) = ZERO
  315    CONTINUE
      GO TO 200
C
C CHEMICAL 3
C
  320 CONTINUE
      DO 325 J = 1,10
         X(J) = ZERO
  325    CONTINUE
C
  200 CONTINUE
C
C COMPUTE MULTIPLE OF  INITIAL POINT
C
      IF(FACTOR.EQ.ONE)GO TO 250
      IF(NPROB.EQ.6.OR.NPROB.GT.14)GO TO 220
      DO 210 J = 1,N
         X(J) = FACTOR*X(J)
  210    CONTINUE
      GO TO 240
  220 CONTINUE
      DO 230 J = 1,N
         X(J) = FACTOR
  230    CONTINUE
  240 CONTINUE
  250 CONTINUE
      RETURN
      END
      SUBROUTINE VECFCN(N,X,FVEC,NPROB)
      INTEGER N,NPROB
      COMMON/CHEM/RFACTR
      DOUBLE PRECISION X(N),FVEC(N)
C     **********
C
C     SUBROUTINE VECFCN
C     THIS SUBROUTINE DEFINES FOURTEEN TEST FUNCTIONS. THE FIRST
C     FIVE TEST FUNCTIONS ARE OF DIMENSIONS 2,4,2,4,3, RESPECTIVELY,
C     WHILE THE REMAINING TEST FUNCTIONS ARE
C     OF VARIABLE DIMENSION N FOR ANY N GREATER THAN OR EQUAL
C     TO 1 (PROBLEMS 6 IS AN EXCEPTION TO THIS, SINCE IT
C     DOES NOT ALLOW N = 1).
C
C     THE SUBROUTINE STATEMENT IS
C       SUBROUTINE VECFCN(N,X,FVEC,NPROB)
C     WHERE
C       N IS A POSITIVE INTEGER VARIABLE.
C
C       X IS A LINEAR ARRAY OF LENGTH N.
C
C       FVEC IS A LINEAR ARRAY OF LENGTH N,ON OUTPUT FVEC
C       CONTAINS THE NPROB FUNCTION VECTOR EVALUATED AT X.
C
C       NPROB IS A POSITIVE INTEGER VARIABLE WHICH DEFINES THE
C         NUMBER OF THE PROBLEM. NPROB MUST NOT EXCEED 14.
C
C     FORTRAN-SUPPLIED ... DATAN,DCOS,DEXP,DSIGN,DSIN,DSQRT
C                          MAX0,MIN0
C
C     MINPACK. VERSION OF JULY 1978.
C     BURTON S. GARBOW, KENNETH E. HILLSTROM, JORGE J. MORE
C     INTEGER I,IVAR,J,K,K1,K2,KP1,NM1,NMJP1,ML,MU
C     INTEGER I,J,K,K1,K2,KP1,NM1,NMJP1,ML,MU
C     DOUBLE PRECISION C1,C2,C3,C4,C5,C6,C7,C8,CONST1,CONST2,
C    1       CONST3,DELTA,FIVE,H,ONE,P5,SUM,SUM1,SUM2,
C    2       TEMP,TEMP1,TEMP2,TEN,TI,TK,TJ,TWO,ZERO
C     DOUBLE PRECISION AUX1(30),AUX2(30),PI(9),PIM1(9),Y(9)
C     DOUBLE PRECISION DFLOAT
      DOUBLE PRECISION C1,C2,C3,C4,C5,C6,C7,C8,C9
      DOUBLE PRECISION EIGHT,FIVE,H,ONE,P,PROD,PSQRT
      DOUBLE PRECISION R,RFACTR,SUM,SUM1,SUM2
      DOUBLE PRECISION TEMP,TEMP1,TEMP2,TEN,THREE,TI,TJ,TK
      DOUBLE PRECISION TOT,TPI,TWO,ZERO
      DATA ZERO,ONE,TWO,THREE,FIVE,EIGHT,TEN
     1    /0.D0,1.D0,2.D0,3.D0,5.D0,8.D0,10.D0/,
     2    C1,C2 /1.D4,10001.D-4/,
     3    C3,C4,C5,C6  /2.D2,202.D-1,198.D-1,180.D0/,
     4    C7 /25.D-2/ ,C8 /5.D-1/, C9 /2.9D1/
C    5     P5 /5.D-1/,
C     DFLOAT(IVAR) = IVAR
C
C     SELECTION OF PROBLEM.
C
      GO TO (10,20,30,40,50,60,120,170,200,220,270,300,330,350,
     1    400,410,420,420,420,420),NPROB
C
C     ROSENBROCK S FUNCTION.
C
   10 CONTINUE
      FVEC(1) = ONE - X(1)
      FVEC(2) = TEN*(X(2) - X(1)**2)
      GO TO 380
C
C     POWELL S SINGULAR FUNCTION.
C
   20 CONTINUE
      FVEC(1) = (X(1) + TEN*X(2))
      FVEC(2) = DSQRT(FIVE)*(X(3) - X(4))
      FVEC(3) = (X(2) - TWO*X(3))**2
      FVEC(4) = DSQRT(TEN)*(X(1) - X(4))**2
      GO TO 380
C
C     POWELL S BADLY SCALED FUNCTION.
C
   30 CONTINUE
      FVEC(1) = C1*X(1)*X(2) - ONE
      FVEC(2) = DEXP(-X(1)) + DEXP(-X(2)) - C2
      GO TO 380
C
C     WOOD
C
   40 CONTINUE
      TEMP1 = X(2) - X(1)**2
      TEMP2 = X(4) - X(3)*X(3)
      FVEC(1) = -C3*X(1)*TEMP1 - (ONE - X(1))
      FVEC(2) = C3*TEMP1 + C4*(X(2) - ONE) + C5*(X(4) - ONE)
      FVEC(3) = - C6*X(3)*TEMP2 - (ONE - X(3))
      FVEC(4) = C6*TEMP2 + C4*(X(4) - ONE) + C5*(X(2) - ONE)
      GO TO 380
C
C     HELICAL VALLEY
C
   50 CONTINUE
      TPI = EIGHT*DATAN(ONE)
      TEMP1 = DSIGN(C7,X(2))
      IF(X(1).GT.ZERO)TEMP1 = DATAN(X(2)/X(1))/TPI
      IF(X(1).LT.ZERO)TEMP1 = DATAN(X(2)/X(1))/TPI+C8
      TEMP2 = DSQRT(X(1)**2+X(2)**2)
      FVEC(1) = TEN*(X(3)-TEN*TEMP1)
      FVEC(2) = TEN*(TEMP2-ONE)
      FVEC(3) = X(3)
      GO TO 380
C
C     WATSON  FUNCTION.
C
   60 CONTINUE
      DO 70 K = 1, N
         FVEC(K) = ZERO
   70    CONTINUE
      DO 110 I = 1, 29
         TI = DBLE(FLOAT(I))/C9
         SUM1 = ZERO
         TEMP = ONE
         DO 80 J = 2,N
            SUM1 = SUM1+DBLE(FLOAT(J-1))*TEMP*X(J)
            TEMP = TI*TEMP
   80       CONTINUE
         SUM2 = ZERO
         TEMP = ONE
         DO 90 J = 1,N
            SUM2 = SUM2+TEMP*X(J)
            TEMP = TI*TEMP
   90       CONTINUE
         TEMP1 = SUM1-SUM2**2-ONE
         TEMP2 = TWO*TI*SUM2
         TEMP = ONE/TI
         DO 100 K = 1,N
            FVEC(K) = FVEC(K)+TEMP*(DBLE(FLOAT(K-1))-TEMP2)*TEMP1
            TEMP = TI*TEMP
  100       CONTINUE
  110    CONTINUE
      TEMP = X(2)-X(1)**2-ONE
      FVEC(1) = FVEC(1)+X(1)*(ONE-TWO*TEMP)
      FVEC(2) = FVEC(2)+TEMP
      GO TO 380
C
C     CHEBYQUAD FUNCTIONS.
C
  120 CONTINUE
      DO 130 K = 1,N
         FVEC(K) = ZERO
  130    CONTINUE
      DO 150 J = 1,N
         TEMP1 = ONE
         TEMP2 = TWO*X(J)-ONE
         TEMP = TWO*TEMP2
         DO 140 I = 1,N
            FVEC(I) = FVEC(I)+TEMP2
            TI = TEMP*TEMP2-TEMP1
            TEMP1 = TEMP2
            TEMP2 = TI
  140       CONTINUE
  150    CONTINUE
      TK = ONE/DBLE(FLOAT(N))
      IEV = -1
      DO 160 K = 1,N
         FVEC(K) = TK*FVEC(K)
         IF(IEV.GT.0)FVEC(K) = FVEC(K)+ONE/(DBLE(FLOAT(K))**2-ONE)
         IEV = -IEV
  160    CONTINUE
      GO TO 380
C
C     BROWN ALMOST LINEAR FUNCTION.
C
  170 CONTINUE
      SUM = -DBLE(FLOAT(N+1))
      PROD = ONE
      DO 180 J = 1,N
         SUM = SUM+X(J)
         PROD = PROD*X(J)
  180    CONTINUE
      DO 190 K = 1,N
         FVEC(K) = X(K)+SUM
  190    CONTINUE
      FVEC(N) = PROD-ONE
      GO TO 380
C
C     DISCRETE BOUNDARY VALUE FUNCTION.
C
  200 CONTINUE
      H = ONE/DBLE(FLOAT(N+1))
      DO 210 K = 1,N
         TEMP = (X(K)+DBLE(FLOAT(K))*H+ONE)**3
         TEMP1 = ZERO
         IF(K.NE.1)TEMP1 = X(K-1)
         TEMP2 = ZERO
         IF(K.NE.N)TEMP2 = X(K+1)
         FVEC(K) = TWO*X(K)-TEMP1-TEMP2+TEMP*H**2/TWO
  210    CONTINUE
      GO TO 380
C
C     DISCRETE INTEGRAL EQUATION FUNCTION.
C
  220 CONTINUE
      H = ONE/DBLE(FLOAT(N+1))
      DO 260 K = 1,N
         TK = DBLE(FLOAT(K))*H
         SUM1 = ZERO
         DO 230 J = 1,K
            TJ = DBLE(FLOAT(J))*H
            TEMP = (X(J)+TJ+ONE)**3
            SUM1 = SUM1+TJ*TEMP
  230       CONTINUE
         SUM2 = ZERO
         KP1 = K+1
         IF(N.LT.N+1)GO TO 250
         DO 240 J = KP1,N
            TJ = DBLE(FLOAT(J))*H
            TEMP = (X(J)+TJ+ONE)**3
            SUM2 = SUM2+(ONE-TJ)*TEMP
  240       CONTINUE
  250    CONTINUE
         FVEC(K) = X(K)+H*((ONE - TK)*SUM1+TK*SUM2)/TWO
  260    CONTINUE
      GO TO 380
C
C     TRIGONOMETRIC FUNCTION .
C
  270 CONTINUE
      SUM = ZERO
      DO 280 J = 1,N
         FVEC(J) = DCOS(X(J))
         SUM = SUM+FVEC(J)
  280    CONTINUE
      DO 290 K = 1,N
      FVEC(K) = DBLE(FLOAT(N+K))-DSIN(X(K))-SUM-DBLE(FLOAT(K))*FVEC(K)
  290    CONTINUE
      GO TO 380
C
C     VARIABLY DIMENSIONED FUNCTION.
C
  300 CONTINUE
      SUM = ZERO
      DO 310 J = 1, N
         SUM = SUM + DBLE(FLOAT(J))*(X(J) - ONE)
  310    CONTINUE
      TEMP = SUM*(ONE+TWO*SUM**2)
      DO 320 K = 1, N
         FVEC(K) = X(K) - ONE + TEMP*DBLE(FLOAT(K))
  320    CONTINUE
      GO TO 380
C
C     BROYDEN TRIDIAGONAL FUNCTION.
C
  330 CONTINUE
      DO 340 K = 1,N
         TEMP = (THREE-TWO*X(K))*X(K)
         TEMP1 = ZERO
         IF(K.NE.1)TEMP1 = X(K-1)
         TEMP2 = ZERO
         IF(K.NE.N)TEMP2 = X(K+1)
         FVEC(K) = TEMP-TEMP1-TWO*TEMP2+ONE
  340    CONTINUE
      GO TO 380
C
C     BROYDEN BANDED FUNCTION.
C
  350 CONTINUE
      ML = 5
      MU = 1
      TEMP = ZERO
      DO 370 K = 1, N
         K1 =MAX0(1, K - ML)
         K2 =MIN0( K + MU,N)
         TEMP = ZERO
         DO 360 J = K1, K2
            IF(J.NE.K)TEMP = TEMP+X(J)*(ONE+X(J))
  360       CONTINUE
         FVEC(K) = X(K)*(TWO+FIVE*X(K)**2)+ONE-TEMP
  370    CONTINUE
      GO TO 380
C
C     CHEMICAL EQUILIBRIUM PROBLEM 1
C
  400 CONTINUE
      FVEC(1) = X(2)-TEN
      FVEC(2) = X(1)-X(2)-5.D4
      GO TO 380
C
C     CHEMICAL EQUILIBRIUM PROBLEM 2
C
  410 CONTINUE
      FVEC(1) = X(1)+X(2)+X(4)-1.D-3
      FVEC(2) = X(5)+X(6)-55.D0
      FVEC(3) = X(1)+X(2)+X(3)+2.D0*X(5)+X(6)-110.001D0
      FVEC(4) = X(1)-.1D0*X(2)
      FVEC(5) = X(1)-1.D4*X(3)*X(4)
      FVEC(6) = X(5)-55.D14*X(3)*X(6)
      GO TO 380
C
C     CHEMICAL EQUILIBRIUM PROBLEM 3
C
  420 CONTINUE
      R = RFACTR
      P = ONE
      PSQRT = DSQRT(P)
      TOT = ZERO
      DO 425 I = 1,10
         TOT = TOT+X(I)
  425    CONTINUE
      FVEC(1) = X(1)+X(4)-THREE
      FVEC(2) = TWO*X(1)+X(4)+X(2)+X(7)+X(8)+X(9)+TWO*X(10)-R
      FVEC(3) = TWO*(X(2)+X(5))+X(6)+X(7)-EIGHT
      FVEC(4) = TWO*X(3)+X(9)-4.D0*R
      FVEC(5) = X(1)*X(5)-1.93D-1*X(2)*X(4)
      FVEC(6) = X(6)*DSQRT(X(1)*P)-2.597D-3*DSQRT(X(2)*X(4)*TOT)
      FVEC(7) = X(7)*DSQRT(X(4)*P)-3.448D-3*DSQRT(X(1)*X(2)*TOT)
      FVEC(8) = X(8)*X(4)*P-1.799D-5*X(1)*TOT
      FVEC(9) = X(9)*X(4)*PSQRT-2.155D-4*X(1)*DSQRT(X(3)*TOT)
      FVEC(10) = X(10)*X(4)*X(4)*P-3.846D-5*X(4)*X(4)*TOT
C
  380 CONTINUE
      RETURN
C
C     LAST CARD OF SUBROUTINE VECFCN.
C
      END
      SUBROUTINE TRANGE (N,X,NPROB,IER)
C
C THE SUBROUTINE  TRANGE  IS DESIGNED TO AVOID FAILURES (E.G. OVERFLOWS)
C DURING THE COMPUTATION OF FUNCTION VALUES IN SOME OF THE MOST
C DANGEROUS FUNCTIONS IN THE SUBROUTINE  VECFCN .
C TRANGE  IS CALLED BY THE SUBROUTINES  VECF  AND  JCBNA  (TEST
C VERSIONS) AND RETURNS A NEGATIVE VALUE OF THE
C ERROR INDICATOR  IER  IF THE POINT  X  DOES NOT BELONG TO AN
C ADMISSIBLE REGION.
C THE ADMISSIBLE REGION IS DEFINED BY
C     FOR ALL FUNCTIONS   ABS(X(I)).LE.XMAX ,   I = 1,...,N
C     FOR FUNCTION 3           X(I).GE.P3MIN ,  I = 1,...,N
C     FOR FUNCTION 7      ABS(X(I)).LE.P7MAX ,  I = 1,...,N
C     FOR FUNCTION 8    ABS(X(1)*X(2)*...*X(I)).LE.P8MAX ,  I = 1,...,N
C     FOR FUNCTION 9      ABS(X(I)).LE.P9MAX ,  I = 1,...,N
C
      DOUBLE PRECISION X,TOT,XMAX,P3MIN,P7MAX,P8MAX,P9MAX
      DIMENSION X(N)
C
C USE THE FIRST VALUE OF XMAX (1D8) OR THE SECOND (1D4) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1D.120 (E.G. UNIVAC 1100), OR BELOW
C 1D.120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA XMAX/1.D8/
C     DATA XMAX/1.D4/
C
C USE THE FIRST VALUE OF P3MIN (-184) IF THE MACHINE DYNAMIC RANGE
C IS ABOVE 1.D280 (E.G. UNIVAC 1100) THE THIRD VALUE (-24)  IF BELOW
C 1.D45 (E.G. DEC VAX, PDP11), AND THE SECOND VALUE (-46) IF
C IN BETWEEN (E.G. AMDAHL 470, IBM 360/370).
      DATA P3MIN/-184.D0/
C     DATA P3MIN/-46.D0/
C     DATA P3MIN/-23.D0/
C
C USE THE FIRST VALUE OF P7MAX (1D4) OR THE SECOND (1D2) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1D.120 (E.G. UNIVAC 1100), OR BELOW
C 1D.120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA P7MAX/1.D4/
C     DATA P7MAX/1.D2/
C
C USE THE FIRST VALUE OF P8MAX (1.D80) IF THE MACHINE DYNAMIC RANGE
C IS ABOVE 1.D280 (E.G. UNIVAC 1100) THE THIRD VALUE (1.D10)  IF BELOW
C 1.D45 (E.G. DEC VAX, PDP11), AND THE SECOND VALUE (1.D20) IF
C IN BETWEEN (E.G. AMDAHL 470, IBM 360/370).
      DATA P8MAX/1.D80/
C     DATA P8MAX/1.D20/
C     DATA P8MAX/1.D10/
C
C USE THE FIRST VALUE OF P9MAX (1D4) OR THE SECOND (1D2) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1D.120 (E.G. UNIVAC 1100), OR BELOW
C 1D.120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA P9MAX/1.D4/
C     DATA P9MAX/1.D2/
C
      IER = 0
      DO 10 I = 1,N
         IF (DABS(X(I)).GT.XMAX) IER = -1
  10     CONTINUE
      GO TO (20,30,40,50,60,70,80,100,130,150,160,170,180,190,200,210,
     1    220), NPROB
C
  20  CONTINUE
      RETURN
  30  CONTINUE
      RETURN
  40  CONTINUE
      IF (DMIN1(X(1),X(2)).LT.P3MIN) IER = -1
      RETURN
  50  CONTINUE
      RETURN
  60  CONTINUE
      RETURN
  70  CONTINUE
      RETURN
  80  CONTINUE
      DO 90 I = 1,N
         IF (DABS(X(I)).GT.P7MAX) IER = -1
  90     CONTINUE
      RETURN
 100  CONTINUE
      TOT = 1.D0
      DO 110 I = 1,N
         TOT = TOT*X(I)
         IF (DABS(TOT).GT.P8MAX) GO TO 120
 110     CONTINUE
      RETURN
 120  CONTINUE
      IER = -1
      RETURN
 130  CONTINUE
      DO 140 I = 1,N
         IF (DABS(X(I)).GT.P9MAX) IER = -1
 140     CONTINUE
      RETURN
 150  CONTINUE
      RETURN
 160  CONTINUE
      RETURN
 170  CONTINUE
      RETURN
 180  CONTINUE
      RETURN
 190  CONTINUE
      RETURN
 200  CONTINUE
      RETURN
 210  CONTINUE
      RETURN
 220  CONTINUE
      TOT = 0.D0
      DO 230 I = 1,10
         TOT = TOT+X(I)
 230     CONTINUE
      IF (X(1).LT.0.D0.OR.X(2)*TOT.LT.0.D0) IER = -1
      IF (X(4).LT.0.D0.OR.X(3)*TOT.LT.0.D0) IER = -1
      RETURN
      END
      SUBROUTINE DAFNE (N,X,FMIN,PRF,IDER,NU,NO,IDA,
     1    NFEV,NJEV,NSAC,NST,IER,
     2    V,FP,VFU,BUF,XTP,WF,FPS,S,C,CS)
C
C THE SUBROUTINE  DAFNE  ATTEMPTS TO SOLVE A SYSTEM
C OF  N  NONLINEAR EQUATIONS
C     FI(X) = 0,  I = 1,...,N,  X = (X1,...,XN)        (1)
C BY MEANS OF THE DIFFERENTIAL-EQUATION APPROACH DESCRIBED
C WITHIN THE SUBROUTINE  EVOLV .
C THE NUMERICAL INTEGRATION AND RESCALING CAPABILITY OF  EVOLV
C IS EXPLOITED BY  DAFNE  TO OBTAIN A TRAJECTORY WHICH IS AT
C TIMES STOPPED TO RESCALE THE FUNCTION AND THEN RESTARTED FROM
C THE POINT REACHED.
C THIS IS IMPLEMENTED BY PERFORMING
C A SEQUENCE OF CALLS TO THE SUBROUTINE  EVOLV ,
C SETTING SUITABLE VALUES FOR THE CALL PARAMETERS,
C EACH TIME USING AS INITIAL POINT THE FINAL POINT
C RETURNED BY THE PRECEDING CALL, AND COMPUTING
C A NEW VECTOR OF WEIGHTS
C     WF = (WF1,WF2,...,WFN)
C USED FOR FUNCTION RESCALING.
C THE RESCALING WEIGHTS ARE COMPUTED BY CALLING
C THE SUBROUTINE  NEWF .
C
C CALL STATEMENT :
C     CALL DAFNE (N,X,FMIN,PRF,IDER,NU,NO,IDA,
C    1    NFEV,NJEV,NSAC,NST,IER,
C    2    B1,B2,B3,B4,B5,B6,B7,B8,C1,C2)
C
C DESCRIPTION OF PARAMETERS OF THE CALL STATEMENT
C     (FORTRAN IMPLICIT TYPE DEFINITION IS USED.
C     ALL NON-INTEGER PARAMETERS ARE DOUBLE-PRECISION)
C
C N  IS THE (INPUT) NUMBER OF NONLINEAR EQUATIONS.
C
C X  IS A VECTOR OF LENGTH  N  CONTAINING THE X-VARIABLES . (INITIAL
C     VALUES ON INPUT, FINAL ON OUTPUT. INITIAL AND FINAL  X  ARE
C     THE INITIAL AND FINAL ESTIMATES OF THE SOLUTION OF SYSTEM (1)).
C
C FMIN  IS  AN INPUT VECTOR OF DIMENSION  N  CONTAINING THE THRESHOLD
C     VALUES FOR THE FUNCTION VALUES  FI(X) , I = 1,...,N .
C     CONVERGENCE IS CLAIMED IF
C            ABS(FI(X)).LT.FMIN(I) , I = 1,...,N
C
C THE FOLLOWING 5 PARAMETERS (PRF,IDER,NU,NO,IDA)
C ARE USED ONLY IN THE COMPUTATION OF THE JACOBIAN MATRIX
C ( SUBROUTINE  JCBN ).
C
C PRF  IS AN (INPUT) ROUGH ESTIMATE OF THE AVERAGE RELATIVE
C     ERROR AFFECTING THE FUNCTION VALUES SUPPLIED BY SUBROUTINE  VECF .
C     IF THIS ERROR IS ESTIMATED TO BE NOT SEVERAL ORDERS OF MAGNITUDE
C     GREATER THAN MACHINE PRECISION, THE USER MAY SIMPLY PUT  PRF = 0 .
C
C IDER  IS AN INPUT CONTROL VARIABLE TO REQUEST THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDER.EQ.0   ANALYTICAL COMPUTATION BY SUBROUTINE  JCBNA
C     IDER.EQ.1   COMPUTATION BY FORWARD FINITE DIFFERENCES
C     IDER.EQ.2   COMPUTATION BY CENTRAL FINITE DIFFERENCES
C
C NU  AND  NO  ARE NUMBERS DEFINING THE (POSSIBLY) BANDED
C STRUCTURE OF THE JACOBIAN MATRIX, USED FOR FINITE-DIFFERENCES
C JACOBIAN COMPUTATION
C     NU IS THE N. OF NON-ZERO JACOBIAN DIAGONALS UNDER MAIN DIAGONAL
C     NO IS THE N. OF NON-ZERO JACOBIAN DIAGONALS ABOVE MAIN DIAGONAL
C     (FOR A NON-BANDED MATRIX  NU = NO = N-1 )
C
C IDA  IS AN (OUTPUT) INDICATOR OF THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDA.EQ.1   THE DERIVATIVES WERE COMPUTED ANALYTICALLY
C                (AS REQUESTED)
C     IDA.EQ.0   THE DERIVATIVES WERE COMPUTED BY FINITE DIFFERENCES.
C                (CENTRAL FINITE DIFFERENCES IF REQUESTED BY  IDER = 2,
C                FORWARD FINITE DIFFERENCES IF REQUESTED BY  IDER = 1,
C                OR IF REQUESTED ANALYTICAL DERIVATIVES WERE NOT CODED
C                IN SUBROUTINE  JCBNA ).
C
C THE FOLLOWING 4 PARAMETERS (NFEV,NJEV,NSAC,NST) ARE OUTPUT COUNTERS
C     NFEV  OF FUNCTION EVALUATIONS (EXCLUDING THOSE USED
C           FOR JACOBIAN EVALUATION).
C     NJEV  OF JACOBIAN EVALUATIONS.
C     NSAC  OF ACCEPTED TIME INTEGRATION STEPS.
C     NST   OF ATTEMPTED (ACCEPTED + REJECTED) TIME INTEGRATION STEPS.
C
C IER  IS AN OUTPUT ERROR INDICATOR AS FOLLOWS
C     IER = 0    SUCCESS
C     IER = 1    FAILURE (MAX. NUMBER OF ATTEMPTS)
C     IER = 2    FATAL FAILURE IN COMPUTING THE FUNCTION VALUES
C                        AT THE INITIAL POINT OR THE DERIVATIVES.
C
C THE LAST 10 PARAMETERS OF THE CALL STATEMENT
C ARE WORKING STORAGE ARRAYS
C     B1,...,B8  OF DIMENSION  N
C     C1,C2  OF DIMENSION  (N,N)
C NOTE THAT IN THE SUBROUTINE DEFINITION STATEMENT
C THESE ARRAYS ARE NAMED ACCORDING TO THEIR INTERNAL USE.
C THE OUTPUT CONTENT OF THESE ARRAYS SHOULD NOT CONCERN
C THE AVERAGE USER.
C THE INTERESTED USER WILL FIND ON OUTPUT IN  B1, B3, B5, B6
C THE FINAL VALUES OF THE VECTORS  V , VFU , XTP , WF
C ( V  AND  VFU  ARE DESCRIBED IN SUBROUTINE  EVOLV , AND  XTP
C IN SUBROUTINE  JCBN ).
C
C CALL TO  DAFNE
C
C THE PROGRAM CALLING  DAFNE  MUST SET THE INPUT
C CALL PARAMETERS  N , X , FMIN , PRF , IDER , NU , NO .
C
C RETURN FROM DAFNE
C
C THE CALLING PROGRAM RECEIVES FROM  DAFNE  THE OUTPUT
C PARAMETERS   X , IDA , NFEV , NJEV , NSAC , NST , IER .
C
C SUBPROGRAMS CALLED  :  EVOLV , NEWF .
C
C THE USER MUST PROVIDE THE FOLLOWING TWO SUBPROGRAMS
C (FORTRAN IMPLICIT TYPE DEFINITION IS USED.
C ALL NON-INTEGER PARAMETERS ARE DOUBLE-PRECISION)
C
C SUBROUTINE  VECF (N,X,VFU,NR,IER)
C THE SUBROUTINE  VECF  MUST RETURN IN THE N-VECTOR  VFU
C THE (NON-RESCALED) VALUES OF THE FUNCTIONS OF SYSTEM (1)
C COMPUTED AT THE POINT  X .
C NR  IS THE (INPUT) NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     MUST BE USED TO INDICATE THE FACT THAT  X  IS OUTSIDE A SAFETY
C     REGION DEFINED BY THE USER).
C
C SUBROUTINE  JCBNA (N,X,C,NR,IDA,IER)
C THE SUBROUTINE  JCBNA  IS DESIGNED TO GIVE THE USER
C THE OPPORTUNITY OF PROVIDING THE ANALYTICAL JACOBIAN.
C IN THIS CASE THE USER MUST CODE THE SUBROUTINE ACCORDING
C TO THE FOLLOWING SPECIFICATIONS.
C THE SUBROUTINE  JCBNA  MUST RETURN IN THE (N,N) ARRAY  C
C THE JACOBIAN OF THE FUNCTIONS OF SYSTEM (1)
C COMPUTED AT THE POINT  X .
C NR  IS THE (INPUT) NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C IDA  MUST BE SET EQUAL TO 1.
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     MUST BE USED TO INDICATE THE FACT THAT  X  IS OUTSIDE A SAFETY
C     REGION DEFINED BY THE USER).
C IF THE USER DOES NOT WANT TO GIVE THE ANALYTICAL JACOBIAN,
C THE SUBROUTINE  JCBNA  MUST ONLY SET  IDA  EQUAL TO 0.
C
      DOUBLE PRECISION AMU,AMU0,BE,BE0,BUF,C,CS,FMIN,FP,FPS,H,H0
      DOUBLE PRECISION S,PRF,V,VFU,V0,WF,X,XTP
      DIMENSION X(N),FMIN(N),V(N),FP(N),VFU(N),BUF(N),XTP(N),WF(N)
      DIMENSION FPS(N),S(N),C(N,N),CS(N,N)
      DIMENSION NSC(8)
C
C NTRIAL  IS THE DIMENSION OF THE ARRAY  NSC  CONTAINING
C     THE MAX. NUMBER OF STEPS FOR THE SUCCESSIVE CALLS TO  EVOLV
      DATA NTRIAL/8/
      DATA NSC(1),NSC(2),NSC(3),NSC(4),NSC(5),NSC(6),NSC(7),NSC(8)
     1    /  1,     2,     4,     8,    16,     32,   64,    128 /
C
C INITIAL VALUE OF TIME INTEGRATION STEPLENGTH
      H0 = 1.D0
C INITIAL VALUES OF MASS AND FRICTION COEFFICIENTS
      AMU0 = 1.D0
      BE0 = 1.D0
C INITIAL VALUE FOR THE VELOCITIES
      V0 = 0.D0
C
C INITIALIZE FUNCTION WEIGHTS
      DO 10 I = 1,N
         WF(I) = 1.D0
  10     CONTINUE
C
C INITIALIZE TIME INTEGRATION STEPLENGTH
      H = H0
C
C INITIALIZE COUNTERS
      NFEV = 0
      NJEV = 0
      NSAC = 0
      NST = 0
C
C START SEQUENCE OF CALL TO  EVOLV
      DO 30 ITRIAL = 1,NTRIAL
C
C        SET MAX. NUMBER OF STEPS FOR NEXT CALL TO  EVOLV
         NIT = NSC(ITRIAL)
C
C        INITIALIZE VELOCITIES, MASS, AND FRICTION COEFFICIENT
         DO 20 I = 1,N
            V(I) = V0
  20        CONTINUE
         AMU = AMU0
         BE = BE0
C
C        CALL THE BASIC EQUATION-SOLVING SUBROUTINE
         CALL EVOLV (N,X,V,AMU,BE,H,FMIN,NIT,NR,II,IER,
     1       PRF,IDER,NU,NO,IDA,NFEV,NJEV,
     2       FP,C,VFU,BUF,XTP,WF,FPS,CS,S)
C
C        UPDATE COUNTERS
         NSAC = NSAC+II
         NST = NST+NR
C
         IF (IER.EQ.0) GO TO 60
         IF (IER.EQ.4.OR.IER.EQ.5) GO TO 40
C
C        COMPUTE FUNCTION WEIGHTS
         NRP = NR+1
         NFEV = NFEV+1
         CALL NEWF (N,X,WF,PRF,IDER,NU,NO,IDA,
     1       XTP,VFU,C,FP,BUF,NRP,NJEV,IER)
         IF (IER.LT.0) GO TO 50
C
  30     CONTINUE
C
C ERROR EXIT (MAX. NUMBER OF ATTEMPTS REACHED)
      IER = 1
      RETURN
C
C ERROR EXIT (FATAL FAILURE IN COMPUTING THE FUNCTION VALUES
C AT THE INITIAL POINT OR THE DERIVATIVES)
  40  CONTINUE
  50  CONTINUE
      IER = 2
      RETURN
C
C NORMAL EXIT
  60  CONTINUE
      IER = 0
C
      RETURN
      END
      SUBROUTINE CHAH (NSR,VF,GA,IER,II,NR)
C
C THE SUBROUTINE  CHAH  COMPUTES THE VARIATION OF THE TIME INTEGRATION
C STEPLENGTH, AND MAY DECIDE TO REJECT THE LAST STEP PERFORMED.
C
C DESCRIPTION OF PARAMETERS
C
C NSR  IS THE NUMBER OF STEPS FOR WHICH TOTAL ENERGY IS REMEMBERED.
C VF  IS AN ARRAY OF LENGTH  NSR  CONTAINING THE LAST VALUES
C     OF THE TOTAL ENERGY.
C GA  IS THE (OUTPUT) MULTIPLICATIVE FACTOR RETURNED TO UPDATE
C     THE TIME INTEGRATION STEPLENGTH.
C IER  IS AN (INPUT OR OUTPUT) ERROR INDICATOR.
C     A NEGATIVE VALUE INDICATES THAT THE LAST STEP IS TO BE
C     REJECTED.
C II  IS THE NUMBER OF TIME INTEGRATION STEPS.
C NR  IS THE TOTAL NUMBER OF ATTEMPTED TIME INTEGRATION STEPS,
C     INCLUDING THE STEPS ATTEMPTED BUT REJECTED.
C
      DOUBLE PRECISION GA,GAB,GAMAX,GAMIN,GAR,ONE,TOL,UGAMAX,UGAR,VF
      DOUBLE PRECISION GAEPS
      DIMENSION VF(NSR)
C
      DATA GAMIN/2.D-2/
      DATA ILB/-20/
      DATA GAMAX /1.D+2/
      DATA GAB /25.D0/
      DATA ONE /1.D0/
      DATA TOL/1.D-10/
      DATA GAR/1.1D0/
      DATA ISR/0/
      DATA GAEPS/1.D-10/
C
      IF (NR.EQ.1) ILB = -20
      IF (NR.EQ.1) ISR = 0
      IF (IER.LT.0) GO TO 80
      UGAMAX = GAMAX
      IF (II.GT.NSR) UGAMAX = GAB
      IF (ILB+5.GT.II) UGAMAX = ONE
      UGAR = DMIN1(GAR,UGAMAX)
      IP = 0
      IM = 0
      DO 10 I = 2,NSR
         IF (IP.EQ.0.AND.VF(I).GT.VF(1)*(1.D0+TOL)) IP = I-1
         IF (IM.EQ.0.AND.VF(I).LT.VF(1)*(1.D0-TOL)) IM = I-1
         IF (IP.NE.0.AND.IM.NE.0) GO TO 20
  10     CONTINUE
  20  CONTINUE
      IF (IP.EQ.1.AND.IM.EQ.0) GO TO 70
      IF (IM.EQ.1) GO TO 30
      IF (IP.EQ.1) GO TO 40
      IF (IP.GT.1.AND.IM.GT.1) GO TO 50
      IF (IP.GT.1) GO TO 60
      GA = UGAR
      GO TO 90
  30  CONTINUE
      IF (IP.EQ.0) GO TO 80
      GA = ONE
      IF ((IP-1)**2.GT.II) GA = GAMIN
      GO TO 90
  40  CONTINUE
      GA = ONE
      IF (II.GT.NSR) GA = UGAR
      IF ((IM-1)**2.GT.II*4) GA = UGAMAX
      GO TO 90
  50  CONTINUE
      GA = UGAMAX
      GO TO 90
  60  CONTINUE
      GA = UGAMAX
      GO TO 90
  70  CONTINUE
      GA = UGAMAX
      GO TO 90
  80  CONTINUE
      ISR = ISR+1
      IER = -3
      GA = GAMIN
      IF (ISR.GE.3) GA = GAEPS
      ILB = II
      RETURN
  90  CONTINUE
      IER = 1
      ISR = 0
      RETURN
      END
      SUBROUTINE CREXTP (N,X,SPRR,XTP)
C
C THE SUBROUTINE  CREXTP  INITIALIZES THE ARRAY  XTP  (SEE BELOW)
C
C DESCRIPTION OF PARAMETERS
C
C N  AND  X  ARE DEFINED IN SUBROUTINE  EVOLV .
C SPRR  IS A VARIABLE WHOSE VALUE IS COMPUTED BY SUBROUTINE
C     JCBN , DEPENDING ON THE ESTIMATED NUMERICAL PRECISION.
C XTP  IS AN ARRAY OF LENGTH  N  USED BY
C     SUBROUTINE  JCBN  TO COMPUTE THE STEPLENGTHS USED
C     FOR COMPUTING THE FINITE-DIFFERENCES JACOBIAN.
C
      DOUBLE PRECISION XTP,SPRR,X
      DIMENSION X(N),XTP(N)
C
      DO 10 I = 1,N
         XTP(I) = DMAX1(SPRR,DABS(X(I)))
  10     CONTINUE
      RETURN
      END
      SUBROUTINE DERIV (N,X,VFU,FP,C,WF,NR,NJEV,IER,
     1    PRF,IDER,NU,NO,IDA,BUF,XTP)
C
C THE SUBROUTINE  DERIV RETURNS IN THE N-VECTOR  FP
C THE GRADIENT OF THE WEIGHTED SUM
C     FW = W1*F1(X)**2 + ... + WN*FN(X)**2
C OF THE SQUARES OF THE LEFT-HAND SIDE FUNCTIONS OF SYSTEM (1)
C (SEE  EVOLV ), AND IN THE (N BY N) MATRIX  C  THE (APPROXIMATE)
C HESSIAN MATRIX
C
C          T
C     2 * J  * W * J
C
C OF  FW , COMPUTED AT THE POINT  X , WHERE  W  IS
C THE DIAGONAL MATRIX OF THE RESCALING
C WEIGHTS, AND  J  IS THE JACOBIAN MATRIX OF THE (NON-RESCALED)
C FUNCTIONS IN (1), OBTAINED BY CALLING THE SUBROUTINE  JCBN .
C VFU  IS AN INPUT ARRAY OF LENGTH  N  CONTAINING THE VALUES OF
C     THE FUNCTIONS IN SYSTEM (1).
C WF  IS THE N-VECTOR CONTAINING THE RESCALING WEIGHTS.
C NR  IS THE (INPUT) NUMBER OF THE TIME INTEGRATION STEP IN
C     PROGRESS (COUNTING ALSO THE REJECTED STEPS).
C NJEV  IS THE COUNTER OF THE JACOBIAN EVALUATIONS
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     INDICATES A FAILURE IN COMPUTING THE JACOBIAN).
C PRF , IDER , NU , NO , IDA  ARE DESCRIBED WITHIN THE
C     SUBROUTINE  DAFNE .
C BUF  IS A WORKING STORAGE ARRAY OF LENGTH  N
C XTP  IS DESCRIBED IN THE SUBROUTINE  JCBN .
C
C SUBPROGRAMS CALLED  :  JCBN
C
      DOUBLE PRECISION BUF,C,FP,PRF,VFU,WF,X,XTP
      DIMENSION X(N),C(N,N),FP(N),VFU(N),BUF(N),XTP(N),WF(N)
C
C COMPUTE THE JACOBIAN MATRIX
      CALL JCBN (N,X,VFU,C,FP,BUF,XTP,PRF,IDER,NU,NO,IDA,
     1    NR,NJEV,IER)
      IF (IER.LT.0) RETURN
C
C COMPUTE THE GRADIENT VECTOR  FP
      DO 20 I = 1,N
         FP(I) = 0.D0
         DO 10 J = 1,N
            FP(I) = FP(I)+VFU(J)*C(J,I)*WF(J)
  10        CONTINUE
         FP(I) = FP(I)*2.D0
  20     CONTINUE
C
      DO 70 K = 1,N
C
C     COMPUTE COLUMN  K  OF THE APPROX. HESSIAN AND
C     OVERWRITE IT ON COLUMN  K  OF THE JACOBIAN.
C
C        COMPUTE THE BELOW-DIAGONAL PART OF COLUMN  K ...
         DO 40 J = K,N
            BUF(J) = 0.D0
            DO 30 I = 1,N
               BUF(J) = BUF(J)+C(I,J)*C(I,K)*WF(I)
  30           CONTINUE
  40        CONTINUE
C
C        ... AND OVERWRITE IT ON COLUMN  K  OF THE JACOBIAN
         DO 50 J = K,N
            C(J,K) = BUF(J)*2.D0
  50        CONTINUE
C
C        FILL BY SYMMETRY THE ABOVE-DIAGONAL PART OF COLUMN  K
         DO 60 J = 1,K
            C(J,K) = C(K,J)
  60        CONTINUE
C
  70     CONTINUE
      RETURN
      END
      SUBROUTINE DLIN (N,FP,C,FPS,CS,S,IER)
C
C THE SUBROUTINE DLIN COMPUTES THE MINIMUM-NORM
C LEAST-SQUARES SOLUTION OF THE LINEAR SYSTEM
C WITH POSSIBLY SINGULAR  N BY N  COEFFICIENT MATRIX  C
C AND RIGHT-HAND SIDE N-VECTOR  FP ,  AND RETURNS
C IN  FP  THE SOLUTION VECTOR (THE CONTENT OF  C
C MAY BE ALTERED).
C FPS , CS , S  ARE WORKING STORAGE ARRAYS OF LENGTH N , (N,N) , N .
C IER IS AN OUTPUT INDICATOR AS FOLLOWS
C     IER = -1       SOLUTION NOT FOUND
C     IER = 0        SOLUTION FOUND WITHOUT REGULARIZATION
C     IER = 1        SOLUTION FOUND BY REGULARIZATION
C
C SUBPROGRAMS CALLED  :  QMS , MOVM , DLINX .
C
      DOUBLE PRECISION C,FP
      DOUBLE PRECISION CS,D,DEPS1,DEPS2
      DOUBLE PRECISION DMEM,DMIN,DMAX,DSTEP,FPS,S,SP,SPV,SR
      DIMENSION FP(N),C(N,N)
      DIMENSION FPS(N),CS(N,N),S(N)
C
C THE REGULARIZATION PARAMETER  D  STARTS FROM THE PRECEDING
C VALUE  DMEM  AND IS REPEATEDLY UPDATED BY A FACTOR  DSTEP
C NOT BEYOND THE BOUNDARY VALUES  DMIN , DMAX .
C
C USE THE FIRST VALUES (-60,60) OR THE SECOND ONES (-30,30) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1.D120 (E.G. UNIVAC 1100), OR BELOW
C 1.D120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA DMIN,DMAX/1.D-60,1.D60/
C     DATA DMIN,DMAX/1.D-30,1.D30/
C
      DATA DSTEP/10.D0/
      DATA DMEM/1.D0/
C
C DEPS1  AND  DEPS2  ARE TOLERANCE FACTORS USED BY SUBROUTINE  DLINX .
      DATA DEPS1,DEPS2/1.D-16,1.D-16/
C
      DATA NMAX/1000/
C
C RESCALE BY NORMALIZING MAIN DIAGONAL ELEMENTS TO ONE .
C
C     COMPUTE SCALE FACTORS
      DO 10 I = 1,N
         S(I) = 1.D0/DSQRT(C(I,I))
  10     CONTINUE
C
      DO 30 I = 1,N
C
C        RESCALE COEFFICIENT MATRIX
         DO 20 J = 1,N
            C(I,J) = C(I,J)*S(I)*S(J)
  20        CONTINUE
C
C        RESCALE RIGHT-HAND-SIDE VECTOR
         FP(I) = FP(I)*S(I)
C
  30     CONTINUE
C
C     SAVE  FP  AND  C  IN  FPS  AND  CS
      CALL MOVM (N,FP,C,FPS,CS)
C
C     TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C     SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C     SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C     WITH TOLERANCE PARAMETER  DEPS1 .
      CALL DLINX (N,FPS,CS,DEPS1,IER)
C
      IF (IER.GT.0) IER = 0
      IF (IER.EQ.0) GO TO 140
C
C     THE LINEAR SYSTEM IS SINGULAR
C     START SEARCH FOR REGULARIZATION FACTOR  D
C     STARTING FROM PRECEDING VALUE  DMEM .
      D = DMEM
C
  40  CONTINUE
C
C DECREASE REGULARIZATION PARAMETER
      D = D/DSTEP
C
      DO 50 I = 1,NMAX
C
C        REGULARIZE THE LINEAR SYSTEM  C ,  FP , WITH REGULARIZATION
C        PARAMETER  D , OBTAINING THE LINEAR SYSTEM  FPS ,  CS .
         CALL QMS (N,FP,C,FPS,CS,D)
C
C        TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C        SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C        SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C        WITH TOLERANCE PARAMETER  DEPS2 .
         CALL DLINX (N,FPS,CS,DEPS2,IER)
C
         IF (IER.GE.0.AND.I.EQ.1.AND.D.GT.DMIN) GO TO 40
         IF (IER.GE.0) GO TO 80
         IF (D.GT.DMAX) GO TO 60
C
C        INCREASE THE REGULARIZATION PARAMETER
         D = D*DSTEP
C
  50     CONTINUE
C
  60  CONTINUE
  70  CONTINUE
C
C ERROR EXIT
      IER = -1
      RETURN
C
  80  CONTINUE
C INCREASE THE REGULARIZATION PARAMETER
      DMEM = D*DSTEP
C
C     THE LINEAR SYSTEM IS ESTIMATED NONSINGULAR BY SUBROUTINE  DLINX
C     LOOK FOR LEAST-SQUARES SOLUTION
      SP = 0.D0
C
      DO 110 I = 1,NMAX
         SPV = SP
         SP = 0.D0
         DO 100 J = 1,N
            SR = 0.D0
            DO 90 K = 1,N
               SR = SR+C(J,K)*FPS(K)
  90           CONTINUE
            SR = SR-FP(J)
            SP = SP+SR*SR
 100        CONTINUE
         IF (SP.GT.SPV.AND.I.GT.1) GO TO 120
C
C        INCREASE THE REGULARIZATION PARAMETER
         D = D*DSTEP
         IF (D.GT.DMAX) GO TO 130
C
C        REGULARIZE THE LINEAR SYSTEM  C ,  FP , WITH REGULARIZATION
C        PARAMETER  D , OBTAINING THE LINEAR SYSTEM  FPS ,  CS .
         CALL QMS (N,FP,C,FPS,CS,D)
C
C        TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C        SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C        SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C        WITH TOLERANCE PARAMETER  DEPS2 .
         CALL DLINX (N,FPS,CS,DEPS2,IER)
C
         IF (IER.LT.0) GO TO 70
C
 110     CONTINUE
C
 120  CONTINUE
      D = D/DSTEP
C
C     REGULARIZE THE LINEAR SYSTEM  C ,  FP , WITH REGULARIZATION
C     PARAMETER  D , OBTAINING THE LINEAR SYSTEM  FPS ,  CS .
      CALL QMS (N,FP,C,FPS,CS,D)
C
C     TRY TO SOLVE THE SUPPOSED NON-SINGULAR LINEAR
C     SYSTEM OF COEFFICIENT MATRIX  CS  AND RIGHT-HAND-
C     SIDE VECTOR  FPS , BY CALLING THE SUBROUTINE  DLINX
C     WITH TOLERANCE PARAMETER  DEPS2 .
      CALL DLINX (N,FPS,CS,DEPS2,IER)
C
 130  CONTINUE
      IER = 1
 140  CONTINUE
C
C REVERSE SCALING FOR VECTOR  FP
      DO 150 I = 1,N
         FP(I) = FPS(I)*S(I)
 150     CONTINUE
C
      RETURN
C
      END
      SUBROUTINE DLINX (N,B,A,DEPS,IER)
C
C THE SUBROUTINE  DLINX  ATTEMPTS TO COMPUTE THE SOLUTION
C OF THE LINEAR SYSTEM WITH  N BY N  COEFFICIENT MATRIX  A
C AND RIGHT-HAND SIDE N-VECTOR  B ,  AND RETURNS
C IN  B  THE SOLUTION VECTOR (THE CONTENT OF  A
C MAY BE ALTERED).
C DEPS  IS AN INPUT TOLERANCE FACTOR USED TO DETECT FAILURE.
C IER  IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE INDICATES FAILURE)
C
      DOUBLE PRECISION A,AM,DEPS,P,PI,B,Z,ERR
C
      DIMENSION A(N,N),B(N)
      IF(N.LE.0)GO TO 230
C
      IER = 0
      IRM = 0
      ICM = 0
      P = 0.D0
      DO 30 IR = 1,N
         DO 10 IC = 1,N
            AM = DABS(A(IR,IC))
            IF(AM.LE.P)GO TO 10
            P = AM
            IRM = IR
            ICM = IC
  10        CONTINUE
  30     CONTINUE
      IF(IRM.EQ.0.OR.ICM.EQ.0)GO TO 240
      ERR = DEPS*P
C
      DO 170 K = 1,N
C
C
         IF(P.LE.0.D0)GO TO 230
         IF(IER.NE.0)GO TO 70
         IF(P.GT.ERR)GO TO 70
         IER = K-1
  70     PI = 1.D0/A(IRM,ICM)
         I = IRM - K
         J = ICM - K
C
         Z = PI*B(IRM)
         B(IRM) = B(K)
         B(K) = Z
C
         IF(K.GE.N)GO TO 180
C
         IF(J.LE.0)GO TO 120
         DO 110 IR = K,N
            Z = A(IR,K)
            A(IR,K) = A(IR,ICM)
 110        A(IR,ICM) = Z
C
 120     DO 130 IC = K,N
            Z = PI*A(IRM,IC)
            A(IRM,IC) = A(K,IC)
 130        A(K,IC) = Z
C
         A(K,K) = J
C
         P = 0.D0
         IRS = K + 1
         DO 160 IR = IRS,N
            PI = -A(IR,K)
            IST = K+1
            DO 150 IC = IST,N
               A(IR,IC) = A(IR,IC)+PI*A(K,IC)
               Z = DABS(A(IR,IC))
               IF(Z.LE.P)GO TO 150
               P = Z
               IRM = IR
               ICM = IC
 150           CONTINUE
 160        B(IR) = B(IR)+PI*B(K)
 170     CONTINUE
C
 180  IF (N-1) 230,220,190
 190  IR = N
      IC = N+1
      DO 210 I = 2,N
         IR = IR - 1
         IC = IC - 1
         L = A(IR,IR)+.5D0
         Z = B(IR)
         LL = IR
         DO 200 K = IC,N
            LL = LL+1
 200        Z = Z-A(IR,K)*B(LL)
         K = IR+L
         B(IR) = B(K)
 210     B(K) = Z
 220  RETURN
C
C     ERROR RETURN
 230  IER = -1
      RETURN
 240  CONTINUE
      IER = -2
      RETURN
      END
      SUBROUTINE EVOLV (N,X,V,AMU,BE,H,FMIN,NIT,NR,II,IER,
     1    PRF,IDER,NU,NO,IDA,NFEV,NJEV,
     2    FP,C,VFU,BUF,XTP,WF,FPS,CS,S)
C
C THE SUBROUTINE  EVOLV  ATTEMPTS TO SOLVE A SYSTEM
C OF  N  NONLINEAR EQUATIONS
C     FI(X) = 0,  I = 1,...,N,  X = (X1,...,XN)        (1)
C THE SOLUTION IS SOUGHT AT A ZERO-VALUED MINIMUM OF THE WEIGHTED SUM
C     FW = W1*F1(X)**2 + ... + WN*FN(X)**2             (2)
C OF THE SQUARES OF THE LEFT-HAND SIDE FUNCTIONS OF SYSTEM (1).
C THE MINIMIZER IS SOUGHT AS THE ASYMPTOTIC VALUE (AS TIME  T  GOES TO
C INFINITY) OF THE SOLUTION TRAJECTORY  X(T) = (X1(T),...,XN(T))  OF
C THE ASSOCIATED DYNAMICAL SYSTEM DEFINED BY THE DIFFERENTIAL EQUATIONS
C     AMU*(DV/DT) = -BE*V -(DFW/DX)                    (3)
C WHERE  V = (V1,...,VN)  AND  (DV/DT)  ARE THE VECTORS OF THE
C VELOCITIES  VI(T) = DXI/DT  AND OF THE ACCELERATIONS  DVI/DT,
C I = 1,...,N,  (DFW/DX)  IS THE GRADIENT OF  FW  W.R.T.  X ,  AND
C AMU  AND  BE  ARE THE (SCALAR) MASS AND FRICTION COEFFICIENT.
C THE NUMERICAL INTEGRATION OF (3) IS PERFORMED WITH THE
C LAMBERT-SIGURDSSON LINEARLY IMPLICIT SCHEME
C WITH VARIABLE TIME INTEGRATION STEPLENGTH AND WITH
C THE POSSIBILITY OF REJECTING A STEP WHICH IS ESTIMATED
C NON SATISFACTORY.
C NOTE THAT THE WEIGHTING OF THE SQUARES OF THE FUNCTIONS IN  (2)
C IS EQUIVALENT TO RESCALING THE FUNCTION VALUES WITH FACTORS SQRT(W1),
C ..., SQRT(WN).
C
C DESCRIPTION OF PARAMETERS
C
C N  IS THE NUMBER OF NONLINEAR EQUATIONS.
C X,V  ARE VECTORS OF LENGTH  N  CONTAINING X-VARIABLES AND VELOCITIES.
C     (INITIAL VALUES ON INPUT, FINAL ON OUTPUT. INITIAL AND FINAL  X
C     ARE THE INITIAL AND FINAL ESTIMATES OF THE SOLUTION OF SYSTEM (1))
C AMU  IS THE VALUE OF THE  MASS COEFFICIENT.
C BE  IS THE FRICTION COEFFICIENT
C     (INITIAL VALUE ON INPUT, FINAL VALUE ON OUTPUT).
C H  IS THE TIME STEPLENGTH FOR THE NUMERICAL INTEGRATION OF SYSTEM (3)
C     (INITIAL VALUE ON INPUT, FINAL VALUE ON OUTPUT).
C FMIN  IS  A VECTOR OF DIMENSION  N  CONTAINING THE THRESHOLD
C     VALUES FOR THE FUNCTION VALUES  FI(X) , I = 1,...,N  (SEE BELOW).
C NIT  IS THE MAXIMUM ALLOWED NUMBER OF TIME INTEGRATION STEPS.
C NR  IS THE TOTAL NUMBER OF ATTEMPTED TIME INTEGRATION STEPS
C     (I.E. INCLUDING STEPS ATTEMPTED BUT REJECTED)
C II  IS THE NUMBER OF ACCEPTED STEPS.
C IER  IS AN OUTPUT INDICATOR OF THE STOPPING CONDITION, AS FOLLOWS
C     IER = 0   CONVERGENCE REACHED  AT STEP  II ,  I.E.
C               ABS(FI(X)).LT.FMIN(I) , I = 1,...,N
C     IER = 1   NUMBER OF ACCEPTED INTEGRATION STEPS  .EQ.NIT
C     IER = 2   ESTIMATED RATE OF CONVERGENCE TOO SLOW
C     IER = 3   ATTEMPT TO REDUCE THE TIME INTEGRATION STEPLENGTH
C               BELOW THE ALLOWED MINIMUM.
C     IER = 4   FATAL FAILURE IN COMPUTING THE DERIVATIVES.
C     IER = 5   FATAL FAILURE IN COMPUTING THE FUNCTION VALUES AT THE
C               INITIAL POINT.
C THE PARAMETERS  PRF , IDER , NU , NO , IDA
C ARE USED ONLY IN THE COMPUTATION OF THE NUMERICAL JACOBIAN
C ( SUBROUTINE  JCBN ), AND ARE DESCRIBED WITHIN THE SUBROUTINE  DAFNE .
C THE FOLLOWING 2 PARAMETERS ( NFEV , NJEV ) ARE OUTPUT COUNTERS
C     NFEV  OF FUNCTION EVALUATIONS (EXCLUDING THOSE USED
C           FOR JACOBIAN EVALUATION).
C     NJEV  OF JACOBIAN EVALUATIONS.
C FP,C  ARE WORKING STORAGE ARRAYS OF DIMENSION  N  AND  (N,N) .
C VFU  IS A VECTOR OF DIMENSION  N  USED FOR INTERNAL STORAGE
C     OF THE VALUES OF THE FUNCTIONS IN (1).
C BUF  AND  XTP  ARE WORKING STORAGE ARRAY OF LENGTH  N .
C WF  IS THE VECTOR OF THE WEIGHTS DESCRIBED IN SUBROUTINE  DAFNE .
C FPS , CS , S  ARE WORKING STORAGE ARRAYS OF DIMENSION  N , (N,N) , N
C     NEEDED BY SUBROUTINE  DLIN .
C
C SUBPROGRAMS CALLED  :  CHAH , FUNZ , DERIV , DLIN .
C
      DOUBLE PRECISION AMU,B,BE,BEMAX,BUF,C,CS,EMAX,EMIN
      DOUBLE PRECISION ET,F,FMIN,FP,FPS,GA,H,HMAX,HMIN
      DOUBLE PRECISION PRF,RAPMIN,S,V,VF,VFU,WF,X,XTP
      DOUBLE PRECISION FUNZ
      DIMENSION X(N),V(N),FMIN(N),FP(N),C(N,N),VFU(N),BUF(N)
      DIMENSION XTP(N),WF(N),FPS(N),CS(N,N),S(N)
      DIMENSION VF(30)
C
      DATA RAPMIN/1.0001D0/
C
C EXTREME VALUES FOR  H  AND  BE
C USE THE FIRST VALUES (-50,50,30) OR THE SECOND ONES (-25,25,17) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1.D120 (E.G. UNIVAC 1100), OR BELOW
C 1.D120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA HMIN,HMAX/1.D-50,1.D+50/,BEMAX/1.D30/
C     DATA HMIN,HMAX/1.D-25,1.D+25/,BEMAX/1.D17/
C
C NSR IS THE NUMBER OF STEPS FOR WHICH TOTAL ENERGY IS REMEMBERED
C NSR MUST BE .LE. DIMENSION OF VF
      DATA NSR/30/
C
      DATA ET/0.D0/
C
C INITIALIZE
      II = 0
      NR = 0
      F = FUNZ(N,X,VFU,NR,NFEV,WF,IER)
      IF (IER.LT.0) GO TO 190
      DO 10 ID = 1,NSR
         VF(ID) = F*DBLE(FLOAT(ID))**4
  10     CONTINUE
      H = DMAX1(H,HMIN)
C
C START NUMERICAL INTEGRATION LOOP  ( NIT  STEPS)
C
      DO 150 IT = 1,NIT
         II = IT
         DO 20 ID = 2,NSR
            IID = NSR+1-ID
            VF(IID+1) = VF(IID)
  20        CONTINUE
C
C        START OF TIME INTEGRATION STEP
C
  30     CONTINUE
         NR = NR+1
C
C        COMPUTE GRADIENT AND HESSIAN
         CALL DERIV (N,X,VFU,FP,C,WF,NR,NJEV,IER,
     1       PRF,IDER,NU,NO,IDA,BUF,XTP)
         IF (IER.LT.0) GO TO 200
C
         DO 50 L = 1,N
            DO 40 J = 1,N
               C(L,J) = C(L,J)*H
  40           CONTINUE
            C(L,L) = C(L,L)+BE+AMU/H
            FP(L) = -FP(L)*H+V(L)*AMU
  50        CONTINUE
C
C        SOLUTION OF LINEAR SYSTEM
         CALL DLIN (N,FP,C,FPS,CS,S,IER)
C
C        SAVE AND UPDATE X-VARIABLES AND VELOCITIES
         DO 60 L = 1,N
            C(1,L) = X(L)
            BUF(L) = VFU(L)
            B = V(L)
            V(L) = FP(L)/H
            X(L) = X(L)+FP(L)
            FP(L) = B
  60        CONTINUE
         IF (IER.LT.0) GO TO 100
C
C        COMPUTE AND SAVE TOTAL ENERGY
         F = FUNZ(N,X,VFU,NR,NFEV,WF,IER)
         IF (IER.LT.0) GO TO 110
         ET = F
         DO 70 J = 1,N
            ET = ET+AMU*0.5D0*V(J)**2
  70        CONTINUE
         VF(1) = ET
C
C
C        END OF TIME INTEGRATION STEP
C
C        TEST CONVERGENCE
         DO 80 J = 1,N
            IF (DABS(VFU(J)).GT.FMIN(J)) GO TO 90
  80        CONTINUE
         GO TO 160
  90     CONTINUE
C
 100     CONTINUE
 110     CONTINUE
C
C        COMPUTE TIME INTEGRATION STEPLENGTH
         CALL CHAH (NSR,VF,GA,IER,II,NR)
         H = H*GA
         IF (H.LT.HMIN) GO TO 170
         H = DMIN1(H,HMAX)
C
C        IF STEP REJECTED RESTORE PRECEDING CONDITIONS
         IF (IER.GE.0) GO TO 130
         DO 120 L = 1,N
            X(L) = C(1,L)
            V(L) = FP(L)
            VFU(L) = BUF(L)
 120        CONTINUE
         GO TO 30
C
 130     CONTINUE
C
C        TEST IF CONVERGENCE TOO SLOW
         EMAX = ET
         EMIN = ET
         DO 140 ID = 2,NSR
            EMAX = DMAX1(VF(ID),EMAX)
            EMIN = DMIN1(VF(ID),EMIN)
 140        CONTINUE
         IF (EMAX.LT.EMIN*RAPMIN) GO TO 180
C
C        COMPUTE FRICTION COEFFICIENT
         IF (II.GT.10) BE = DMIN1(BEMAX,BE+BE)
C
 150     CONTINUE
C
C END OF NUMERICAL INTEGRATION LOOP  ( NIT  STEPS)
C
C EXIT CONDITIONS
C
C MAX  NUMBER OF TIME INTEGRATION STEPS
      IER = 1
      RETURN
C
C CONVERGENCE
 160  CONTINUE
      IER = 0
      RETURN
C
C ATTEMPT TO REDUCE THE TIME INTEGRATION STEP BELOW THE ALLOWED MINIMUM
 170  CONTINUE
      IER = 3
      RETURN
C
C ESTIMATED RATE OF CONVERGENCE TOO SLOW
 180  CONTINUE
      IER = 2
      RETURN
C
C FATAL FAILURE IN COMPUTING FUNCTION VALUES AT INITIAL POINT.
 190  CONTINUE
      IER = 5
      RETURN
C
C FATAL FAILURE IN COMPUTING DERIVATIVES.
 200  CONTINUE
      IER = 4
      RETURN
C
      END
      DOUBLE PRECISION FUNCTION FUNZ (N,X,VFU,NR,NFEV,WF,IER)
C
C THE FUNCTION  FUNZ  RETURNS AS ITS
C FUNCTION VALUE THE WEIGHTED SUM
C     FW = W1*F1(X)**2 + ... + WN*FN(X)**2
C OF THE SQUARES OF THE LEFT-HAND SIDE FUNCTIONS OF SYSTEM (1)
C (SEE SUBROUTINE  EVOLV ), WHICH AMOUNTS TO RESCALE THE
C FUNCTION VALUES WITH FACTORS  SQRT(W1),..., SQRT(WN),
C AND RETURNS IN THE VECTOR  VFU  OF DIMENSION  N  THE VALUES
C OF THE (NON-RESCALED) FUNCTIONS IN (1), COMPUTED AT  X
C BY CALLING THE SUBROUTINE  VECF .
C NR  IS THE NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C NFEV  IS THE COUNTER OF FUNCTION EVALUATIONS (EXCLUDING
C     THOSE USED FOR JACOBIAN EVALUATION.
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     INDICATES A FAILURE IN COMPUTING A FUNCTION VALUE).
C THE WEIGHTED SUM OF THE SQUARES IS COMPUTED BY MEANS OF
C THE WEIGHTS  W1,...,WN  CONTAINED IN THE N-VECTOR  WF .
C
C SUBPROGRAMS CALLED  :  VECF
C   ( VECF  MUST BE SUPPLIED BY THE USER, SEE SUBROUTINE  DAFNE ).
C
      DOUBLE PRECISION F,VFU,WF,X
      DIMENSION X(N),VFU(N),WF(N)
C
C UPDATE FUNCTION EVALUATIONS COUNTER
      NFEV = NFEV+1
C
C COMPUTE FUNCTION VALUES
      CALL VECF (N,X,VFU,NR,IER)
      IF (IER.LT.0) RETURN
C
C COMPUTE THE WEIGHTED SUM OF THE SQUARES
      F = 0.D0
      DO 10 I = 1,N
         F = F+WF(I)*VFU(I)**2
  10     CONTINUE
C
      FUNZ = F
      RETURN
      END
      SUBROUTINE JCBN (N,X,VFU,C,FP,BUF,XTP,PRF,IDER,NU,NO,IDA,
     1    NR,NJEV,IER)
C
C THE SUBROUTINE  JCBN  COMPUTES THE JACOBIAN MATRIX
C OF THE SYSTEM (1) (SEE SUBROUTINE  EVOLV ).
C
C DESCRIPTION OF PARAMETERS
C
C THE PARAMETERS  N , X , VFU , FP , BUF  ARE DESCRIBED WITHIN THE
C SUBROUTINE  EVOLV .
C
C C  IS AN OUTPUT N BY N ARRAY CONTAINING THE COMPUTED JACOBIAN
C
C XTP  IS AN N-VECTOR INITIALIZED BY SUBROUTINE  CREXTP ,
C     UPDATED BY SUBROUTINE  NEWXTP , AND USED HERE FOR
C     COMPUTING THE FINITE-DIFFERENCES JACOBIAN.
C
C PRF  IS AN (INPUT) ROUGH ESTIMATE OF THE AVERAGE RELATIVE
C     ERROR AFFECTING THE FUNCTION VALUES SUPPLIED BY SUBROUTINE  VECF .
C     IF THIS ERROR IS ESTIMATED TO BE NOT SEVERAL ORDERS OF MAGNITUDE
C     GREATER THAN MACHINE PRECISION, THE USER MAY SIMPLY PUT  PRF = 0 .
C
C IDER  IS AN INPUT CONTROL VARIABLE TO REQUEST THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDER.EQ.0   ANALYTICAL COMPUTATION BY SUBROUTINE  JCBNA
C     IDER.EQ.1   COMPUTATION BY FORWARD FINITE DIFFERENCES
C     IDER.EQ.2   COMPUTATION BY CENTRAL FINITE DIFFERENCES
C
C NU  AND  NO  ARE NUMBERS DEFINING THE (POSSIBLY) BANDED
C STRUCTURE OF THE JACOBIAN MATRIX, USED FOR FINITE-DIFFERENCES
C JACOBIAN COMPUTATION
C     NU IS THE N. OF NON-ZERO JACOBIAN DIAGONALS UNDER MAIN DIAGONAL
C     NO IS THE N. OF NON-ZERO JACOBIAN DIAGONALS ABOVE MAIN DIAGONAL
C     (FOR A NON-BANDED MATRIX  NU = NO = N-1 )
C
C IDA  IS AN (OUTPUT) INDICATOR OF THE CALCULATION MODE
C     FOR THE DERIVATIVES, AS FOLLOWS
C     IDA.EQ.1   THE DERIVATIVES WERE COMPUTED ANALYTICALLY
C                (AS REQUESTED)
C     IDA.EQ.0   THE DERIVATIVES WERE COMPUTED BY FINITE DIFFERENCES.
C                (CENTRAL FINITE DIFFERENCES IF REQUESTED BY  IDER = 2,
C                FORWARD FINITE DIFFERENCES IF REQUESTED BY  IDER = 1,
C                OR IF REQUESTED ANALYTICAL DERIVATIVES WERE NOT CODED
C                IN SUBROUTINE  JCBNA ).
C
C NR  IS THE NUMBER OF THE TIME INTEGRATION STEP IN PROGRESS
C     (COUNTING ALSO THE REJECTED STEPS).
C NJEV  IS THE COUNTER OF THE JACOBIAN EVALUATIONS
C IER IS AN OUTPUT ERROR INDICATOR (A NEGATIVE VALUE
C     INDICATES A FAILURE IN COMPUTING A FUNCTION VALUE).
C
C SUBPROGRAMS CALLED  :  PRMACH , CREXTP , NEWXTP , JCBNA
C ( JCBNA  MUST BE SUPPLIED BY THE USER, SEE SUBROUTINE  DAFNE ).
C
      DOUBLE PRECISION BUF,C,E,EE,XTP,FP
      DOUBLE PRECISION PR,PRF,PRM,PRR,SPRR
      DOUBLE PRECISION VFU,X,XX
      DIMENSION X(N),FP(N),C(N,N),BUF(N),VFU(N),XTP(N)
      DATA PRM/0.D0/
C
C UPDATE JACOBIAN EVALUATIONS COUNTER
      NJEV = NJEV+1
C
      IDA = 0
      IER = 0
C
C TRY TO COMPUTE ANALYTICAL JACOBIAN, IF REQUIRED
      IF (IDER.EQ.0) CALL JCBNA (N,X,C,NR,IDA,IER)
C
C ERROR EXIT
      IF (IER.LT.0) RETURN
C
C EXIT DUE TO SUCCESSFUL COMPUTATION OF THE
C ANALYTICAL JACOBIAN
      IF (IDA.EQ.1) RETURN
C
C COMPUTE FINITE-DIFFERENCE JACOBIAN (IF REQUIRED,
C OR DUE TO AN UNSUCCESSFUL ATTEMPT TO COMPUTE ANALYTICAL
C JACOBIAN).
C
      IDERM = IDER
      IF (IDER.GT.0) GO TO 10
C
C COMPUTE FORWARD FINITE-DIFFERENCE JACOBIAN IF
C ANALYTICAL JACOBIAN WAS NOT AVAILABLE.
      IDER = 1
C
  10  CONTINUE
C
C INITIALIZATION (PERFORMED ONLY IN THE FIRST
C TIME INTEGRATION STEP).
C
      IF (NR.GT.1) GO TO 20
C
C ESTIMATE MACHINE PRECISION  PRM
      IF (PRM.LE.0.D0) CALL PRMACH (PRM)
C
C ESTIMATE ACTUAL WORKING PRECISION  PR
      PR = DMAX1(PRF,PRM)
      PRR = PR
      IF (IDER.EQ.2) PRR = PR**(2.D0/3.D0)
      SPRR = DSQRT(PRR)
C
C INITIALIZE VECTOR  XTP
      CALL CREXTP (N,X,SPRR,XTP)
C
  20  CONTINUE
C END OF INITIALIZATION
C
      NRP = NR+1
C
C TEST APPLICABILITY OF BANDED JACOBIAN ALGORITHM
      ND = NO+NU+1
      IF (ND.LT.N) GO TO 70
C
C NON-BANDED JACOBIAN ALGORITHM
      DO 60 I = 1,N
C
C        COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
         E = SPRR*DMAX1(DABS(X(I)),XTP(I))
         E = DSIGN(E,X(I))
C
         XX = X(I)
         X(I) = X(I)+E
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,X,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         EE = 1.D0/E
C
C        COMPUTE FORWARD FINITE-DIFFERENCE DERIVATIVES
         DO 30 J = 1,N
            C(J,I) = (BUF(J)-VFU(J))*EE
  30        CONTINUE
C
         IF (IDER.EQ.1) GO TO 50
C
C        COMPUTATION FOR CENTRAL DIFFERENCE DERIVATIVES
         X(I) = XX-E
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,X,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         DO 40 J = 1,N
            C(J,I) = .5D0*(C(J,I)-(BUF(J)-VFU(J))*EE)
  40        CONTINUE
C
  50     CONTINUE
         X(I) = XX
C
  60     CONTINUE
C END OF COMPUTATIONS FOR CENTRAL DIFFERENCE DERIVATIVES
C
      GO TO 160
C
C ALGORITHM FOR BANDED JACOBIAN
  70  CONTINUE
C
      DO 80 I = 1,N
         FP(I) = X(I)
  80     CONTINUE
C
      DO 150 I = 1,ND
C
         DO 90 J = I,N,ND
C
C           COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
            E = SPRR*DMAX1(DABS(X(J)),XTP(J))
            E = DSIGN(E,X(J))
C
            FP(J) = X(J)+E
  90        CONTINUE
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,FP,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         DO 110 J = I,N,ND
C
C           COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
            E = SPRR*DMAX1(DABS(X(J)),XTP(J))
            E = DSIGN(E,X(J))
C
            FP(J) = X(J)-E*DBLE(FLOAT(IDER-1))
            EE = 1.D0/E
            K1 = MAX0(1,J-NO)
            K2 = MIN0(N,J+NU)
C
C           FORWARD FINITE DIFFERENCES
            DO 100 K = K1,K2
               C(K,J) = (BUF(K)-VFU(K))*EE
 100           CONTINUE
C
 110        CONTINUE
C
         IF (IDER.EQ.1) GO TO 140
C
C        CENTRAL FINITE DIFFERENCES
C
C        COMPUTE FUNCTION VALUES
         CALL VECF (N,FP,BUF,NRP,IER)
         IF (IER.LT.0) RETURN
C
         DO 130 J = 1,N,ND
C
C           COMPUTE STEPLENGTH FOR FINITE-DIFFERENCE DERIVATIVES
            E = SPRR*DMAX1(DABS(X(J)),XTP(J))
            E = DSIGN(E,X(J))
C
            FP(J) = X(J)
            EE = 1.D0/E
            K1 = MAX0(1,J-NO)
            K2 = MIN0(N,J+NU)
C
            DO 120 K = K1,K2
               C(K,J) = .5D0*(C(K,J)-(BUF(K)-VFU(K))*EE)
 120           CONTINUE
 130        CONTINUE
 140     CONTINUE
 150     CONTINUE
 160  CONTINUE
C
C UPDATE VECTOR  XTP
      CALL NEWXTP (N,X,VFU,C,XTP)
C
      IDER = IDERM
      RETURN
      END
      SUBROUTINE MOVM (N,FP,C,FPS,CS)
C
C THE SUBROUTINE  MOVM  SAVES THE ARRAYS  FP  AND  C
C (OF DIMENSION  N  AND  (N,N) ) IN THE ARRAYS  FPS  AND  CS .
C
      DOUBLE PRECISION FP,C,FPS,CS
      DIMENSION FP(N),C(N,N),FPS(N),CS(N,N)
C
      DO 20 I = 1,N
         FPS(I) = FP(I)
         DO 10 J = 1,N
            CS(I,J) = C(I,J)
  10        CONTINUE
  20     CONTINUE
      RETURN
      END
      SUBROUTINE NEWF (N,X,WF,PRF,IDER,NU,NO,IDA,
     1    XTP,VFU,C,FP,BUF,NRP,NJEV,IER)
C
C THE SUBROUTINE NEWF COMPUTES THE NEW WEIGHTS FOR FUNCTION RESCALING
C THE PARAMETERS ARE DESCRIBED WITHIN THE SUBROUTINE  DAFNE .
C
      DOUBLE PRECISION X,VFU,FP,BUF,C,WF,PRF,XTP
      DOUBLE PRECISION RWMAX,RWMIN
      DOUBLE PRECISION BMAX,BMIN
      DOUBLE PRECISION SB,B,RW
      DIMENSION X(N),VFU(N),FP(N),BUF(N),C(N,N),WF(N),XTP(N)
C
C RWMAX  AND  RWMIN  ARE THE MAX AND MIN VALUES OF THE
C MULTIPLICATIVE FACTOR USED TO UPDATE THE WEIGHTS
      DATA RWMAX,RWMIN/1.D+3,1.D-3/
C
C BMAX  AND  BMIN  DEFINE THE ALLOWED RANGE FOR THE SQUARED NORM  B .
C USE THE FIRST VALUES (40,-40) OR THE SECOND ONES (30,-30) IF THE
C MACHINE DYNAMIC RANGE IS ABOVE 1.D120 (E.G. UNIVAC 1100), OR BELOW
C 1.D120 (E.G. AMDAHL 470, IBM 360/370, DEC VAX)
      DATA BMAX,BMIN/1.D+40,1.D-40/
C     DATA BMAX,BMIN/1.D+30,1.D-30/
C
C COMPUTE FUNCTION VALUES
      CALL VECF (N,X,VFU,NRP,IER)
      IF (IER.LT.0) RETURN
C
C COMPUTE JACOBIAN
      CALL JCBN (N,X,VFU,C,FP,BUF,XTP,PRF,IDER,NU,NO,IDA,
     1    NRP,NJEV,IER)
C
      SB = 1.D0
      DO 20 I = 1,N
C
C        COMPUTE SQUARED NORM OF GRADIENT OF FUNCTION  I
         B = 0.D0
         DO 10 J = 1,N
            B = B+C(I,J)**2
  10        CONTINUE
C
         B = DMAX1(B,BMIN)
         B = DMIN1(B,BMAX)
         BUF(I) = WF(I)*B
         SB = SB*BUF(I)**(1.D0/DBLE(FLOAT(N)))
C
  20     CONTINUE
C
C COMPUTE WEIGHTS
      DO 30 I = 1,N
         RW = SB/BUF(I)
         RW = DMIN1(RW,RWMAX)
         RW = DMAX1(RW,RWMIN)
         WF(I) = WF(I)*RW
  30     CONTINUE
C
      RETURN
      END
      SUBROUTINE NEWXTP (N,X,VFU,C,XTP)
C
C THE SUBROUTINE  NEWXTP  UPDATES THE ARRAY  XTP
C USED BY SUBROUTINE  JCBN  TO COMPUTE THE STEPLENGTHS USED
C FOR COMPUTING THE FINITE-DIFFERENCES JACOBIAN.
C
C THE PARAMETERS ARE DEFINED IN THE SUBROUTINE  EVOLV .
C
      DOUBLE PRECISION C,ECIJ,XTP,RM,RP,VFI,VFU,X
      DIMENSION X(N),VFU(N),C(N,N),XTP(N)
C
C RP  AND  RM  ARE THE VALUES OF THE MULTIPLICATIVE
C     FACTOR USED TO UPDATE THE VECTOR  XTP .
      DATA RP,RM/10.D0,.1D0/
C
      DO 20 I = 1,N
         IM = 0
         VFI = DABS(VFU(I))
         DO 10 J = 1,N
            ECIJ = DABS(C(I,J)*XTP(J))
            IC = 1
            IF (VFI.LT.ECIJ*RP) IC = 2
            IF (VFI.LT.ECIJ*RM) IC = 3
            IM = MAX0(IM,IC)
  10        CONTINUE
         IF (IM.EQ.1) GO TO 50
         IF (IM.EQ.2) XTP(I) = -DABS(XTP(I))
         IF (IM.EQ.3) XTP(I) = DABS(XTP(I))
  20     CONTINUE
      DO 40 J = 1,N
         IM = 0
         DO 30 I = 1,N
            VFI = DABS(VFU(I))
            ECIJ = DABS(C(I,J)*XTP(J))
            IC = 1
            IF (VFI.LT.ECIJ*RP) IC = 2
            IF (VFI.LT.ECIJ*RM) IC = 3
            IF (IC.EQ.2.AND.XTP(I).LT.0.D0) IC = 4
            IM = MAX0(IM,IC)
  30        CONTINUE
         IF (IM.EQ.1) XTP(J) = XTP(J)*RP
         IF (IM.EQ.3) XTP(J) = XTP(J)*RM
  40     CONTINUE
      GO TO 70
  50  CONTINUE
      DO 60 J = 1,N
         XTP(J) = XTP(J)*RP
  60     CONTINUE
  70  CONTINUE
      DO 80 J = 1,N
         XTP(J) = DMAX1(DABS(XTP(J)),DABS(X(J)))
  80     CONTINUE
      RETURN
      END
      SUBROUTINE PRMACH (PRM)
C
C THE SUBROUTINE  PRMACH  ESTIMATES THE MACHINE PRECISION  PRM
C (SOFTWARE PRECISION, ESTIMATED ON SOME INTRINSIC
C FORTRAN FUNCTIONS).
C
      DOUBLE PRECISION A,E,P,PRM
C
      P = 0.D0
      DO 10 I = 1,64
         A = DBLE(FLOAT(I))
         E = 0.D0
         E = E+(DSQRT(A*A)-A)**2
         E = E+(DSQRT(A)**2-A)**2
         E = E+(DEXP(DLOG(A))-A)**2
         E = E+(DLOG(DEXP(A))-A)**2
         E = E+(A*DSIN(A)**2+A*DCOS(A)**2-A)**2
         P = P+E/(5.D0*A**2)
  10     CONTINUE
      P = DSQRT(P/64.D0)
      PRM = P
C
      RETURN
      END
      SUBROUTINE QMS (N,FP,C,FPS,CS,D)
C
C THE SUBROUTINE  QMS  REGULARIZES THE LINEAR SYSTEM
C WITH  N BY N  COEFFICIENT MATRIX  C  AND RIGHT-HAND-SIDE
C N-VECTOR  FP  WITH REGULARIZATION PARAMETER  D
C OBTAINING THE LINEAR SYSTEM  CS , FPS .
C
      DOUBLE PRECISION FP,C,FPS,CS,D
      DIMENSION FP(N),C(N,N),FPS(N),CS(N,N)
C
      DO 30 I = 1,N
         FPS(I) = 0.D0
         DO 20 J = 1,N
            FPS(I) = FPS(I)+C(J,I)*FP(J)
            CS(I,J) = 0.D0
            DO 10 K = 1,N
               CS(I,J) = CS(I,J)+C(K,I)*C(K,J)
  10           CONTINUE
  20        CONTINUE
         CS(I,I) = CS(I,I)+D
  30     CONTINUE
C
      RETURN
      END
         1         ROSENBROCK FUNCTION.
    1    2                        . N
    7    3                        . NRIP
    8    1                        . NRS1
    9    3                        . NRS2
    0
         2         POWELL SINGULAR FUNCTION.
    1    4
    0
         3         POWELL BADLY SCALED FUNCTION.
    1    2
    0
         4         WOOD FUNCTION.
    1    4
    0
         5         HELICAL VALLEY FUNCTION.
    1    3
    0
         6         WATSON FUNCTION.
    1    6
    0
         6         WATSON FUNCTION.
    1    9
    0
         7         CHEBYQUAD FUNCTION.
    1    5
    0
         7         CHEBYQUAD FUNCTION.
    1    6
    0
         7         CHEBYQUAD FUNCTION.
    1    7
    0
         7         CHEBYQUAD FUNCTION.
    1    8
    0
         7         CHEBYQUAD FUNCTION.
    1    9
    0
         8         BROWN ALMOST LINEAR FUNCTION.
    1   10
    0
         8         BROWN ALMOST LINEAR FUNCTION.
    1   30
    0
         8         BROWN ALMOST LINEAR FUNCTION.
    1   40
    0
         9         DISCRETE BOUNDARY VALUE FUNCTION.
    1   10
    0
        10         DISCRETE INTEGRAL EQUATION FUNCTION.
    1    1
    0
        10         DISCRETE INTEGRAL EQUATION FUNCTION.
    1   10
    0
        11         TRIGONOMETRIC FUNCTION.
    1   10
    0
        12         VARIABLY DIMENSIONED FUNCTION.
    1   10
    0
        13         BROYDEN TRIDIAGONAL FUNCTION.
    1   10
    0
        14         BROYDEN BANDED FUNCTION.
    1   10
    0
         0