      SUBROUTINE CMTRIX(N,M,RJ,F,RJTJ,RJTF)                             
CDC      ALLOCLEVEL RJ,F,RJTJ,RJTF                                      
      INTEGER I,J,K,M,N                                                 
      REAL*8 RJ(M,N),F(M),RJTJ(N,N),RJTF(N)                             
C                                                                       
C       THIS ROUTINE FORMS THE MATRIX                                   
C     (JACOBIAN TRANSPOSE) * JACOBIAN  AND THE VECTOR                   
C     (JACOBIAN TRANSPOSE) * F.                                         
C                                                                       
C     ON INPUT                                                          
C                                                                       
C        N  IS THE DIMENSION OF THE VECTOR  RJTF, THE COLUMN            
C           DIMENSION OF THE ARRAY  RJ  AND THE DIMENSION OF THE SQUARE 
C           ARRAY  RJTJ.                                                
C                                                                       
C        M  IS THE DIMENSION OF THE VECTOR  F  AND THE ROW DIMENSION    
C           OF THE ARRAY  RJ.                                           
C                                                                       
C        RJ  CONTAINS THE ELEMENTS OF THE JACOBIAN MATRIX.              
C                                                                       
C        F  CONTAINS THE ELEMENTS OF THE FUNCTION VECTOR.               
C                                                                       
C     ON OUTPUT                                                         
C                                                                       
C        RJTJ  CONTAINS THE LOWER TRIANGLE OF THE MATRIX                
C           (JACOBIAN TRANSPOSE) * JACOBIAN  STORED BY ROWS.            
C                                                                       
C        RJTF  CONTAINS THE ELEMENTS OF THE VECTOR                      
C           (JACOBIAN TRANSPOSE) * F.                                   
C                                                                       
  100 DO 200 I = 1, N                                                   
      DO 200 J = 1, I                                                   
      RJTJ(I,J) = 0.0D0                                                 
C                                                                       
      DO 200 K = 1, M                                                   
      RJTJ(I,J) = RJTJ(I,J) + RJ(K,I) * RJ(K,J)                         
  200 CONTINUE                                                          
C                                                                       
      DO 300 I = 1, N                                                   
      RJTF(I) = 0.0D0                                                   
C                                                                       
      DO 300 K = 1, M                                                   
      RJTF(I) = RJTF(I) + F(K) * RJ(K,I)                                
  300 CONTINUE                                                          
C                                                                       
      RETURN                                                            
C     ----- LAST CARD OF SUBROUTINE  CMTRIX  -----                      
      END  
      
      FUNCTION FCUBIC(F0,G0,F1,G1,XLAM)                                 
C     CUBIC INTERPOLATION (ACTON VARIANT) FOR MINIMIZERS                
C     GIVEN FUNCTION (F0,F1) AND DERIVATIVE (G0,G1) AT TWO POINTS       
C     WITH XLAM=X2-X1 (X2.GT.X1) FITS A CUBIC TO THES 4 ITEMS OF        
C     INFORMATION AND RETURNS MINIMUM PREDICTED BY THIS CUBIC           
C     IF WSQ.LT.0 CUBIC HAS NO MINIMUM AND QUADRATIC INTERPOLATION      
C     IS USED INSTEAD. IN EITHER CASE FCUBIC=X(MIN)-X1                  
C     NO TESTS ARE MADE TO ENSURE THAT VALUE RETURNED IS SENSIBLE       
C     I.E. 0.LT.FCUBIC.LT.XLAM                                          
C     THIS SHOULD BE DONE IN CALLING PROGRAM                            
C     MHM JANUARY 4 1976                                                
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
C                                                                       
      F=(F1-F0)/XLAM                                                    
      G=G1-G0                                                           
      C=G-2.0D0*(F-G0)                                                      
      D=G-3.0D0*C                                                           
      WSQ=D*D-12.0D0*C*G0                                                   
      IF(WSQ.LT.0.0D0) GO TO 1                                              
      W=DSQRT(WSQ)                                                      
      Y=-2.0D0*G0/(D+W)                                                     
      GO TO 2                                                           
    1 B=G1-F                                                            
      Y=0.5D0*(1.0D0-F/B)                                                      
    2 FCUBIC=Y*XLAM                                                     
      RETURN                                                            
      END                                                               
C                                                                       
C     ------------------------------------------------------------------
C                                                                       
      SUBROUTINE LMCHOL(NVAR,NFUNS,X,MODE,MAXF,S,LMCFUN,TOL,TOL2,FVECT,
     * RJACOB,RJINV,IERR,NFCALL,NJCALL,SCRV1,SQSCAL,DJTJ,VJTF,SCRV2,
     * HISTRY,MONIT,IDEB )                                             
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,S,FVECT,RJACOB,RJINV,SCRV1,SQSCAL,DJTJ,VJTF,SCRV2 
C                                                                       
CK      INTEGER I,ICONV,IERR,MAXF,MODE,NFCALL,NFUNS,NJCALL,NVAR           
CK      REAL*8 X(NVAR),S(NVAR),FVECT(NFUNS),RJACOB(NFUNS,NVAR),RJINV(NVAR
CK     *,NVAR),SCRV1(NVAR),SQSCAL(NVAR),DJTJ(NVAR),VJTF(NVAR),SCRV2(NFUNS)
      DIMENSION X(NVAR),S(NVAR),FVECT(NFUNS),RJACOB(NFUNS,NVAR),
     * RJINV(NVAR,NVAR),SCRV1(NVAR),SQSCAL(NVAR),DJTJ(NVAR),VJTF(NVAR),
     * SCRV2(NFUNS)
CK      REAL*8 TOL,HISTRY                                                 
CK      REAL*8 ACTRED,DABS,DIRDRV,DSQRT,FSQ,FSQP,PRERED,QCPARM,QPARM,Z    
C                                                                       
C     LEAST SQUARES MINIMIZATION WITH ANALYTIC DERIVATIVES              
C                                                                       
C     THIS ROUTINE HAS BEEN SLIGHTLY MODIFIED BY S. PIEPER:             
C       1)  LMCFUN AND EVLJAC HAVE BEEN COMBINED                        
C       2)  TOL2 PARAMETER ADDED                                        
C       3)  MONIT AND IDEB ADDED                                        
C       4)  IMPLICIT STATEMENT ADDED.                                   
C       5)  A STOP SEARCH BASED ON SMALL CHANGES IN THE SUM OF SQUARES  
C           WAS ADDED.                                                  
C       6)  RJINV IS SET TO AN ESTIMATE OF THE COVARIANCE MATRIX AT     
C           THE END OF THE SEARCH.                                      
C       7)  THE TEST FOR CONVERGENCE DUE TO SMALL STEPS IS MADE         
C           BEFORE EACH LMCFUN CALL.  THIS SAVES AN UNNECESSARY         
C           LMCFUN CALL AFTER CONVERGENCE HAS BEEN ACHIEVED.            
C           ALSO THE FINAL JACOBIAN IS ALWAYS REQUESTED AT THE END      
C           OF AN ITERATION.                                            
C       8)  IF LMCFUN RETURNS  IERR = 1,  A STEP OUT-OF-BOUNDS IS       
C           IS ASSUMED; THE MARQUARDT PARAMETER IS MULTIPLIED BY 4      
C           AND A NEW STEP IS TRIED.                                    
C                                                                       
C       THIS ROUTINE IS A MODULARIZED MODIFICATION OF THE COMPUTER      
C     ALGORITHM  VA07A  BY R. FLETCHER IN HARWELL REPORT  AERE - R.6799.
C                                                                       
C       THIS SUBROUTINE MINIMIZES A SUM OF SQUARES  FSQ(X)  OF          
C     FUNCTIONS  F(J)(X(1),X(2),...,X(NVAR)), J = 1,2,...,NFUNS, USING  
C     A MODIFICATION OF THE LEVENBERG-MARQUARDT ALGORITHM. THE USER     
C     MUST PROVIDE A SUBROUTINE TO CALCULATE THE FUNCTIONS              
C     F(J)  AND THE ELEMENTS OF THE JACOBIAN                            
C     (FIRST PARTIAL DERIVATIVE) MATRIX.                                
C                                                                       
C     TERMINOLOGY:                                                      
C                                                                       
C         QUANTITY TO BE MINIMIZED:                                     
C            FSQ  =  SUM( 1 =< J =< NFUNS )  F(J)**2                    
C                                                                       
C         JACOBIAN:                                                     
C            JACOBIAN(J,N)  =  D(F(J)) / D(X(N)) ,                      
C              1 =< J =< NFUNS ,                                        
C              1 =< N =< NVAR .                                         
C                                                                       
C     ******************************************************************
C       THIS ROUTINE INVOKES THE  MINPACK  PACKAGE MODULES  CMTRIX,     
C     QMTRIX, UPQPAR  AND  SIP  AND THE USER SUPPLIED SUBROUTINES       
C     LMCFUN  AND  MONIT. THE  ANL  LIBRARY ROUTINE  ANL F152S, SYMINV  
C     IS ALSO INVOKED.                                                  
C     ******************************************************************
C                                                                       
C     ON INPUT                                                          
C                                                                       
C        NVAR  IS THE NUMBER OF VARIABLES. IT IS ALSO THE DIMENSION OF  
C           THE VECTORS  X, S, SCRV1, SQSCQL, DJTJ  AND  VJTF  THE      
C           SQUARE ARRAY  RJINV  AND IS THE COLUMN DIMENSION OF THE     
C           ARRAY  RJACOB.                                              
C                                                                       
C        NFUNS  IS THE NUMBER OF FUNCTIONS  F(J). IT IS ALSO THE        
C           DIMENSION OF THE VECTOR  FVECT  AND  SCRV2  AND IS THE ROW  
C           DIMENSION OF THE ARRAY  RJACOB. A DIAGNOSTIC IS RETURNED FOR
C           NFUNS .LT. NVAR.                                            
C                                                                       
C        X  CONTAINS AN ESTIMATE OF THE SOLUTION VECTOR                 
C           (X(1),X(2),...,X(NVAR)).                                    
C                                                                       
C        MODE  IS A PARAMETER SET EQUAL TO  K  WHICH GOVERNS THE        
C           METHOD OF SCALING THE VARIABLES THROUGH THE INPUT PARAMETER 
C           S. THE OPTIONS ARE                                          
C                                                                       
C           K = 1 - THE  S(I)  ARE SET BY  LMCHOL  TO THE DIAGONAL      
C                   ELEMENTS OF  (JACOBIAN TRANSPOSE) * JACOBIAN, IF    
C                   POSITIVE, AND TO ONE OTHERWISE,                     
C           K = 2 - THE  S(I)  ARE SET BY  LMCHOL  TO ONE FOR ALL  I,   
C           K = 3 - THE  S(I)  ARE SET BY THE USER THROUGH THE          
C                   PARAMETER  S,                                       
C                                                                       
C           IF  MODE  IS SET EQUAL TO ANY OTHER VALUE IT IS RESET TO    
C           ONE BY  LMCHOL.                                             
C                                                                       
C        MAXF  IS THE LIMIT ON THE NUMBER OF CALLS TO  LMCFUN.          
C                                                                       
C        S  CONTAINS THE MULTIPLIERS WHICH CONTROL SCALING OF THE       
C           VARIABLES IF  MODE  IS SET TO 3. OTHERWISE IT IS            
C           INITIALIZED INTERNALLY.                                     
C                                                                       
C        LMCFUN  IS A USER SUPPLIED SUBROUTINE TO CALCULATE FUNCTIONS   
C           F(J)  IN  FVECT  AT THE ESTIMATE  X(I),I = 1,2,...,NVAR. THE
C           USER MAY FORCE AN EXIT FROM  LMCHOL  BY SETTING  IERR  TO A 
C           NEGATIVE INTEGER IN  LMCFUN. LMCFUN, MUST BE DECLARED IN AN 
C           EXTERNAL  STATEMENT IN THE USERS CALLING PROGRAM, AND       
C           SHOULD BE WRITTEN AS FOLLOWS                                
C                                                                       
C           SUBROUTINE LMCFUN ( ITYPE, NVAR,NFUNS,X,FVECT, RJACOB,IERR) 
C           REAL*8  X(NVAR),FVECT(NFUNS)                                
C           REAL*8  RJACOB(NFUNS,NVAR)                                  
C           INTEGER  NVAR,NFUNS,IERR                                    
C           -----------------------------                               
C                                                                       
C       STATEMENTS TO EVALUATE                                          
C         1)  FVECT = F(J)       IF ITYPE = 1 OR 3.                     
C         2)  RJACOB = JACOBIAN   IF ITYPE = 2 OR 3.                    
C         3)  SET IERR IF DESIRED TO STOP SEARCH.                       
C     OR                                                                
C         1)  SET IERR = +1 TO INDICATE A STEP OUT OF BOUNDS.           
C                                                                       
C     IF THE USER FINDS IT MORE CONVIENENT (OR OF LITTLE EXTRA COST),   
C     TO COMPUTE BOTH QUANTITIES WHEN ITYPE IS 1 OR 2, HE MAY DO SO AND 
C     SHOULD THEN SET ITYPE TO 3.  IN MANY CASES, COMPUTING THE GRADIENT
C     WHEN ONLY THE FUNCTION IS REQUESTED WILL SAVE A SUBSEQUENT REQUEST
C     FOR THE GRADIENT AT THE SAME PLACE ( HOWEVER, SOMETIMES THE       
C     GRADIENT WILL NEVER BE USED).                                     
C                                                                       
C           -----------------------------                               
C           RETURN                                                      
C           END                                                         
C                                                                       
C                                                                       
C        TOL & TOL2  IS THE ACCURACY REQUIRED IN THE SOLUTION, I.E. A   
C           NORMAL RETURN FROM THE ROUTINE OCCURS WHEN THE DIFFERENCES  
C           BETWEEN THE COMPONENTS OF TWO SUCCESSIVE ESTIMATES OF THE   
C           SOLUTION ARE NOT GREATER THAN  MAX(TOL*ABS(X(I)),TOL2)  FOR 
C           ALL  I. NOTE THAT A NORMAL RETURN WILL ALSO OCCUR IF THE    
C           PREDICTED REDUCTION IN THE SUM OF SQUARES IS .LE. 0.        
C          NORMAL RETURN IS ALSO MADE IF                                
C             FSQ(I-2) - FSQ(I)  <  TOL*FSQ(I)                          
C          OR IF                                                        
C             FSQ(I) < TOL2 .                                           
C          THIS CRITERION IS DESINGED TO PREVENT MINIMIZATION BASED ON  
C          INSIGNIFICANT (OR POSSIBLY NOISY) CHANGES IN THE SUM OF THE  
C          SQUARES.  NOTE THE COMPARASON OF STEP I WITH THAT FOR        
C          I-2 TO ALLOW A MOMENTARY PAUSE IN THE FUNCTION IMPROVEMENT.  
C                                                                       
C        HISTRY  IS A DUMMY PARAMETER.                                  
C                                                                       
C     MONIT IS A SUBROUTINE THAT WILL BE CALLED AT THE END OF           
C        EACH ITERATION IF IDEB IS 1.  SEE BELOW FOR ITS FORM.          
C                                                                       
C     IDEB = 0  MEANS DON'T EVER CALL MONIT.                            
C          = 1  MEANS CALL MONIT AT THE END OF EACH ITERATION.          
C                                                                       
C                                                                       
C     ON OUTPUT                                                         
C                                                                       
C        X  CONTAINS THE BEST AVAILABLE ESTIMATE OF THE SOLUTION VECTOR.
C                                                                       
C        FVECT  CONTAINS THE FUNCTION VALUES  F(J)  CORRESPONDING TO    
C           THE  X  VECTOR WITH THE EXCEPTION THAT THE CONTENTS ARE     
C           MEANINGLESS WHEN  IERR  .EQ. 10.                            
C                                                                       
C        RJACOB  CONTAINS THE ELEMENTS OF THE JACOBIAN MATRIX           
C           CORRESPONDING TO THE  X  VECTOR WITH THE EXCEPTION THAT THE 
C           CONTENTS ARE MEANINGLESS WHEN  IERR .EQ. 10.                
C                                                                       
C        RJINV  CONTAINS AN ESTIMATE OF THE NVAR X NVAR  COVARIANCE     
C            MATRIX.  THE ESTIMATE IS JUST                              
C               INVERSE( (JACOBIAN TRANSPOSE) * JACOBIAN )              
C                                                                       
C        IERR  IS A PARAMETER SET EQUAL TO  K  WHICH GIVES THE FOLLOWING
C           TERMINATION INDICATIONS                                     
C                                                                       
C           NORMAL TERMINATION,                                         
C             K = 0,                                                    
C           INTERMEDIATE TERMINATION,                                   
C             K = -N(N ANY INTEGER) - USER TERMINATION,                 
C             K =  1 - FAILURE TO CONVERGE IN  MAXF  CALLS OF  LMCFUN,  
C             K =  2 - DIRECTIONAL DERIVATIVE IS .GE. 0,                
C           INITIAL TERMINATION,                                        
C             K = 10 - NFUNS .LT. NVAR.                                 
C                                                                       
C        NFCALL IS  THE NUMBER OF CALLS TO  LMCFUN.                     
C                                                                       
C        NJCALL IS  THE NUMBER OF CALLS TO  EVLJAC.                     
C                                                                       
C        SCRV1  IS A SCRATCH VECTOR.                                    
C                                                                       
C        SQSCAL  CONTAINS THE SQUARE ROOTS OF THE ELEMENTS OF THE       
C           SCALE CONTROL VECTOR  S.                                    
C                                                                       
C        DJTJ  CONTAINS THE DIAGONAL ELEMENTS OF                        
C           (JACOBIAN TRANSPOSE) * JACOBIAN.                            
C                                                                       
C        VJTF  CONTAINS THE VECTOR  (JACOBIAN TRANSPOSE) * FVECT.       
C                                                                       
C        SCRV2  IS A SCRATCH VECTOR.                                    
C                                                                       
C     WRITTEN BY  K. E. HILLSTROM, SEPTEMBER, 1975.                     
C     10/31/75 - SLIGHT MODIFICATIONS - S.P.                            
C     1/1/76   - COMPUTE THE COVARIANCE MATRIX IN RJINV - S.P.          
C     2/6/76 - CHECK FOR STEP < TOL BEFORE LMCFUN - S. P.               
C     11/25/76 - ALLOW STEPS OUT OF BOUNDS - S.P.                       
C     2/7/77 - DO NOT STOP IMMEADIATELY AFTER AN OUT-OF-BOUNDS          
C                                                                       
C     ------------------------------------------------------------------
C                                                                       
C     TEST THE INPUT PARAMETERS  NFUNS  AND  NVAR  FOR CONSISTENCY      
C                                                                       
  100 IERR = 0                                                          
      IF (NFUNS .LT. NVAR) IERR = 10                                    
      IF (IERR .NE. 0) GO TO 325                                        
      IF (MODE .LT. 1 .OR. MODE .GT. 3) MODE = 1                        
C                                                                       
C     INITIALIZE THE FOLLOWING PARAMETERS                               
C                                                                       
C     QPARM  - THE MARQUARDT ITERATION CONTROL PARAMETER                
C     QCPARM - THE MARQUARDT PARAMETER CUT-OFF VALUE                    
C     NFCALL - THE NUMBER OF CALLS TO  LMCFUN                           
C     NJCALL - THE NUMBER OF CALLS TO  EVLJAC                           
C                                                                       
      QPARM = 0.0D0                                                     
      QCPARM = 1.0D0                                                    
      NFCALL = 1                                                        
      NJCALL = 1                                                        
      ITCNT = -1                                                        
C                                                                       
C     PREPARE FOR THE FIRST ITERATION BY INITIALIZING THE FUNCTION      
C       VALUES , THE SUM OF SQUARES  FSQ  AND THE JACOBIAN MATRIX       
C                                                                       
C     ******************************************************************
      ITYPE = 3                                                         
      CALL LMCFUN ( ITYPE, NVAR,NFUNS,X,FVECT, RJACOB, IERR )           
C     ******************************************************************
C                                                                       
      IF (IERR .LT. 0) GO TO 325                                        
C                                                                       
CXXXX C     ************************************************************
CXXXX       CALL EVLJAC(NVAR,NFUNS,X,RJACOB,IERR)                       
CXXXX C     ************************************************************
CXXXX       IF (IERR .LT. 0) GO TO 325                                  
C                                                                       
C     ******************************************************************
      CALL SIP(FVECT,1,2,FVECT,1,2,0.0D0,FSQ,NFUNS)                     
C     ******************************************************************
C     CALL  CMTRIX  TO FORM THE MATRIX                                  
C       JTJ  = (JACOBIAN TRANSPOSE) * JACOBIAN                          
C       AND THE VECTOR                                                  
C       JTF  = (JACOBIAN TRANSPOSE) * FVECT.                            
C                                                                       
C       IN THIS SECTION                                                 
C                                                                       
C       RJINV  - CONTAINS THE LOWER TRIANGLE OF THE SYMMETRIC MATRIX    
C                JTJ                                                    
C                                                                       
C       ALSO, FROM THIS SECTION ON                                      
C                                                                       
C       VJTF(I)  - CONTAINS THE VECTOR  JTF                             
C                                                                       
C     ******************************************************************
      CALL CMTRIX(NVAR,NFUNS,RJACOB,FVECT,RJINV,VJTF)                   
C     ******************************************************************
C                                                                       
C     INITIALIZE THE SCALE CONTROL VECTOR  S  ACCORDING TO THE VALUE    
C       OF  MODE. FROM THIS SECTION ON                                  
C                                                                       
C     SQSCAL(I) - CONTAINS THE SQUARE ROOTS OF THE ELEMENTS OF  S       
C                                                                       
      IF (MODE .EQ. 3) GO TO 110                                        
      IF (MODE .EQ. 2) GO TO 120                                        
C                                                                       
C     MODE = 1, SET  S(I) = JTJ(I,I) OR 1                               
C                                                                       
      DO 105 I = 1, NVAR                                                
      S(I) = RJINV(I,I)                                                 
      IF (S(I) .LE. 0.0D0) S(I) = 1.0D0                                 
  105 CONTINUE                                                          
C                                                                       
C     MODE = 3, S(I) HAS BEEN SET BY THE USER THROUGH THE PARAMETER  S  
C                                                                       
  110 DO 115 I = 1, NVAR                                                
      SQSCAL(I) = DSQRT(S(I))                                           
  115 CONTINUE                                                          
C                                                                       
      GO TO 200                                                         
C                                                                       
C     MODE = 2 SET S(I) = 1                                             
C                                                                       
  120 DO 125 I = 1, NVAR                                                
      SQSCAL(I) = 1.0D0                                                 
      S(I) = 1.0D0                                                      
  125 CONTINUE                                                          
C                                                                       
C     BEGIN AN ITERATION BY PLACING DIAGONAL ELEMENTS OF  JTJ  IN       
C       DJTJ(I)                                                         
C                                                                       
C     ALLOW USER TO PRINT STUFF EACH ITERATION                          
C                                                                       
 200  ITCNT = ITCNT + 1                                                 
      IF ( IDEB .EQ. 1 ) CALL MONIT ( ITCNT, NFCALL, NJCALL, NVAR,  NFUN
     *S, X, FSQ, FVECT )                                                
      NFLAST = NFCALL                                                   
C                                                                       
      DO 205 I = 1, NVAR                                                
      DJTJ(I) = RJINV(I,I)                                              
  205 CONTINUE                                                          
C                                                                       
C     CALL  QMTRIX  TO SOLVE THE MARQUARDT MATRIX EQUATION              
C                                                                       
C       (JTJ + QPARM * S) * DELTA = - JTF.                              
C                                                                       
C       NOTE THAT  QMTRIX  INVOKES THE  ANL  LIBRARY ROUTINE            
C       ANL F152S, SYMINV  TO SOLVE THE MATRIX EQUATION.                
C                                                                       
C       WHEN FINISHED                                                   
C                                                                       
C       RJINV    - CONTAINS THE UPPER TRIANGLE OF                       
C                  ((JTJ + QPARM * S) INVERSE)  AND THE STRICT LOWER    
C                  TRIANGLE OF  JTJ                                     
C       SCRV2(I) - CONTAINS THE NEGATIVE OF THE CORRECTION VECTOR       
C                  DELTA                                                
C                                                                       
C                                                                       
C     ******************************************************************
  210 CALL QMTRIX(NVAR,RJINV,DJTJ,VJTF,S,SCRV2,QPARM)                   
C     ******************************************************************
C                                                                       
C     ERROR RETURN IF THE DIRECTIONAL DERIVATIVE ALONG THE CORRECTION,  
C       (JTF TRANSPOSE) * DELTA, IS .GE. 0.                             
C                                                                       
C     ******************************************************************
      CALL SIP(VJTF,1,2,SCRV2,1,2,0.0D0,DIRDRV,NVAR)                    
C     ******************************************************************
C                                                                       
      IF (DIRDRV .LE. 0.0D0) GO TO 290                                  
C                                                                       
C     COMPUTE  PRERED  THE PREDICTED REDUCTION IN THE SUM OF SQUARES    
C                                                                       
      DO 230 I = 1, NVAR                                                
      Z = DJTJ(I) * SCRV2(I)                                            
C                                                                       
C     ******************************************************************
      IF (I .GT. 1) CALL SIP(RJINV(I,1),1,NVAR+1,SCRV2,1,2,Z,Z,  I-1)   
C     ******************************************************************
      IF (I.LT.NVAR) CALL SIP(RJINV(I,I),2,3,SCRV2(I),2,3,Z,Z,NVAR-I)                                                                 
C     ******************************************************************
C                                                                       
      SCRV1(I) = 2.0D0 * VJTF(I) - Z                                    
  230 CONTINUE                                                          
C                                                                       
C     ******************************************************************
      CALL SIP(SCRV1,1,2,SCRV2,1,2,0.0D0,PRERED,NVAR)                   
C     ******************************************************************
C                                                                       
C     ICONV  = 0 IF  DABS(DELTA) .LE. MAX(TOL*DABS(X(I)),TOL)  FOR ALL  
C       I. FROM THIS SECTION ON                                         
C                                                                       
C       SCRV1(I)  -  CONTAINS THE ESTIMATE (XK + DELTA)                 
C                                                                       
      ICONV = 0                                                         
C                                                                       
      DO 235 I = 1, NVAR                                                
      Z = TOL * DABS(X(I))                                              
      IF (Z .LE. TOL2) Z = TOL2                                         
      IF (DABS(SCRV2(I)) .GT. Z) ICONV = 1                              
      SCRV1(I) = X(I) - SCRV2(I)                                        
  235 CONTINUE                                                          
C                                                                       
C     NORMAL EXIT IF  ICONV = 0  OR THE PREDICTED REDUCTION IS .LE. 0.  
C                                                                       
C     DO NOT STOP IF THE LAST STEP WAS OUT OF BOUNDS                    
C                                                                       
      IF ( IERR .EQ. 1 ) GO TO 240                                      
      IF ((ICONV .EQ. 0) .OR. (PRERED .LE. 0.0D0)) GO TO 310            
C                                                                       
C     TEST WHETHER THERE HAVE BEEN  MAXFCN  CALLS OF  LMCFUN            
C                                                                       
 240  IF (NFCALL .GE. MAXF) GO TO 295                                   
C                                                                       
C     COMPUTE  FSQP  THE SUM OF SQUARES AT THE TRIAL ESTIMATE. FROM     
C       THIS SECTION ON                                                 
C                                                                       
C       SCRV2(I)  -  CONTAINS THE FUNCTIONS AT THE TRIAL ESTIMATE       
C                                                                       
C     ******************************************************************
      ITYPE = 1                                                         
      CALL LMCFUN ( ITYPE, NVAR,NFUNS,SCRV1,SCRV2, RJACOB, IERR )       
C     ******************************************************************
C                                                                       
      NFCALL = NFCALL + 1                                               
      IF (IERR .LT. 0) GO TO 310                                        
      IF ( IERR .NE. 1 ) GO TO 260                                      
C                                                                       
C     STEP WAS OUT OF BOUNDS                                            
C                                                                       
      QPARM = 4.0D0*QPARM                                                   
      IF ( QPARM .LT. 0.01D0 ) QPARM = 0.01D0                                 
      QCPARM = 0.01D0*QPARM                                                
      GO TO 210                                                         
C                                                                       
C     ******************************************************************
 260  CALL SIP(SCRV2,1,2,SCRV2,1,2,0.0D0,FSQP,NFUNS)                    
C     ******************************************************************
C                                                                       
C     COMPUTE  ACTRED  THE ACTUAL REDUCTION IN THE SUM OF SQUARES AND   
C       TEST WHETHER TRIAL SUM OF SQUARES IS LEAST                      
C                                                                       
      ACTRED = FSQ - FSQP                                               
C                                                                       
C     CALL  UPQPAR  TO UPDATE THE MARQUARDT PARAMETER                   
C                                                                       
C     ******************************************************************
      CALL UPQPAR(NVAR,SQSCAL,RJINV,ACTRED,PRERED,DIRDRV,QPARM,QCPARM)                                                                  
C     ******************************************************************
C                                                                       
C     TEST WHETHER TRIAL SUM OF SQUARES IS LEAST                        
C                                                                       
      IF (ACTRED .LE. 0.0D0) GO TO 210                                  
C                                                                       
C     THE TRIAL SUM OF SQUARES IS LEAST. SET FSQ,  X  AND FVECT  TO     
C       BEST CALCULATED VALUES.                                         
C                                                                       
      DO 275 I = 1, NVAR                                                
      X(I) = SCRV1(I)                                                   
  275 CONTINUE                                                          
C                                                                       
      DO 280 I = 1, NFUNS                                               
      FVECT(I) = SCRV2(I)                                               
  280 CONTINUE                                                          
C                                                                       
C     THE USER MAY HAVE ELECTED TO COMPUTE THE JACOBIAN ON THE PREVIOUS 
C     FUNCTION REQUEST.                                                 
C                                                                       
      IF ( ITYPE .EQ. 3 ) GO TO 285                                     
C                                                                       
C     ******************************************************************
      ITYPE = 2                                                         
      CALL LMCFUN ( ITYPE, NVAR,NFUNS,X, FVECT, RJACOB, IERR )          
C     ******************************************************************
C                                                                       
 285  NJCALL = NJCALL + 1                                               
      IF (IERR .LT. 0) GO TO 310                                        
C                                                                       
C     SHOULD WE STOP DUE TO INSIGNIFICANT CHANGES IN FSQ.               
C                                                                       
      IF ( ITCNT .LE. 2 ) GO TO 270                                     
      IF ( FSQP .LT. TOL2 ) GO TO 300                                   
      IF ( FSQM1-FSQP .LT. TOL*FSQP ) GO TO 300                         
C                                                                       
C     KEEP GOING                                                        
C                                                                       
 270  FSQM1 = FSQ                                                       
      FSQ = FSQP                                                        
C                                                                       
C     CALL  CMTRIX  TO FORM NEW  JTJ  AND  JTF  AND START A NEW         
C       ITERATION                                                       
C                                                                       
C     ******************************************************************
      CALL CMTRIX(NVAR,NFUNS,RJACOB,FVECT,RJINV,VJTF)                   
C     ******************************************************************
C                                                                       
      GO TO 200                                                         
C                                                                       
C     ERROR RETURN BECAUSE DIRECTIONAL DERIVATIVE IS .GE. 0             
C                                                                       
  290 IERR = 2                                                          
      GO TO 325                                                         
C                                                                       
C     ERROR RETURN BECAUSE THERE HAVE BEEN  MAXF  CALLS OF  LMCFUN      
C                                                                       
  295 IERR = 1                                                          
      GO TO 310                                                         
C                                                                       
C     CONVERGENCE CRITERION HAS BEEN SATISFIED OR THE PREDICTED         
C       REDUCTION IN THE SUM OF SQUARES IS .LE. 0 OR THE FUNCTION CALL  
C       LIMIT HAS BEEN REACHED. X  AND  FVECT ARE UPDATED ONLY IF THE   
C       FINAL REDUCTION IN  FSQ  IS POSITIVE                            
C                                                                       
  300 ITCNT = ITCNT + 1                                                 
C                                                                       
C     ALLOW USER TO PRINT STUFF EACH ITERATION                          
C                                                                       
 310  IF ( IDEB .EQ. 1 .AND. NFCALL .GT. NFLAST )  CALL MONIT ( ITCNT,
     * NFCALL, NJCALL, NVAR,  NFUNS, X, FSQP, FVECT )                     
C                                                                       
C     COMPUTE  INVERSE(JTJ)  AS AN ESTIMATE OF THE COVARIANCE MATRIX    
C     AT THE FINAL VALUE.                                               
C                                                                       
 325  CALL QMTRIX ( NVAR, RJINV, DJTJ, VJTF, S, SCRV2, 0.D0 )           
C                                                                       
      DO 339 I = 1, NVAR                                                
      DO 339 J = 1, I                                                   
      RJINV(I,J) = RJINV(J,I)                                           
 339  CONTINUE                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE LMCFUN ( ITYPE, MPARA, MPTS, X, FS, FJACOB, IERR )     
C                                                                       
C     CHI. SQR. FUNCTION FOR  LMCHOL                                    
C                                                                       
C     THIS IS AN INTERFACE BETWEEN LMCHOL AND CALFUN.                   
C                                                                       
C     10/31/75 - FIRST VERSION - S. P.                                  
C     11/25/76 - RETURN IERR = 1 FOR STEPS OUT OF BOUNDS                
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,FS,FJACOB                                         
C                                                                       
CK      REAL*8 X(MPARA), FS(MPTS), FJACOB(MPTS, MPARA)                    
      DIMENSION X(MPARA), FS(MPTS), FJACOB(MPTS, MPARA)                    
C                                                                       
      COMMON/LOCFIT/IDATA,LPL,LSIGMA,LSMAT,IESETS,IMASK,IEQUAL,NPARA,
     *IPARA,ITEMP,LSCSIG, LEPARM, LFS, NPTS, LEQNAM, MAXEQL,MAXVS, 
     *LVDERV, LSDELS, LDVTYP, LDSDP, LDFDX,  LFGRAD, LTGRAD
C                                                                       
      COMMON /ALLOCS/ ALLOC(1)                                          
CDC      ALLOCLEVEL  ALLOC                                              
C                                                                       
      IERR = 0                                                          
C
C      GO TO ( 10, 20, 20 ), ITYPE                                   
C                                                                   
C 10   CALL CALFUN ( MPTS, MPARA, FS, X )                            
C      RETURN                                                        
C                                                                       
C     IT APPEARS TO BE MOST ECONOMICAL TO ALWAYS COMPUTE EVERYTHING     
C     SINCE LMCHOL SUCCEEDS MOST OF THE TIME AND THEN WANTS THE JACOBIAN
C     ANYWAY.                                                           
C                                                                       
 20   F = CALGRA ( MPTS, MPARA, FS, X, FJACOB, ALLOC(LTGRAD) )          
      ITYPE = 3                                                         
      IF ( F .GT. 1.D20 ) IERR = 1                                      
C                                                                       
      RETURN                                                            
      END                                                               
C     ------------------------------------------------------------------
C                                                                       
      SUBROUTINE MINFIT(NM,M,N,A,W,IP,B,IERR,RV1)                       
C                                                                       
CDC      ALLOCLEVEL A,W,B,RV1                                           
      INTEGER I,J,K,L,M,N,II,IP,I1,KK,K1,LL,L1,M1,NM,ITS,IERR           
      REAL*8 A(NM,N),W(N),B(NM,IP),RV1(N)                               
      REAL*8 C,F,G,H,S,X,Y,Z,EPS,SCALE,MACHEP                           
C      REAL*8 DSQRT,DMAX1,DABS,DSIGN                                     
C                                                                       
C     THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE MINFIT,   
C     NUM. MATH. 14, 403-420(1970) BY GOLUB AND REINSCH.                
C     HANDBOOK FOR AUTO. COMP., VOL II-LINEAR ALGEBRA, 134-151(1971).   
C                                                                       
C     THIS SUBROUTINE DETERMINES, TOWARDS THE SOLUTION OF THE LINEAR    
C                                                        T              
C     SYSTEM AX=B, THE SINGULAR VALUE DECOMPOSITION A=USV  OF A REAL    
C                                         T                             
C     M BY N RECTANGULAR MATRIX, FORMING U B RATHER THAN U.  HOUSEHOLDER
C     BIDIAGONALIZATION AND A VARIANT OF THE QR ALGORITHM ARE USED.     
C                                                                       
C     ON INPUT:                                                         
C                                                                       
C        NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL         
C          ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM          
C          DIMENSION STATEMENT.  NOTE THAT NM MUST BE AT LEAST          
C          AS LARGE AS THE MAXIMUM OF M AND N;                          
C                                                                       
C        M IS THE NUMBER OF ROWS OF A AND B;                            
C                                                                       
C        N IS THE NUMBER OF COLUMNS OF A AND THE ORDER OF V;            
C                                                                       
C        A CONTAINS THE RECTANGULAR COEFFICIENT MATRIX OF THE SYSTEM;   
C                                                                       
C        IP IS THE NUMBER OF COLUMNS OF B.  IP CAN BE ZERO;             
C                                                                       
C        B CONTAINS THE CONSTANT COLUMN MATRIX OF THE SYSTEM            
C          IF IP IS NOT ZERO.  OTHERWISE B IS NOT REFERENCED.           
C                                                                       
C     ON OUTPUT:                                                        
C                                                                       
C        A HAS BEEN OVERWRITTEN BY THE MATRIX V (ORTHOGONAL) OF THE     
C          DECOMPOSITION IN ITS FIRST N ROWS AND COLUMNS.  IF AN        
C          ERROR EXIT IS MADE, THE COLUMNS OF V CORRESPONDING TO        
C          INDICES OF CORRECT SINGULAR VALUES SHOULD BE CORRECT;        
C                                                                       
C        W CONTAINS THE N (NON-NEGATIVE) SINGULAR VALUES OF A (THE      
C          DIAGONAL ELEMENTS OF S).  THEY ARE UNORDERED.  IF AN         
C          ERROR EXIT IS MADE, THE SINGULAR VALUES SHOULD BE CORRECT    
C          FOR INDICES IERR+1,IERR+2,...,N;                             
C                                                                       
C                                   T                                   
C        B HAS BEEN OVERWRITTEN BY U B.  IF AN ERROR EXIT IS MADE,      
C                       T                                               
C          THE ROWS OF U B CORRESPONDING TO INDICES OF CORRECT          
C          SINGULAR VALUES SHOULD BE CORRECT;                           
C                                                                       
C        IERR IS SET TO                                                 
C          ZERO       FOR NORMAL RETURN,                                
C          K          IF THE K-TH SINGULAR VALUE HAS NOT BEEN           
C                     DETERMINED AFTER 30 ITERATIONS;                   
C                                                                       
C        RV1 IS A TEMPORARY STORAGE ARRAY.                              
C                                                                       
C     QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW,        
C     APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY         
C                                                                       
C     12/28/79 - STANDARD MORTRAN - RPG                                 
C                                                                       
C     ------------------------------------------------------------------
C                                                                       
C     ;;;;;; MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING         
C                THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC.   
C                MACHEP = 16.0D0**(-13) FOR LONG FORM ARITHMETIC        
C                ON S360                                                
CIB      DATA MACHEP/Z3410000000000000/                                 
C                                                                       
      DATA MACHEP / 1.D-14 /                                            
C                                                                       
      IERR = 0                                                          
C                                                                       
C           ;;;;;; HOUSEHOLDER REDUCTION TO BIDIAGONAL FORM             
C                                                                       
      G = 0.0D0                                                         
      SCALE = 0.0D0                                                     
      X = 0.0D0                                                         
C                                                                       
      DO 300 I = 1, N                                                   
      L = I + 1                                                         
      RV1(I) = SCALE * G                                                
      G = 0.0D0                                                         
      S = 0.0D0                                                         
      SCALE = 0.0D0                                                     
      IF (I .GT. M) GO TO 210                                           
C                                                                       
      DO 120 K = I, M                                                   
  120 SCALE = SCALE + DABS(A(K,I))                                      
C                                                                       
      IF (SCALE .EQ. 0.0D0) GO TO 210                                   
C                                                                       
      DO 130 K = I, M                                                   
      A(K,I) = A(K,I) / SCALE                                           
      S = S + A(K,I)**2                                                 
  130 CONTINUE                                                          
C                                                                       
      F = A(I,I)                                                        
      G = -DSIGN(DSQRT(S),F)                                            
      H = F * G - S                                                     
      A(I,I) = F - G                                                    
      IF (I .EQ. N) GO TO 160                                           
C                                                                       
      DO 150 J = L, N                                                   
      S = 0.0D0                                                         
C                                                                       
      DO 140 K = I, M                                                   
  140 S = S + A(K,I) * A(K,J)                                           
C                                                                       
      F = S / H                                                         
C                                                                       
      DO 150 K = I, M                                                   
      A(K,J) = A(K,J) + F * A(K,I)                                      
  150 CONTINUE                                                          
C                                                                       
  160 IF (IP .EQ. 0) GO TO 190                                          
C                                                                       
      DO 180 J = 1, IP                                                  
      S = 0.0D0                                                         
C                                                                       
      DO 170 K = I, M                                                   
  170 S = S + A(K,I) * B(K,J)                                           
C                                                                       
      F = S / H                                                         
C                                                                       
      DO 180 K = I, M                                                   
      B(K,J) = B(K,J) + F * A(K,I)                                      
  180 CONTINUE                                                          
C                                                                       
  190 DO 200 K = I, M                                                   
  200 A(K,I) = SCALE * A(K,I)                                           
C                                                                       
  210 W(I) = SCALE * G                                                  
      G = 0.0D0                                                         
      S = 0.0D0                                                         
      SCALE = 0.0D0                                                     
      IF (I .GT. M .OR. I .EQ. N) GO TO 290                             
C                                                                       
      DO 220 K = L, N                                                   
  220 SCALE = SCALE + DABS(A(I,K))                                      
C                                                                       
      IF (SCALE .EQ. 0.0D0) GO TO 290                                   
C                                                                       
      DO 230 K = L, N                                                   
      A(I,K) = A(I,K) / SCALE                                           
      S = S + A(I,K)**2                                                 
  230 CONTINUE                                                          
C                                                                       
      F = A(I,L)                                                        
      G = -DSIGN(DSQRT(S),F)                                            
      H = F * G - S                                                     
      A(I,L) = F - G                                                    
C                                                                       
      DO 240 K = L, N                                                   
  240 RV1(K) = A(I,K) / H                                               
C                                                                       
      IF (I .EQ. M) GO TO 270                                           
C                                                                       
      DO 260 J = L, M                                                   
      S = 0.0D0                                                         
C                                                                       
      DO 250 K = L, N                                                   
  250 S = S + A(J,K) * A(I,K)                                           
C                                                                       
      DO 260 K = L, N                                                   
      A(J,K) = A(J,K) + S * RV1(K)                                      
  260 CONTINUE                                                          
C                                                                       
  270 DO 280 K = L, N                                                   
  280 A(I,K) = SCALE * A(I,K)                                           
C                                                                       
  290 X = DMAX1(X,DABS(W(I))+DABS(RV1(I)))                              
  300 CONTINUE                                                          
C                                                                       
C           ;;;;;; ACCUMULATION OF RIGHT-HAND TRANSFORMATIONS.          
C                FOR I=N STEP -1 UNTIL 1 DO --                          
C                                                                       
      DO 400 II = 1, N                                                  
      I = N + 1 - II                                                    
      IF (I .EQ. N) GO TO 390                                           
      IF (G .EQ. 0.0D0) GO TO 360                                       
C                                                                       
      DO 320 J = L, N                                                   
C                                                                       
C                                                                       
  320 A(J,I) = (A(I,J) / A(I,L)) / G                                    
C                                                                       
      DO 350 J = L, N                                                   
      S = 0.0D0                                                         
C                                                                       
      DO 340 K = L, N                                                   
  340 S = S + A(I,K) * A(K,J)                                           
C                                                                       
      DO 350 K = L, N                                                   
      A(K,J) = A(K,J) + S * A(K,I)                                      
  350 CONTINUE                                                          
C                                                                       
  360 DO 380 J = L, N                                                   
      A(I,J) = 0.0D0                                                    
      A(J,I) = 0.0D0                                                    
  380 CONTINUE                                                          
C                                                                       
  390 A(I,I) = 1.0D0                                                    
      G = RV1(I)                                                        
      L = I                                                             
  400 CONTINUE                                                          
C                                                                       
      IF (M .GE. N .OR. IP .EQ. 0) GO TO 510                            
      M1 = M + 1                                                        
C                                                                       
      DO 500 I = M1, N                                                  
C                                                                       
      DO 500 J = 1, IP                                                  
      B(I,J) = 0.0D0                                                    
  500 CONTINUE                                                          
C                                                                       
C        ;;;;;; DIAGONALIZATION OF THE BIDIAGONAL FORM                  
C                                                                       
  510 EPS = MACHEP * X                                                  
C                                                                       
C                                                                       
      DO 700 KK = 1, N                                                  
      K1 = N - KK                                                       
      K = K1 + 1                                                        
      ITS = 0                                                           
C                                                                       
C       ;;;; TEST FOR SPLITTING.                                    0023
C                FOR L=K STEP -1 UNTIL 1 DO --                          
  520 DO 530 LL = 1, K                                                  
      L1 = K - LL                                                       
      L = L1 + 1                                                        
      IF (DABS(RV1(L)) .LE. EPS) GO TO 565                              
C                                                                       
C       ;;;;; RV1(1) IS ALWAYS ZERO, SO THERE IS NO EXIT             002
C                THROUGH THE BOTTOM OF THE LOOP                         
      IF (DABS(W(L1)) .LE. EPS) GO TO 540                               
  530 CONTINUE                                                          
C                                                                       
C      ;;;;; CANCELLATION OF RV1(L) IF L GREATER THAN 1                 
  540 C = 0.0D0                                                         
      S = 1.0D0                                                         
C                                                                       
      DO 560 I = L, K                                                   
      F = S * RV1(I)                                                    
      RV1(I) = C * RV1(I)                                               
      IF (DABS(F) .LE. EPS) GO TO 565                                   
      G = W(I)                                                          
      H = DSQRT(F*F+G*G)                                                
      W(I) = H                                                          
      C = G / H                                                         
      S = -F / H                                                        
      IF (IP .EQ. 0) GO TO 560                                          
C                                                                       
      DO 550 J = 1, IP                                                  
      Y = B(L1,J)                                                       
      Z = B(I,J)                                                        
      B(L1,J) = Y * C + Z * S                                           
      B(I,J) = -Y * S + Z * C                                           
  550 CONTINUE                                                          
C                                                                       
  560 CONTINUE                                                          
C                                                                       
C      ;;;;; TEST FOR CONVERGENCE                                       
  565 Z = W(K)                                                          
      IF (L .EQ. K) GO TO 650                                           
C                                                                       
C    ;;;;; SHIFT FROM BOTTOM 2 BY 2 MINOR                               
      IF (ITS .EQ. 30) GO TO 1000                                       
      ITS = ITS + 1                                                     
      X = W(L)                                                          
      Y = W(K1)                                                         
      G = RV1(K1)                                                       
      H = RV1(K)                                                        
      F = ((Y - Z) * (Y + Z) + (G - H) * (G + H)) / (2.0D0 * H * Y)     
      G = DSQRT(F*F+1.0D0)                                              
      F = ((X - Z) * (X + Z) + H * (Y / (F + DSIGN(G,F)) - H)) / X      
C                                                                       
C          ;;;;; NEXT QR TRANSFORMATION                                 
      C = 1.0D0                                                         
      S = 1.0D0                                                         
C                                                                       
      DO 600 I1 = L, K1                                                 
      I = I1 + 1                                                        
      G = RV1(I)                                                        
      Y = W(I)                                                          
      H = S * G                                                         
      G = C * G                                                         
      Z = DSQRT(F*F+H*H)                                                
      RV1(I1) = Z                                                       
      C = F / Z                                                         
      S = H / Z                                                         
      F = X * C + G * S                                                 
      G = -X * S + G * C                                                
      H = Y * S                                                         
      Y = Y * C                                                         
C                                                                       
      DO 570 J = 1, N                                                   
      X = A(J,I1)                                                       
      Z = A(J,I)                                                        
      A(J,I1) = X * C + Z * S                                           
      A(J,I) = -X * S + Z * C                                           
  570 CONTINUE                                                          
C                                                                       
      Z = DSQRT(F*F+H*H)                                                
      W(I1) = Z                                                         
C                                                                       
C         ;;;;; ROTATION CAN BE ARBITRARY IF Z IS ZERO                  
      IF (Z .EQ. 0.0D0) GO TO 580                                       
      C = F / Z                                                         
      S = H / Z                                                         
  580 F = C * G + S * Y                                                 
      X = -S * G + C * Y                                                
      IF (IP .EQ. 0) GO TO 600                                          
C                                                                       
      DO 590 J = 1, IP                                                  
      Y = B(I1,J)                                                       
      Z = B(I,J)                                                        
      B(I1,J) = Y * C + Z * S                                           
      B(I,J) = -Y * S + Z * C                                           
  590 CONTINUE                                                          
C                                                                       
  600 CONTINUE                                                          
C                                                                       
      RV1(L) = 0.0D0                                                    
      RV1(K) = F                                                        
      W(K) = X                                                          
      GO TO 520                                                         
C                                                                       
C      ;;;;; CONVERGENCE ;;;;;;;;;:                                 0023
  650 IF (Z .GE. 0.0D0) GO TO 700                                       
C                                                                       
C       ;;;;; W(K) IS MADE NON-NEGATIVE                                 
      W(K) = -Z                                                         
C                                                                       
      DO 690 J = 1, N                                                   
  690 A(J,K) = -A(J,K)                                                  
C                                                                       
  700 CONTINUE                                                          
C                                                                       
      GO TO 1001                                                        
C                                                                       
C       ;;;;; SET ERROR -- NO CONVERGENCE TO A                       002
C                SINGULAR VALUE AFTER 30 ITERATIONS                     
 1000 IERR = K                                                          
 1001 RETURN                                                            
C                                                                       
C         ;;;;; LAST CARD OF MINFIT                                     
      END                                                               

      SUBROUTINE MINIM(DIAG,EPS,X,XP,G,GP,S,SIG,Y,ETA,RHO,H,FS,FJACOB,
     * NA,N,IPRINT,MAXFUN,ISIG)                                           
C     VARIABLE-METRIC MINIMIZATION ; DX = CONST * (-H * G)              
C     WHERE G IS LOCAL DERIVATIVE ,H IS INVERSE SECOND-DERIVATIVE MATRIX
C     H UPDATED BY DAVIDON OR FLETCHER PRESCRIPTION                     
C                                                                       
C     MHM DECEMBER 29 1975                                              
C                                                                       
C     DIAG : IF GE.0 MINIMUM NOT ALWAYS LOCATED ALONG -H*G              
C            IF LT.0 MINIMUM LOCATED ALONG -H*G AT EACH STEP            
C            /DIAG/ CONTROLS METRIC MODIFICATION AND STEP VATIATION     
C            IF /DIAG/.GT.1 FLETCHER WITH FACTOR=1/(/DIAG/)             
C            IF /DIAG/.LT.1 DAVIDON WITH FACTOR=/DIAG/                  
C            IF /DIAG/.EQ.0  OR 1 FACTOR=1/3                            
C     EPS  : RELATIVE ACCURACY REQUESTED                                
C     X    : VECTOR OF VARIABLE-VALUES  ( DIMENSION N)                  
C            INITIALLY CONTAINS GUESS AT POSITION OF MINIMUM            
C            RETURNS VALUES AT LAST ITERATION ; THESE WILL BE EITHER    
C            VALUES AT MIN TO ACCURACY EPS, OR A BEST APPROACH TO THESE 
C            IN MAXFUN CALLS OF SUBROUTINE CALGRA                       
C     XP,G,GP,S,SIG,Y,ETA,RHO : 8 WORK ARRAYS OF LENGTH N               
C     H    : N BY N MATRIX - RETURNS VARIANCE-COVARIANCE MATRIX         
C     FS   : VECTOR (DIMENSION NA) OF FUNCTIONS IN SUM-OF-SQUARES       
C     FJACOB: WORK ARRAY OF DIMENSION NA BY N                           
C     N    : NO OF PARAMETERS                                           
C     NA   : NO OF TERMS IN CHISQ SUM                                   
C     IPRINT: PRINT CONTROL (USES MINMON)                               
C     MAXFUN : UPPER LIMIT ON NUMBER OF FUNCTION CALLS                  
C     ISIG : AT RETURN ISIG = 0 MEANS NORMAL COMPLETION                 
C                      ISIG = 1 - STOPPED AFTER MAXFUN FUNCTION CALLS   
C                      ISIG =-5 - ABNORMAL CONDITIONS ENCOUNTERED       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,SIG,Y,ETA,RHO,H,FS,FJACOB               
C                                                                       
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),SIG(N),Y(N),ETA(N),RHO(N), 
     * H(N,N),FS(NA),FJACOB(NA,N)                                          
C                                                                       
C     INITIAL FUNCTION-CALL                                             
C                                                                       
      IF(DIAG.GE.0) ND=0                                                
      IF(DIAG.LT.0) ND=1                                                
      DG=DABS(DIAG)                                                     
      IF(DG.LT.1.0D0) W=DG                                                  
      IF((DG.EQ.0.0D0).OR.(DG.EQ.1.0D0)) W=1.0D0/3.0D0                                
      IF(DG.GT.1.0D0) W=1.0D0/DG                                                
      ITYP=0                                                            
      IF(DG.LT.1.0D0) ITYP=1                                                
      WRITE(6,3985)
 3985 FORMAT(1H0,2X,'++ VARIABLE-METRIC MINIMIZATION ++')               
      IF((ND.EQ.0).AND.(ITYP.EQ.0)) WRITE(6,4050)W                        
 4050 FORMAT(1H ,' FLETCHER METRIC,INCOMPLETE LINEAR SEARCH, W=',G12.6) 
      IF((ND.EQ.0).AND.(ITYP.EQ.1)) WRITE(6,4450)W                        
 4450 FORMAT(1H ,' DAVIDON METRIC,INCOMPLETE LINEAR SEARCH, W=',G12.6)  
      IF((ND.EQ.1).AND.(ITYP.EQ.1)) WRITE(6,4650)W                        
 4650 FORMAT(1H ,' DAVIDON METRIC,COMPLETE LINEAR SEARCH, W=',G12.6)    
      IF((ND.EQ.1).AND.(ITYP.EQ.0)) WRITE(6,4850)W                        
 4850 FORMAT(1H ,' FLETCHER METRIC,COMPLETE LINEAR SEARCH, W=',G12.6)   
      ISCREW=1                                                          
      ISIG=10                                                           
      NIT=1                                                             
      AA=1.0D0                                                              
      UPRLIM=1.D+18                                                     
      NCALL=1                                                           
      DO 5 I=1,N                                                        
    5 XP(I)=X(I)                                                        
      CALL CALGRA(NA,N,FS,XP,FJACOB,GP)                                 
      FXP=0.0D0                                                             
      DO 10 I=1,NA                                                      
   10 FXP=FXP+FS(I)*FS(I)                                               
      FX1=FXP                                                           
      FX2=0.0D0                                                             
      FX3=0.0D0                                                             
C                                                                       
C     INITIAL ESTIMATE OF METRIC MATRIX                                 
C                                                                       
   30 DO 50 I=1,N                                                       
      DO 50 J=I,N                                                       
      H(I,J)=0.0D0                                                          
      DO 40 IA=1,NA                                                     
   40 H(I,J)=H(I,J)+FJACOB(IA,I)*FJACOB(IA,J)                           
      H(I,J)=2.0D0*H(I,J)                                                   
   50 H(J,I)=H(I,J)                                                     
      CALL SYMINV(H,N,Y,0,DET,N)                                        
      DO 60 I=2,N                                                       
      JTOP=I-1                                                          
      DO 60 J=1,JTOP                                                    
   60 H(I,J)=H(J,I)                                                     
      GO TO 90                                                          
C                                                                       
C     PRINT-OUT OF INFORMATION AT START OF ITERATION                    
C                                                                       
   70 ISIG=0                                                            
      IF(IPRINT.EQ.1) CALL MINMON(NIT,NCALL,NCALL,N,NA,X,FX,FJACOB)     
C                                                                       
C     MAIN ITERATIVE SECTION  : INCOMPLETE LINEAR SEARCH                
C                                                                       
  132 ALPHA=AA                                                          
      IF((AA.LT.1.0D0).AND.(ND.EQ.0)) AA=W                                  
      NTURN=1                                                           
  120 SIGX=ALPHA*SXNEW                                                  
      IF(SIGX.GE.EPS/10.0D0) GO TO 1120                                     
      ISIG=-5                                                           
      WRITE(6,7000)SIGX,NTURN,NIT                                         
 7000 FORMAT(1H0,2X'-- MINIM TRIED TO TAKE STEP OF LENGTH ',G12.6/ 3X,'A
     *T THE ',I2,'TH TURN OF ITERATION ',I3,'--')                       
      GO TO 250                                                         
 1120 DO 130 I=1,N                                                      
      SIG(I)=ALPHA*S(I)                                                 
  130 XP(I)=X(I)+SIG(I)                                                 
      IF ( NCALL .GE. MAXFUN ) GO TO 215                                
      CALL CALGRA(NA,N,FS,XP,FJACOB,GP)                                 
      NCALL=NCALL+1                                                     
  360 G1=0                                                              
      DO 282 I=1,N                                                      
  282 G1=G1+S(I)*GP(I)                                                  
  292 FXP=0.0D0                                                             
      DO 140 I=1,NA                                                     
  140 FXP=FXP+FS(I)*FS(I)                                               
      IF(FXP.GT.UPRLIM) GO TO 157                                       
C                                                                       
C     FIND SATISFACTORY LENGTH FOR STEP ALONG VECTOR S                  
C                                                                       
  138 IF(ND.NE.0) GO TO 4138                                            
      IF(FXP.LT.FX) GO TO 159                                           
      IF((G1.GT.0.0D0).AND.(NTURN.EQ.1)) GO TO 147                          
  157 ALPHA=ALPHA*W                                                     
      AA=AA*W                                                           
      GO TO 148                                                         
  147 ALPHA=FCUBIC(FX,G0,FXP,G1,ALPHA)                                  
  148 NTURN=NTURN+1                                                     
      GO TO 120                                                         
  159 IF(G1.GE.0.0D0) GO TO 150                                             
      AA=AA/W                                                           
      AA=DMIN1(1.D0,AA)                                                 
      GO TO 150                                                         
C                                                                       
C     MAIN ITERATIVE SECTION : LINEAR SEARCH TO MINIMUM                 
C                                                                       
 4138 IF(FXP.LT.FX) GO TO 4159                                          
      IF(G1.LT.0.0D0) GO TO 4869                                            
      IKEY=1                                                            
      GO TO 4140                                                        
 4150 IF(FZ.GE.FX) GO TO 4869                                           
      GO TO 4160                                                        
 4869 ALPHA=ALPHA*W                                                     
      AA=AA*W                                                           
      NTURN=NTURN+1                                                     
      GO TO 120                                                         
 4159 IF(G1.GE.0.0D0) GO TO 4871                                            
      AA=AA/W                                                           
      AA=DMIN1(1.D0,AA)                                                 
      GO TO 150                                                         
 4871 IKEY=2                                                            
      GO TO 4140                                                        
 4170 IF(FZ.GE.FXP) GO TO 4877                                          
 4160 FXP=FZ                                                            
      DO 4867 I=1,N                                                     
      XP(I)=ETA(I)                                                      
 4867 GP(I)=RHO(I)                                                      
      GO TO 150                                                         
 4877 AA=AA*W                                                           
      GO TO 150                                                         
 4140 ALPHA=FCUBIC(FX,G0,FXP,G1,ALPHA)                                  
      DO 4873 I=1,N                                                     
 4873 ETA(I)=X(I)+ALPHA*X(I)                                            
      CALL CALGRA(NA,N,FS,ETA,FJACOB,RHO)                               
      NCALL=NCALL+1                                                     
      FZ=0                                                              
      DO 4875 I=1,NA                                                    
 4875 FZ=FZ+FS(I)*FS(I)                                                 
      IF(FZ.LT.UPRLIM) GO TO 4220                                       
      AA=AA*W                                                           
      GO TO 4869                                                        
 4220 GO TO (4150,4170),IKEY                                            
C                                                                       
C     ITERATION COMPLETE                                                
C     GO ON TO CALCULATE S,UPDATE METRIC AND TEST FOR CONVERGENCE       
C                                                                       
  150 SX=SXNEW                                                          
      T=TNEW                                                            
      DO 160 I=1,N                                                      
  160 Y(I)=GP(I)-G(I)                                                   
      SSX=SX*ALPHA                                                      
      FX3=FX2                                                           
      FX2=FX1                                                           
      FX1=FXP                                                           
C                                                                       
C     UPDATE METRIC MATRIX                                              
C                                                                       
      SY=0.0D0                                                              
      EY=0.0D0                                                              
      DO 180 I=1,N                                                      
      ETA(I)=0.0D0                                                          
      DO 170 J=1,N                                                      
  170 ETA(I)=ETA(I)+H(I,J)*Y(J)                                         
      SY=SY+SIG(I)*Y(I)                                                 
  180 EY=EY+ETA(I)*Y(I)                                                 
      IF(ITYP.EQ.0) GO TO 2280                                          
      DO 280 I=1,N                                                      
      DO 280 J=1,N                                                      
  280 H(I,J)=H(I,J)+SIG(I)*SIG(J)/SY-ETA(I)*ETA(J)/EY                   
      GO TO 90                                                          
 2280 Z=1.0D0+EY/SY                                                         
      DO 430 I=1,N                                                      
      DO 430 J=1,N                                                      
  430 H(I,J)=H(I,J)-(SIG(I)*ETA(J)+ETA(I)*SIG(J)-Z*SIG(I)*SIG(J))/SY    
C                                                                       
C     CALCULATES DERIVATIVE IN DIRECTION S=-H*G                         
C     CHECKS THAT THIS IS NEGATIVE (H POSITIVE DEFINITE)                
C     IF NOT RESETS H TO UNIT MATRIX OR QUITS(SECOND TIME)              
C                                                                       
   90 SXNEW=0.0D0                                                           
      TNEW=0.0D0                                                            
      G0=0.0D0                                                              
      DO 136 I=1,N                                                      
      S(I)=0.0D0                                                            
      DO 100 J=1,N                                                      
  100 S(I)=S(I)-H(I,J)*GP(J)                                            
      TES=DABS(S(I))                                                    
      IF(TES.LE.SXNEW) GO TO 136                                        
      SXNEW=TES                                                         
      TNEW=DABS(XP(I))                                                  
  136 G0=G0+S(I)*GP(I)                                                  
      IF(G0.LE.0.0D0) GO TO 300                                             
  824 IF(ISCREW.GT.2) GO TO 850                                         
      DO 855 I=1,N                                                      
      DO 855 J=1,N                                                      
      H(I,J)=0.0D0                                                          
      IF(I.NE.J) GO TO 855                                              
      H(I,I)=1.0D0                                                          
  855 CONTINUE                                                          
      ISCREW=ISCREW+1                                                   
      GO TO 90                                                          
  850 ISIG=-5                                                           
      WRITE(6,2900)NIT,G0                                                 
 2900 FORMAT(1H0,2X,'--- IT ',I4,' G0 = ',G12.6/ 3X,'    METRIC CANNOT B
     *E POSITIVE DEFINITE ___')                                         
      GO TO 250                                                         
C                                                                       
C      CONVERGENCE TEST                                                 
C      APPLIED ONLY AT 3 RD AND LATER ITERATIONS                        
C                                                                       
  300 DO 210 I=1,N                                                      
      G(I)=GP(I)                                                        
  210 X(I)=XP(I)                                                        
      FX=FXP                                                            
      IF(G0.EQ.0.0D0) GO TO 250                                             
      IF(ISIG.EQ.10) GO TO 70                                           
      IF(NIT.LE.3) GO TO 220                                            
      T1=DABS(FX3-FX1)                                                  
      T2=EPS*DABS(FX1)                                                  
      IF(((SX.LE.EPS*T).AND.(SXNEW.LE.EPS*TNEW)).OR.(T1.LE.T2))GO TO 250
      IF(NCALL.LT.MAXFUN) GO TO 220                                     
 215  ISIG=1                                                            
      GO TO 250                                                         
  220 NIT=NIT+1                                                         
      GO TO 70                                                          
C                                                                       
  250 DO 444 I=1,N                                                      
      DO 444 J=1,N                                                      
  444 H(I,J)=2.0D0*H(I,J)                                                   
      RETURN                                                            
      END                                                               
C                                                                       
C     -------------------------------------------------------------     
C                                                                       
      SUBROUTINE QMTRIX(N,RJTJ,D,RJTF,S,DELTA,QPARM)                    
C                                                                       
CDC      ALLOCLEVEL RJTJ,D,RJTF,S,DELTA                                 
      INTEGER I,J,N                                                     
      REAL*8 RJTJ(N,N),D(N),RJTF(N),S(N),DELTA(N)                       
      REAL*8 QPARM                                                      
      REAL*8 DETERM                                                     
C                                                                       
C       THIS ROUTINE SOLVES THE MARQUARDT MATRIX EQUATION               
C     ((J TRANSPOSE) * J + QPARM * S) * DELTA = -(J  TRANSPOSE) * F     
C     WHERE  J  IS THE JACOBIAN MATRIX, S  IS THE SCALE CONTROL VECTOR, 
C     DELTA  IS THE CORRECTION VECTOR AND  F  IS THE FUNCTION VECTOR.   
C     THE MARQUARDT PARAMETER  QPARM  IS INCREASED IN VALUE IF THE      
C     ORIGINAL MARQUARDT MATRIX IS SINGULAR.                            
C                                                                       
C       THIS ROUTINE REQUIRES THE SUBROUTINE  SYMINV  OR ITS EQUIVALENT.
C     SYMINV  SOLVES THE MATRIX EQUATION  A * X = B  BY GAUSSIAN        
C     ELIMINATION. THE MATRIX  A  IS SYMMETRIC. THE DETERMINANT AND THE 
C     INVERSE OF  A  ARE ALSO OBTAINED. ONLY THE UPPER TRIANGULAR       
C     PORTION OF  A  IS USED AND ONLY THE UPPER TRIANGULAR PORTION OF   
C     (A INVERSE)  IS STORED.                                           
C                                                                       
C     ON INPUT                                                          
C                                                                       
C        N  IS THE DIMENSION OF THE VECTORS  D, DELTA, RJTF  AND  S     
C           AND THE SQUARE MATRIX IN  RJTJ.                             
C                                                                       
C        RJTJ  CONTAINS THE LOWER TRIANGLE OF  (J TRANSPOSE) * J        
C           STORED BY ROWS                                              
C                                                                       
C        D  CONTAINS THE DIAGONAL ELEMENTS OF  RJTJ.                    
C                                                                       
C        RJTF  CONTAINS THE ELEMENTS OF  (J TRANSPOSE) * F.             
C                                                                       
C        S  CONTAINS THE MARQUARDT MULTIPLIERS.                         
C                                                                       
C        QPARM  IS THE MARQUARDT PARAMETER.                             
C                                                                       
C     ON OUTPUT                                                         
C                                                                       
C        RJTJ  CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF THE        
C           MARQUARDT MATRIX STORED BY ROWS. THE ELEMENTS IN  RJTJ      
C           BELOW THE MAIN DIAGONAL HAVE NOT BEEN ALTERED BY THE        
C           ROUTINE.                                                    
C                                                                       
C        DELTA  CONTAINS THE NEGATIVE CORRECTION VECTOR.                
C                                                                       
C        QPARM  IS THE UPDATED MARQUARDT PARAMETER.                     
C                                                                       
C     -------------------------------------------------------------     
C                                                                       
C     ASSEMBLE THE MARQUARDT MATRIX. IN THIS SECTION THE DIAGONAL       
C       ELEMENTS OF  (J  TRANSPOSE) * J  ARE AUGMENTED BY  QPARM * S    
C       SO THAT WHEN FINISHED                                           
C                                                                       
C       RJTJ  -  CONTAINS THE LOWER TRIANGLE OF THE MARQUARDT MATRIX    
C                                                                       
C       ALSO                                                            
C                                                                       
C       DELTA(I)  -  CONTAINS THE VECTOR  RJTF                          
C                                                                       
  100 DO 110 I = 1, N                                                   
      DELTA(I) = RJTF(I)                                                
      RJTJ(I,I) = D(I) + QPARM * S(I)                                   
  110 CONTINUE                                                          
C                                                                       
C     SOLVE THE MARQUARDT MATRIX EQUATION. WHEN FINISHED                
C                                                                       
C       RJTJ     - CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF THE    
C                  MARQUARDT MATRIX. THE ELEMENTS BELOW THE MAIN        
C                  DIAGONAL HAVE NOT BEEN ALTERED                       
C       DELTA(I) - CONTAINS THE NEGATIVE CORRECTION VECTOR              
C                                                                       
      DO 120 I = 1, N                                                   
      DO 120 J = I, N                                                   
      RJTJ(I,J) = RJTJ(J,I)                                             
  120 CONTINUE                                                          
C                                                                       
      CALL SYMINV(RJTJ,N,DELTA,1,DETERM,N)                              
C                                                                       
C     IF THE DETERMINANT OF THE MARQUARDT MATRIX IS ZERO, QPARM  IS     
C       INCREASED IN VALUE AND THE MARQUARDT MATRIX IS REEVALUATED      
C                                                                       
      IF (DETERM .NE. 0.0D0) GO TO 300                                  
      QPARM = 2.0D0 * QPARM                                             
      IF (QPARM .EQ. 0.0D0) QPARM = 1.0D0                               
      GO TO 100                                                         
C                                                                       
  300 RETURN                                                            
C     ----- LAST CARD OF SUBROUTINE  QMTRIX  ----                       
      END 
      
      SUBROUTINE QUAVER(DIAG,EPS,TOL,X,XP,G,GP,S,SIG,Y,ETA,RHO,H,FS, 
     * FJACOB,NA,N,IPRINT,MAXFUN,ISIG)                                      

C     QUASI-NEWTON MINIMIZATION ; DX = CONST * (-H * G)                 
C     WHERE G IS LOCAL DERIVATIVE ,H IS INVERSE SECOND-DERIVATIVE MATRIX
C     H IS APPROXIMATED BY JT*J WITH J JACOBIAN OF RESIDUALS            
C                                                                       
C     MHM DECEMBER 29 1975                                              
C                                                                       
C     MODIFIED TO WORK IN TERMS OF PSEUDO-INVERSES                      
C     MHM JANUARY 29 1976                                               
C                                                                       
C     DIAG : STEPSIZE CONTROL                                           
C            W = /DIAG/ OR /DIAG/**-1 WHICHEVER IS LESS THAN 1          
C            IF DIAG=0 OR /DIAG/=1 W=1/3                                
C            WHEN CHISQ INCREASES REPLACE STEP BY W*STEP                
C     EPS  : RELATIVE ACCURACY REQUESTED                                
C     TOL  : SINGULAR VALUES OF J SET 0 IF LT. TOL*(MAXIMUM VAL)        
C     X    : VECTOR OF VARIABLE-VALUES  ( DIMENSION N)                  
C            INITIALLY CONTAINS GUESS AT POSITION OF MINIMUM            
C            RETURNS VALUES AT LAST ITERATION ; THESE WILL BE EITHER    
C            VALUES AT MIN TO ACCURACY EPS, OR A BEST APPROACH TO THESE 
C            IN MAXFUN CALLS OF SUBROUTINE CALGRA                       
C     XP,G,GP,S,SIG,Y,RHO : 7 WORK ARRAYS OF LENGTH N                   
C     ETA  : VECTOR (DIMENSION N) OF SINGULAR VALUES OF JACOBIAN        
C     H    : N BY N MATRIX - RETURNS VARIANCE-COVARIANCE MATRIX         
C     FS   : VECTOR (DIMENSION NA) OF WEIGHTED RESIDUALS                
C     FJACOB: NA BY N JACOBIAN MATRIX DF(I)/DX(J) OF RESIDUALS          
C     N    : NO OF PARAMETERS                                           
C     NA   : NO OF TERMS IN CHISQ SUM                                   
C     IPRINT: PRINT CONTROL (USES MINMON)                               
C     MAXFUN : UPPER LIMIT ON NUMBER OF FUNCTION CALLS                  
C     ISIG : AT RETURN ISIG = 0 MEANS NORMAL COMPLETION                 
C                      ISIG = 1 - STOPPED AFTER MAXFUN FUNCTION CALLS   
C                      ISIG =-5 - ABNORMAL CONDITIONS ENCOUNTERED       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,SIG,Y,ETA,RHO,H,FS,FJACOB               
C                                                                       
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),SIG(N),Y(N),ETA(N),RHO(N), 
     * H(N,N),FS(NA),FJACOB(NA,N)                                          
C                                                                       
C     INITIAL FUNCTION-CALL                                             
C                                                                       
      DG=DABS(DIAG)                                                     
      IF(DG.LT.1.0D0) W=DG                                                  
      IF((DG.EQ.0.0D0).OR.(DG.EQ.1.0D0)) W=1.0D0/3.0D0                                
      IF(DG.GT.1.0D0) W=1.0D0/DG                                                
      WRITE(6,4050)
 4050 FORMAT(1H0,2X,'++ QUASI-NEWTON MINIMIZATION (PSEUDO-INVERSE) ++') 
      WRITE(6,4450)W,TOL                                                  
 4450 FORMAT(1H ,5X,'STEP-REDUCTION FACTOR = ',G12.6/ 6X,'MAXIMUM RATIO 
     *OF LARGEST TO SMALLEST SING.VAL. = ',G12.6)                       
      ISIG=10                                                           
      NIT=0                                                             
      AA=1.0D0                                                              
      UPRLIM=1.D+18                                                     
      NCALL=1                                                           
      DO 5 I=1,N                                                        
    5 XP(I)=X(I)                                                        
      FXP=CALGRA(NA,N,FS,XP,FJACOB,GP)                                  
      FX1=0.0D0                                                             
      FX2=0.0D0                                                             
      FX3=0.0D0                                                             
C                                                                       
C     DETERMINE DIRECTION FOR NEXT STEP                                 
C     BY SOLUTION OF LINEAR LSQ PROBLEM WITH JACOBIAN AND RESIDUALS     
C                                                                       
   30 DO 50 I=1,N                                                       
      DO 50 IA=1,NA                                                     
   50 FJACOB(IA,I)=FJACOB(IA,I)*XP(I)                                   
      CALL MINFIT(NA,NA,N,FJACOB,ETA,1,FS,IERR,RHO)                     
      IF(IERR.EQ.0) GO TO 32                                            
      WRITE(6,4500)NIT,IERR                                               
 4500 FORMAT(1H0,2X,'-- AT IT NO.',I3,' TROUBLE WITH ',I2, 'TH SINGULAR 
     *VAL OF J --')                                                     
      ISIG=-5                                                           
      GO TO 250                                                         
   32 TOP=0.0D0                                                             
      DO 55 J=1,N                                                       
   55 IF(ETA(J).GT.TOP) TOP=ETA(J)                                      
      THRESH=TOL*TOP                                                    
      NCHOP=0                                                           
      DO 34 I=1,N                                                       
      S(I)=0.0D0                                                            
      Y(I)=ETA(I)                                                       
      IF(ETA(I).GE.THRESH) GO TO 34                                     
      Y(I)=0.0D0                                                            
      NCHOP=NCHOP+1                                                     
   34 CONTINUE                                                          
      DO 60 K=1,N                                                       
      IF(Y(K).EQ.0) GO TO 60                                            
      DO 36 I=1,N                                                       
   36 S(I)=S(I)-FJACOB(I,K)*FS(K)/Y(K)                                  
   60 CONTINUE                                                          
      DO 80 J=1,N                                                       
   80 S(J)=XP(J)*S(J)                                                   
      IF(NCHOP.GT.0) WRITE(6,4570)NCHOP                                   
 4570 FORMAT(1H,2X,'THE ',I2,' SMALLEST SINGULAR VALUES ','SET TO ZERO')                                                                
      GO TO 90                                                          
C                                                                       
C     PRINT-OUT OF INFORMATION AT START OF ITERATION                    
C                                                                       
   70 ISIG=0                                                            
      IF(IPRINT.GE.1) CALL MINMON(NIT,NCALL,NCALL,N,NA,X,FX,FJACOB)     
      IF(IPRINT.EQ.2) WRITE(6,4550)(ETA(I),I=1,N)                         
 4550 FORMAT(1H0,2X,'SINGULAR VALUES OF J '/ ( 3X, 10G13.6 ) )          
C                                                                       
C     MAIN ITERATIVE SECTION  : INCOMPLETE LINEAR SEARCH                
C                                                                       
  132 IF(AA.LT.1.0D0) AA=W                                                  
      ALPHA=AA                                                          
      NTURN=1                                                           
  120 DO 130 I=1,N                                                      
      SIG(I)=ALPHA*S(I)                                                 
  130 XP(I)=X(I)+SIG(I)                                                 
      FXP=CALGRA(NA,N,FS,XP,FJACOB,GP)                                  
      IF(FXP.GT.UPRLIM) GO TO 157                                       
      NCALL=NCALL+1                                                     
  360 G1=0.0D0                                                              
      DO 282 I=1,N                                                      
  282 G1=G1+S(I)*GP(I)                                                  
C                                                                       
C     FIND SATISFACTORY LENGTH FOR STEP ALONG VECTOR S                  
C                                                                       
  138 IF(FXP.LT.FX) GO TO 159                                           
      IF((G1.GE.0.0D0).AND.(NTURN.EQ.1)) GO TO 147                          
  157 ALPHA=ALPHA*W                                                     
      AA=AA*W                                                           
      GO TO 148                                                         
  147 ALPHA=FCUBIC(FX,G0,FXP,G1,ALPHA)                                  
      AA=AA*W                                                           
  148 NTURN=NTURN+1                                                     
      IF(NTURN.GT.5) GO TO 149                                          
      GO TO 120                                                         
  159 IF((G1.GE.0.0D0).OR.(NTURN.GT.1)) GO TO 30                            
      AA=AA/W                                                           
      AA=DMIN1(1.D0,AA)                                                 
      GO TO 30                                                          
  149 ISIG=-5                                                           
      WRITE(6,3200)
 3200 FORMAT(1H0,2X,'--- F DOES NOT DECREASE ENOUGH IN 5 CHOPS ---')    
      GO TO 250                                                         
C                                                                       
C     ITERATION COMPLETE                                                
C     GO ON TO TEST FOR CONVERGENCE                                     
C                                                                       
   90 FX3=FX2                                                           
      FX2=FX1                                                           
      FX1=FXP                                                           
      SX=0.0D0                                                              
      T=0.0D0                                                               
      G0=0.0D0                                                              
      DO 136 I=1,N                                                      
      TES=DABS(S(I))                                                    
      IF(TES.LE.SX) GO TO 136                                           
      SX=TES                                                            
      T=DABS(XP(I))                                                     
  136 G0=G0+S(I)*GP(I)                                                  
      IF(G0.LE.0.0D0) GO TO 300                                             
  850 ISIG=-5                                                           
      WRITE(6,2900)NIT,G0                                                 
 2900 FORMAT(1H0,2X,'--- IT ',I4,' G0 = ',G12.6/ 3X,'    METRIC CANNOT B
     *E POSITIVE DEFINITE ___')                                         
      GO TO 250                                                         
C                                                                       
C     CONVERGENCE TEST                                                  
C     NOT APPLIED ON FIRST TWO PASSES                                   
C                                                                       
  300 DO 210 I=1,N                                                      
      G(I)=GP(I)                                                        
  210 X(I)=XP(I)                                                        
      FX=FXP                                                            
      IF(G0.EQ.0.0D0) GO TO 250                                             
      IF(ISIG.EQ.10) GO TO 70                                           
      IF(NIT.LT.1) GO TO 220                                            
      IF((SX.LE.EPS*T).OR.(DABS(FX3-FX1).LE.EPS*DABS(FX1)))GO TO 250    
      IF(NCALL.LT.MAXFUN) GO TO 220                                     
      ISIG=1                                                            
      GO TO 250                                                         
  220 NIT=NIT+1                                                         
      GO TO 70                                                          
C                                                                       
C     CALCULATE FINAL VARIANCE-COVARIANCE MATRIX AND QUIT               
C                                                                       
  250 DO 65 I=1,N                                                       
      DO 65 J=1,N                                                       
      H(I,J)=0.0D0                                                          
      DO 40 K=1,N                                                       
      IF(ETA(K).LT.1.D-14) GO TO 40                                     
      H(I,J)=H(I,J)+FJACOB(I,K)*FJACOB(J,K)/(ETA(K)*ETA(K))             
   40 CONTINUE                                                          
      H(I,J)=H(I,J)*X(I)*X(J)                                           
      IF(I.EQ.J) GO TO 65                                               
      H(J,I)=H(I,J)                                                     
   65 CONTINUE                                                          
      NIT=NIT+1                                                         
      IF(IPRINT.GE.1) CALL MINMON(NIT,NCALL,NCALL,N,NA,X,FX,FJACOB)     
      IF(IPRINT.EQ.2) WRITE(6,4550)(ETA(I),I=1,N)                         
      RETURN                                                            
      END 
      
      SUBROUTINE SDIAG (NM, N, EPS, ILIMIT, A, D, X, E)                 
C                                                                       
C     DIAGONALIZATION OF SYMMETRIC MATRICES                             
C                                                                       
C     A - (NM X NM) - SYMMETRIC INPUT MATRIX IN EXPANDED FORM           
C     D - (N) - WILL GET THE ARRAY OF N EIGENVALUES                     
C     X - (NM X NM) - WILL GET THE N X N MATRIX OF EIGEVECTORS.         
C         X MAY BE THE SAME AS A.                                       
C     E - (N) - WORK ARRAY.                                             
C                                                                       
C     2/10/77 - SPECIAL CASE FOR 1 X 1 MATRICES - S. P.                 
C     12/28/79 - CFR%F - RPG                                            
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL A,D,X,E                                             
      DIMENSION A(NM,NM),D(NM),X(NM,NM),E(NM)                           
C                                                                       
      IF ( N .LE. 1 ) GO TO 500                                         
      DO 100 I = 1,N                                                    
      DO 100 J = 1,I                                                    
  100 X(I,J) = A(I,J)                                                   
C ----HOUSEHOLDER REDUCTION                                             
      NM1 = N - 1                                                       
      DO 101 NN = 1,NM1                                                 
      I = N +1 - NN                                                     
      M = I - 2                                                         
      F = X(I,I -1)                                                     
      G = F                                                             
      S = 0.0D0                                                             
      IF (M.LT.1) GOTO 99                                               
      DO 102 K = 1,M                                                    
  102 S = S + X(I,K)*X(I,K)                                             
   99 D(I) = S                                                          
      IF (S.EQ.0.0D0) GOTO 103                                             
      M = M + 1                                                         
      S = S + F * F                                                     
      G = DSQRT (S)                                                     
      IF (F.GE.0.0D0) G=-G                                                
      RH = S - F * G                                                    
      X(I,I-1) = F - G                                                  
      F = 0.0D0                                                             
      IF (M .LT. 1) GOTO 103                                            
      DO 107 J = 1,M                                                    
      X(J,I) = X(I,J)/ RH                                               
      S = 0.0D0                                                             
      DO 108 K = 1,J                                                    
  108 S = S + X(J,K) * X(I,K)                                           
      JP1 = J + 1                                                       
      IF (M.LT.JP1) GOTO 98                                             
      DO 109 K = JP1, M                                                 
  109 S = S + X(K,J) * X(I,K)                                           
   98 E(J) = S/RH                                                       
  107 F = F + S * X(J,I)                                                
   89 F = F / (RH + RH)                                                 
      DO 110 J = 1,M                                                    
  110 E(J)= E(J) - F * X(I,J)                                           
      DO 111 J = 1,M                                                    
      DO 111 K = 1,J                                                    
  111 X(J,K) = X(J,K) - X(I,J) * E(K) -X(I,K) * E(J)                    
C-----END  S NOT EQUAL TO ZERO                                          
  103 E(I-1) = G                                                        
  101 CONTINUE                                                          
C-----END I LOOP                                                        
C-----ACCUMULATION  OF TRANSFORMATION  MATRICES                         
      D(1) = 0.0D0                                                          
      E(N)= 0.0D0                                                           
      B = 0.0D0                                                             
      F = 0.0D0                                                             
      DO 200 I = 1,N                                                    
      M = I - 1                                                         
      IF (D(I).EQ.0.0D0)GOTO 201                                           
      DO 202 J = 1,M                                                    
      S = 0.0D0                                                             
      DO 203 K = 1,M                                                    
  203 S = S + X(I,K) * X(K,J)                                           
      DO 202 K = 1,M                                                    
C-----THE FOLLOWING NONSENSE STATEMENT IS NEEDED TO COMBAT COMPILER BUG 
      L=K+I-K                                                           
  202 X(K,J) = X(K,J) - S * X(K,L)                                      
C-----END  J LOOP                                                       
  201 D(I) = X(I,I)                                                     
      X(I,I) = 1.0D0                                                      
      IF (M.LT. 1 ) GOTO 200                                            
      DO 204 J = 1,M                                                    
      X(I,J) = 0.0D0                                                        
  204 X(J,I) = 0.0D0                                                        
  200 CONTINUE                                                          
C-----END I                                                             
C-----DIAGONALIZATION  OF THE TRIDIAGONAL MATRIX                        
      DO 205 L = 1,N                                                    
      J = 0                                                             
      RH = EPS * (DABS (D(L)) + DABS (E(L)))                            
      IF (RH .LE. B) GOTO 206                                           
      B = RH                                                            
C-----TEST FOR SPLITTING                                                
  206 DO 207 MM = L,N                                                   
      M = MM                                                            
      IF (DABS (E(M)) .LT. B) GOTO 208                                  
  207 CONTINUE                                                          
C-----LABEL 208 IS TEST FOR CONVERGENCE                                 
  208 IF (M .EQ. L) GOTO 209                                            
C-----LABEL 209 IS CONVERGENCE                                          
C-----LABEL 210 IS SHIFT                                                
  210 P =0.5D0*(D(L+1)- D(L))                                             
      R = DSQRT (P * P + E(L)* E(L))                                    
      RH = D(L) + P - R                                                 
      IF (P .LT.0.0D0) RH=RH+R+R                                           
      DO 213 I = L,N                                                    
  213 D(I) = D(I) - RH                                                  
      F = F + RH                                                        
C-----QR TRANSFORMATION                                                 
      J = J + 1                                                         
      P = D(M)                                                          
      C = 1.0D0                                                           
      S = 0.0D0                                                             
      MML = M - L                                                       
      DO 214 II = 1,MML                                                 
      I = M - II                                                        
      R = DSQRT (P*P + E(I) * E(I))                                     
      G = C * E(I)                                                      
      RH = C * P                                                        
      E(I+1) = S*R                                                      
      C = P/R                                                           
      S = E(I)/ R                                                       
      D(I+1) = RH + S * (C * G + S * D(I))                              
      DO 214 K = 1,N                                                    
      RH = X(K,I+1)                                                     
      P = C * D(I) - S*G                                                
      X(K,I+1) = X(K,I) * S + RH * C                                    
  214 X(K,I) = X(K,I) * C - RH * S                                      
C-----END K LOOP AND END I LOOP                                         
  215 E(L) = S * P                                                      
      D(L)= C * P                                                       
      IF (DABS(E(L)) .LT. B) GOTO 209                                   
      IF (J .LT. ILIMIT) GOTO 210                                       
      GOTO 1000                                                         
C-----1000 IS ERROR RETURN INDICATING FAILURE OF QR TO CONVERGE         
  209 D(L) = D(L)+ F                                                    
  205 CONTINUE                                                          
C-----END L LOOP                                                        
C-----ORDERING OF EIGENVALUES                                           
      DO 300 I = 1,N                                                    
      K = I                                                             
      P = D(I)                                                          
      IP1 = I + 1                                                       
      IF (N .LT. IP1) GOTO 301                                          
      DO 302 J = IP1 , N                                                
      IF (D(J) .GE. P) GOTO 302                                         
      K = J                                                             
      P = D(J)                                                          
  302 CONTINUE                                                          
C-----END J LOOP                                                        
  301 IF (K.EQ.I) GOTO 300                                              
      D(K) = D(I)                                                       
      D(I) = P                                                          
      DO 304 J = 1,N                                                    
      P = X(J,I)                                                        
      X(J,I) = X(J,K)                                                   
  304 X(J,K) = P                                                        
  300 CONTINUE                                                          
C-----END J,END K NOT EQUAL I, I, SDIAG                                 
      GOTO 1001                                                         
C                                                                       
C     SPECIAL CASE FOR 1 X 1                                            
C                                                                       
 500  D(1) = A(1,1)                                                     
      X(1,1) = 1                                                        
      RETURN                                                            
C                                                                       
 1000 WRITE (6,1002)                                                    
 1002 FORMAT (1H0,25HFAILURE OF QR TO CONVERGE)                         
 1001 RETURN                                                            
      END                                                               
C                                                                       
C     ----------------------------------------------------------------- 
C                                                                       
      SUBROUTINE SIP(A,I,J,B,K,L,ALPHA,SUM,M)                           
C                                                                       
CDC      ALLOCLEVEL A,B                                                 
      INTEGER I,J,K,L,M,N                                               
      REAL*8 A(1),B(1)                                                  
      REAL*8 ALPHA,SUM                                                  
C                                                                       
C       THIS ROUTINE SUMS THE PRODUCTS OF ELEMENTS  A(I), A(J),         
C     A(I+(N-1)*(J-I))  IN  A  AND THE CORRESPONDING ELEMENTS  B(K),    
C     B(L),  B(K+(N-1)*(L-K))  IN  B, FOR N = 3,4,...,M, WITH  ALPHA    
C                                                                       
C       THE ARRAYS  A  AND  B  MUST BE OF DIMENSION AT LEAST            
C     I+(M-1)*(J-I)  AND  K+(M-1)*(L-K)  RESPECTIVELY.                  
C                                                                       
C     ----------------------------------------------------------------  
C                                                                       
      SUM = ALPHA + A(I) * B(K)                                         
      IF (M .LE. 1) GO TO 110                                           
      SUM = SUM + A(J) * B(L)                                           
      IF (M .LE. 2) GO TO 110                                           
C                                                                       
      DO 100 N = 3, M                                                   
      SUM = SUM + A(I+(N-1)*(J-I))*B(K+(N-1)*(L-K))                     
  100 CONTINUE                                                          
C                                                                       
  110 RETURN                                                            
C     ------ LAST CARD OF SUBROUTINE  SIP  ------                       
      END 
      
      SUBROUTINE SV01A (M,N,VAR,IVDIM,IVPAR, W)                         
C                                                                       
C     10/11/75 - W ADDED TO ARG. LIST - S.P.                            
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL VAR,W                                               
      DIMENSION VAR(IVDIM,IVDIM), W(1)                                  
      MN = M + N                                                        
      IMIN = MN + N                                                     
      GO TO (1,5), IVPAR                                                
  1   IMAX = N*(MN+2)                                                   
      DO 2 I = 1,N                                                      
      IMIN = IMIN + 1                                                   
      DO 3 J = I,N                                                      
      KI = J - I                                                        
      SUM = 0.0D0                                                         
      DO 4 K = IMIN,IMAX,MN                                             
      KK = K + KI                                                       
      SUM = SUM + W(K)*W(KK)                                            
  4   CONTINUE                                                          
      VAR(I,J) = SUM                                                    
      VAR(J,I) = SUM                                                    
  3   CONTINUE                                                          
  2   CONTINUE                                                          
      RETURN                                                            
  5   IBAS = (N+1) * (MN + 1)                                           
      DO 6 I = 1,N                                                      
      JCUR = IMIN + I                                                   
      DO 7 J = 1,N                                                      
      W(J) = 0.0D0                                                        
  7   CONTINUE                                                          
      IDEL = IBAS                                                       
      DO 8 J = 1,N                                                      
      W(J) = W(J) + W(IDEL)*W(JCUR)                                     
      IF (J.EQ.N) GO TO 8                                               
      KMIN = J + 1                                                      
      KCUR = JCUR                                                       
      DO 9 K = KMIN,N                                                   
      IDEL = IDEL + 1                                                   
      KCUR = KCUR + MN                                                  
      W(J) = W(J) + W(IDEL) * W(KCUR)                                   
      W(K) = W(K) + W(IDEL) * W(JCUR)                                   
 9    CONTINUE                                                          
      IDEL = IDEL + 1                                                   
      JCUR = JCUR + MN                                                  
  8   CONTINUE                                                          
      DO 10 J = I,N                                                     
      KCUR = IMIN + J                                                   
      SUM = 0.0D0                                                         
      DO 11 K = 1,N                                                     
      SUM = SUM + W(K) * W(KCUR)                                        
      KCUR = KCUR + MN                                                  
 11   CONTINUE                                                          
      VAR(I,J) = SUM                                                    
      VAR(J,I) = SUM                                                    
  10  CONTINUE                                                          
  6   CONTINUE                                                          
      RETURN                                                            
      END 
      
      SUBROUTINE SYMINV(A,N,B,M,DETERM,NMAX)                            
C                                                                       
C     SYMMETRIC MATRIX INVERSION (OPERATES ON UPPER TRIANGLE)           
C                                                                       
C     PTOLEMY VERSION FOR USE BY LMCHOL                                 
C                                                                       
C       1)  ONLY ONE INHOMOGENEOUS TERM.                                
C       2)  PIVOT MADE A LOGICAL*1 ARRAY.                               
C       3)  DETERMINANT CALCULATION REMOVED; LAST PIVOT ELEMENT IS      
C           RETURNED AS DETERMINANT.                                    
C       4)  LIMITED TO  25 X 25  MATRICES.                              
C                                                                       
C     11/5/75 - ABOVE CHANGES FOR PTOLEMY                               
C     12/28/79 - CDC TO CNI - RPG                                       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL A,B                                                 
C                                                                       
CIB      LOGICAL*1 PIVOT                                                
      LOGICAL PIVOT                                                     
      DIMENSION A(NMAX,N), B(NMAX)                                      
      COMMON /F152/ P(25), Q(25), PIVOT(25)                             
C                                                                       
C     INITIALIZATION                                                    
C                                                                       
      DO 20 I=1,N                                                       
      PIVOT(I)=.FALSE.                                                  
   20 CONTINUE                                                          
C                                                                       
C     GRAND LOOP                                                        
C                                                                       
      DO 550 I=1,N                                                      
C                                                                       
C     SEARCH FOR PIVOT ELEMENT                                          
C                                                                       
      DETERM=0.0D0                                                        
      DO 100 J=1,N                                                      
      IF (PIVOT(J)) GO TO 100                                           
      TEMP=DABS(A(J,J))                                                 
      IF (TEMP.LT.DETERM) GO TO 100                                     
      L=J                                                               
      DETERM=TEMP                                                       
  100 CONTINUE                                                          
      DETERM=A(L,L)                                                     
C                                                                       
C     IF PIVOT ELEMENT IS ZERO, PROGRAM FAILS                           
C                                                                       
      IF (DETERM.EQ.0.0D0) RETURN                                         
C                                                                       
C     PREPARATION OF ELIMINATION STEP                                   
C                                                                       
      PIVOT(L)= .TRUE.                                                  
      P(L)=1.0D0                                                          
      Q(L)=1.0D0/DETERM                                                    
      A(L,L)=0.0D0                                                        
      QB = B(L)                                                         
      B(L)=0.0D0                                                          
  120 CONTINUE                                                          
C                                                                       
  130 IF (L.EQ.1) GO TO 230                                             
      L1=L-1                                                            
      DO 200 J=1,L1                                                     
      P(J)=A(J,L)                                                       
      Q(J)=P(J)*Q(L)                                                    
      IF ( .NOT. PIVOT(J) ) Q(J)=-Q(J)                                  
  200 A(J,L)=0.0D0                                                        
C                                                                       
  230 IF (L.EQ.N) GO TO 400                                             
      L1=L+1                                                            
      DO 300 K=L1,N                                                     
      P(K)=-A(L,K)                                                      
      Q(K)=P(K)*Q(L)                                                    
      IF ( .NOT. PIVOT(K) ) P(K)=-P(K)                                  
  300 A(L,K)=0.0D0                                                        
C                                                                       
C     ELIMINATION PROPER                                                
C                                                                       
  400 DO 550 J=1,N                                                      
      DO 450 K=J,N                                                      
  450 A(J,K)=A(J,K)+P(J)*Q(K)                                           
  500 B(J)=B(J)+Q(J)*QB                                                 
C                                                                       
  550 CONTINUE                                                          
C                                                                       
      RETURN                                                            
      END                                                               
C                                                                       
C     -------------------------------------------------------------     
C                                                                       
      SUBROUTINE UPQPAR(N,SQS,RJTJ,ACTRED,PRERED,DIRDRV,QPARM,QCPARM)   
C                                                                       
CDC      ALLOCLEVEL SQS,RJTJ                                            
      INTEGER I,II,J,N                                                  
      REAL*8 SQS(N),RJTJ(N,N)                                           
      REAL*8 ACTRED,PRERED,DIRDRV,QPARM,QCPARM                          
      REAL*8 DABS,FNORM,TRACE,Y,Z                                       
C                                                                       
C       THIS ROUTINE UPDATES THE MARQUARDT PARAMETER  QPARM  BASED ON   
C     THE DEGREE OF AGREEMENT BETWEEN THE ACTUAL AND PREDICTED REDUCTION
C     IN THE SUM OF SQUARES  FSQ.                                       
C                                                                       
C     ON INPUT                                                          
C                                                                       
C        N  IS THE DIMENSION OF THE VECTOR  SQS  AND THE SQUARE  MATRIX 
C           IN  RJTJ.                                                   
C                                                                       
C        SQS  CONTAINS  SQRT(S(I)), I = 1,2,...,N, WHERE THE  S(I)  ARE 
C           ELEMENTS OF THE MARQUARDT SCALE CONTROL VECTOR.             
C                                                                       
C        RJTJ  CONTAINS THE UPPER TRIANGLE OF THE INVERSE OF THE        
C           MARQUARDT MATRIX  ((J TRANSPOSE) * J + QPARM * S).          
C                                                                       
C        ACTRED  IS THE ACTUAL REDUCTION IN THE SUM OF SQUARES  FSQ.    
C                                                                       
C        PRERED  IS THE PREDICTED REDUCTION IN THE SUM OF SQUARES  FSQ. 
C                                                                       
C        DIRDRV  IS THE DIRECTIONAL DERIVATIVE ALONG THE CURRENT        
C           CORRECTION VECTOR.                                          
C                                                                       
C        QPARM  IS THE MARQUARDT PARAMETER.                             
C                                                                       
C        QCPARM  IS THE LOWER CUT-OFF LIMIT ON THE MARQUARDT            
C           PARAMETER.                                                  
C                                                                       
C     ON OUTPUT                                                         
C                                                                       
C        QPARM  IS THE UPDATED MARQUARDT PARAMETER.                     
C                                                                       
C        QCPARM  IS THE UPDATED LOWER CUT-OFF LIMIT ON THE              
C           MARQUARDT PARAMETER.                                        
C                                                                       
C     -------------------------------------------------------------     
C                                                                       
C     TEST WHETHER THE RATIO OF THE ACTUAL REDUCTION TO THE PREDICTED   
C       REDUCTION IS .GE. 0.25                                          
C                                                                       
  100 IF (ACTRED .GE.0.25D0*PRERED) GO TO 160                          
C                                                                       
C     ACTRED / PRERED  .LT. 0.25. CALCULATE  NU  SUCH THAT              
C       2 .LE.  NU  .LE. 10. IN THIS SECTION                            
C                                                                       
C       Y  -  CONTAINS 1 / NU                                           
C                                                                       
      Y = 0.5D0                                                           
      Z = 2.0D0 * DIRDRV - ACTRED                                       
      IF (Z .GT. 0.0D0) Y = DIRDRV / Z                                  
      IF (Y .GT. 0.5D0) Y = 0.5D0                                           
      IF (Y .LT. 0.1D0) Y = 0.1D0                                       
C                                                                       
C     TEST WHETHER  QPARM = 0                                           
C                                                                       
      IF (QPARM .NE. 0.0D0) GO TO 150                                   
C                                                                       
C     QPARM = 0. CALCULATE  QCPARM = MAX(1/FNORM),(1/TRACE))  AND SET   
C       QPARM = QCPARM,  Y = 2 * Y. IN THIS SECTION                     
C                                                                       
C       FNORM  -  CONTAINS A SCALED INFINITY NORM OF  RJTJ  INVERSE     
C       TRACE  -  CONTAINS A SCALED TRACE OF  RJTJ  INVERSE             
C                                                                       
      Y = 2.0D0 * Y                                                     
      QPARM = 0.0D0                                                     
      TRACE = 0.0D0                                                     
C                                                                       
      DO 140 I = 1, N                                                   
      TRACE = TRACE + RJTJ(I,I) * SQS(I) * SQS(I)                       
      FNORM = 0.0D0                                                     
C                                                                       
      DO 110 J = 1, I                                                   
      FNORM = FNORM + DABS(RJTJ(J,I)) * SQS(J)                          
  110 CONTINUE                                                          
C                                                                       
      IF (I .EQ. N) GO TO 130                                           
      II = I + 1                                                        
C                                                                       
      DO 120 J = II, N                                                  
      FNORM = FNORM + DABS(RJTJ(I,J)) * SQS(J)                          
  120 CONTINUE                                                          
C                                                                       
  130 FNORM = FNORM * SQS(I)                                            
      IF (FNORM .GT. QPARM) QPARM = FNORM                               
  140 CONTINUE                                                          
C                                                                       
      IF (TRACE .LT. QPARM) QPARM = TRACE                               
      QPARM = 1.0D0 / QPARM                                             
      QCPARM = QPARM                                                    
C                                                                       
C     QPARM = NU * QPARM                                                
C                                                                       
  150 QPARM = QPARM / Y                                                 
      GO TO 200                                                         
C                                                                       
C     ACTRED / PRERED  .GE.  0.25. TEST WHETHER  ACTRED / PRERED        
C       .LE.  0.75                                                      
  160 IF(ACTRED .LE.0.75D0*PRERED) GO TO 200                           
C                                                                       
C     ACTRED / PRERED  .GT.  0.75. SET  QPARM = QPARM / 2. THEN SET     
C       QPARM = 0 IF  QPARM  .LT.  QCPARM                               
C                                                                       
      QPARM = QPARM * 0.5D0                                               
      IF (QPARM .LT. QCPARM) QPARM = 0.0D0                              
  200 RETURN                                                            
C     ----- LAST CARD OF SUBROUTINE  UPQPAR  ----                       
      END 
      
      SUBROUTINE VA02A ( M, N, F, X, TOL, TOL2, ESCALE, IPRINT, MAXFUN, 
     * IRET, W )                                                        
C                                                                       
C     POWELL 1965 LEAST SQUARES FITTER FROM  AMD ARGONNE                
C                                                                       
C     M = NUMBER OF DATA POINTS TO BE FIT TO.                           
C     N = NUMBER OF PARAMETERS TO VARY.                                 
C     F = COMPUTED ARRAY OF FUNCTION VALUES FOR EACH DATA POINT.        
C         DIMENSION IS M.                                               
C     X = PARAMETER ARRAY (DIMENSION N).                                
C     TOL AND TOL2 DETERMIN THE ACCURACY WITH WHICH THE MINIMUM IS      
C     FOUND:  THE FITTER STOPS WHENEVER A STEP IS SUCH THAT             
C     1)  A STEP IS SUCH THAT  DELTA(X(I)) < MAX( TOL*\X(I)\, TOL2 )    
C         FOR EVERY I;  OR                                              
C     2)  THE SUM OF SQUARES FOR 4 CONSEQUTIVE ITERATIONS HAS           
C         CHANGED BY LESS THAN TOL:                                     
C              FSQ(K-4) - FSQ(K) < TOL*FSQ(K)                           
C     3)  THE SUM OF SQUARES BECOMES LESS THAN TOL2.                    
C                                                                       
C     ESCALE - THE MAXIMUM SINGLE STEP IN X(I) WILL BE ESCALE*DELTA.    
C           WHERE   DELTA = MAX( TOL*\X(I)\, TOL2 )                     
C     IPRINT -  0 MEANS DON'T CALL MINMON FOR EACH ITERATION            
C               1 MEANS DO CALL MINMON FOR EACH ITERATION               
C     MAXFUN - MAXIMUM ALLOWED NUMBER OF CALLS TO CALFUN.               
C     IRET - 0  FIT HAS COMPLETED O.K.                                  
C          - 1  RAN OUT OF FUNCTION REFS.                               
C     W - A WORK ARRAY WHOSE DIMENSION MUST BE ATLEAST                  
C           (N+1)*(M+N+1) + (N*(N+1))/2 .                               
C                                                                       
C      8/75  IRET ADDED, A FEW PRINT STATEMSNTS CHANGED. D.G.           
C     10/11/75 - W ADDED TO ARG. LIST - S.P.                            
C     10/17/75 - MORE PRINTOUT AT ST. NO. 5; MORE COMMENTS - S.P.       
C     11/6/75 - USE MINMON FOR PRINTING                                 
C     1/22/76 - CHANGE IRET; ADD TOL AND TOL2; CALFUN IS A FUNCTION     
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,F,W                                               
      DIMENSION F(M), X(N), W(1)                                        
C                                                                       
      ITC = 0                                                           
      IRET=0                                                            
      MPLUSN=M+N                                                        
      KST=N+MPLUSN                                                      
      NPLUS=N+1                                                         
      KINV=NPLUS*(MPLUSN+1)                                             
      KSTORE=KINV-MPLUSN-1                                              
      FF = CALFUN (M,N,F,X)                                             
      MC = 1                                                            
      NN=N+N                                                            
      K=NN                                                              
C                                                                       
      DO 1 I=1,M                                                        
      K=K+1                                                             
      W(K)=F(I)                                                         
    1 CONTINUE                                                          
C                                                                       
      IINV=2                                                            
      K=KST                                                             
C                                                                       
C     THIS IS THE START OF A LOOP THAT IS TERMINATED AT ST. NO. 24      
C                                                                       
      I=1                                                               
  2   DELTA = DMAX1( TOL*DABS(X(I)), TOL2 )                             
      X(I)=X(I)+DELTA                                                   
      FC = CALFUN (M,N,F,X)                                             
      MC = MC+1                                                         
      X(I)=X(I)-DELTA                                                   
      DO 3 J=1,N                                                        
      K=K+1                                                             
      W(K)=0.0D0                                                           
      W(J)=0.0D0                                                           
    3 CONTINUE                                                          
      SUM=0.0D0                                                            
      KK=NN                                                             
      DO 4 J=1,M                                                        
      KK=KK+1                                                           
      F(J)=F(J)-W(KK)                                                   
      SUM=SUM+F(J)*F(J)                                                 
    4 CONTINUE                                                          
      IF (SUM) 5,5,6                                                    
5001  FORMAT('0**** ACCURACY FOR THE',I3,'TH VARIABLE IS UNREASONABLY',
     * ' SMALL' / ' X(I), DELTA =',G20.10,G15.5,10X,'SUM =', G15.5)                                                                 
5     WRITE(6,5001) I, X(I), DELTA, SUM                                 
      DO 8 J=1,M                                                        
      NN=NN+1                                                           
      F(J)=W(NN)                                                        
    8 CONTINUE                                                          
      GO TO 10                                                          
    6 SUM=1.0D0/DSQRT(SUM)                                                 
      J=K-N+I                                                           
      W(J)=DELTA*SUM                                                    
      DO 9 J=1,M                                                        
      K=K+1                                                             
      W(K)=F(J)*SUM                                                     
      KK=NN+J                                                           
      DO 11 II=1,I                                                      
      KK=KK+MPLUSN                                                      
      W(II)=W(II)+W(KK)*W(K)                                            
   11 CONTINUE                                                          
    9 CONTINUE                                                          
      ILESS=I-1                                                         
      IGAMAX=N+I-1                                                      
      INCINV=N-ILESS                                                    
      INCINP=INCINV+1                                                   
   12 IF (ILESS.NE.0) GO TO 14                                          
      W(KINV)=1.0D0                                                        
      GO TO 15                                                          
   14 B=1.0D0                                                              
      DO 16 J=NPLUS,IGAMAX                                              
      W(J)=0.0D0                                                           
   16 CONTINUE                                                          
      KK=KINV                                                           
      DO 17 II=1,ILESS                                                  
      IIP=II+N                                                          
      W(IIP)=W(IIP)+W(KK)*W(II)                                         
      JL=II+1                                                           
      IF (JL-ILESS) 18,18,19                                            
   18 DO 20 JJ=JL,ILESS                                                 
      KK=KK+1                                                           
      JJP=JJ+N                                                          
      W(IIP)=W(IIP)+W(KK)*W(JJ)                                         
      W(JJP)=W(JJP)+W(KK)*W(II)                                         
   20 CONTINUE                                                          
   19 B=B-W(II)*W(IIP)                                                  
      KK=KK+INCINP                                                      
   17 CONTINUE                                                          
      B=1.0D0/B                                                            
      KK=KINV                                                           
      DO 21 II=NPLUS,IGAMAX                                             
      BB=-B*W(II)                                                       
      DO 22 JJ=II,IGAMAX                                                
      W(KK)=W(KK)-BB*W(JJ)                                              
      KK=KK+1                                                           
   22 CONTINUE                                                          
      W(KK)=BB                                                          
      KK=KK+INCINV                                                      
   21 CONTINUE                                                          
      W(KK)=B                                                           
   15 GO TO (27,24),IINV                                                
   24 I=I+1                                                             
      IF (I-N) 2,2,25                                                   
C                                                                       
C     END OF THE LOOP THAT STARTED AT STNO 2                            
C                                                                       
   25 IINV=1                                                            
      FF=0.0D0                                                             
      KL=NN                                                             
      DO 26 I=1,M                                                       
      KL=KL+1                                                           
      F(I)=W(KL)                                                        
      FF=FF+F(I)*F(I)                                                   
   26 CONTINUE                                                          
      ISS=1                                                             
C                                                                       
C     ALLOW THE USER TO PRINT STUFF FOR EACH ITERATION                  
C                                                                       
 27   IF ( IPRINT .EQ. 1 ) CALL MINMON ( ITC, MC, 0, N, M, X, FF, F )   
C                                                                       
C     ARE WE DONE                                                       
C     SEE THE INITIAL COMMENTS FOR THE STOPPING CRITERIA                
C                                                                       
      IF ( ITC .LT. 4 ) GO TO 34                                        
      IF ( CHANGE .LE. 1 ) GO TO 33                                     
      IF ( F4-FF .GT. TOL*FF .AND. FF .GT. TOL2 ) GO TO 34              
   10 CONTINUE                                                          
   33 RETURN                                                            
C                                                                       
C     SAVE FUNCTION VALUES FROM FOUR PAST ITERATIONS                    
C                                                                       
 34   F4 = F3                                                           
      F3 = F2                                                           
      F2 = F1                                                           
      F1 = FF                                                           
C                                                                       
C                                                                       
C     START OF A NEW ITERATION                                          
C                                                                       
      ITC=ITC+1                                                         
      K=N                                                               
      KK=KST                                                            
      DO 39 I=1,N                                                       
      K=K+1                                                             
      W(K)=0.0D0                                                           
      KK=KK+N                                                           
      W(I)=0.0D0                                                           
      DO 40 J=1,M                                                       
      KK=KK+1                                                           
      W(I)=W(I)+W(KK)*F(J)                                              
   40 CONTINUE                                                          
   39 CONTINUE                                                          
      DM=0.0D0                                                             
      K=KINV                                                            
      DO 41 II=1,N                                                      
      IIP=II+N                                                          
      W(IIP)=W(IIP)+W(K)*W(II)                                          
      JL=II+1                                                           
      IF (JL-N) 42,42,43                                                
   42 DO 44 JJ=JL,N                                                     
      JJP=JJ+N                                                          
      K=K+1                                                             
      W(IIP)=W(IIP)+W(K)*W(JJ)                                          
      W(JJP)=W(JJP)+W(K)*W(II)                                          
   44 CONTINUE                                                          
      K=K+1                                                             
   43 IF (DM-DABS(W(II)*W(IIP))) 45,41,41                               
   45 DM=DABS(W(II)*W(IIP))                                             
      KL=II                                                             
   41 CONTINUE                                                          
C                                                                       
C     ACCUMULATE THE CHANGE VECTOR IN  W(I)  1 =< I =< N                
C     AND THE MAXIMUM (OVER I) CHANGE/TOLERANCE                         
C                                                                       
      II=N+MPLUSN*KL                                                    
      CHANGE=0.0D0                                                         
      DO 46 I=1,N                                                       
      JL=N+I                                                            
      W(I)=0.0D0                                                           
      DO 47 J=NPLUS,NN                                                  
      JL=JL+MPLUSN                                                      
      W(I)=W(I)+W(J)*W(JL)                                              
   47 CONTINUE                                                          
      II=II+1                                                           
      W(II)=W(JL)                                                       
      W(JL)=X(I)                                                        
      DELTA = DMAX1( TOL*DABS(X(I)), TOL2 )                             
      CHANGE = DMAX1( CHANGE, DABS(W(I))/DELTA )                        
   46 CONTINUE                                                          
      DO 49 I=1,M                                                       
      II=II+1                                                           
      JL=JL+1                                                           
      W(II)=W(JL)                                                       
      W(JL)=F(I)                                                        
   49 CONTINUE                                                          
      FC=FF                                                             
      ACC=0.1D0/CHANGE                                                    
      IT=3                                                              
      XC=0.0D0                                                             
      XL=0.0D0                                                             
      IS=3                                                              
      XSTEP=-DMIN1(0.50D0,ESCALE/CHANGE)                                
C                                                                       
   51 CALL VD01A (IT,XC,FC,6,ACC,0.10D0,XSTEP)                          
      GO TO (52,53,53,53),IT                                            
C                                                                       
   52 MC=MC+1                                                           
      IF (MC-MAXFUN) 54,54,55                                           
55    IRET=1                                                            
      ISS=2                                                             
      GO TO 53                                                          
   54 XL=XC-XL                                                          
      DO 57 J=1,N                                                       
      X(J)=X(J)+XL*W(J)                                                 
   57 CONTINUE                                                          
      XL=XC                                                             
      FC = CALFUN (M,N,F,X)                                             
      GO TO (59,59,60),IS                                               
   60 K=N                                                               
      IF (FC-FF) 61,51,62                                               
   61 IS=2                                                              
      FMIN=FC                                                           
      FSEC=FF                                                           
      GO TO 63                                                          
   62 IS=1                                                              
      FMIN=FF                                                           
      FSEC=FC                                                           
      GO TO 63                                                          
   59 IF (FC-FSEC) 64,51,51                                             
   64 K=KSTORE                                                          
      GO TO (75,74),IS                                                  
   75 K=N                                                               
   74 IF (FC-FMIN) 65,51,66                                             
   66 FSEC=FC                                                           
      GO TO 63                                                          
   65 IS=3-IS                                                           
      FSEC=FMIN                                                         
      FMIN=FC                                                           
   63 DO 67 J=1,N                                                       
      K=K+1                                                             
      W(K)=X(J)                                                         
   67 CONTINUE                                                          
      DO 68 J=1,M                                                       
      K=K+1                                                             
      W(K)=F(J)                                                         
   68 CONTINUE                                                          
      GO TO 51                                                          
   53 K=KSTORE                                                          
      KK=N                                                              
      GO TO (69,70,69),IS                                               
   70 K=N                                                               
      KK=KSTORE                                                         
   69 SUM=0.0D0                                                            
      DM=0.0D0                                                             
      JJ=KSTORE                                                         
      DO 71 J=1,N                                                       
      K=K+1                                                             
      KK=KK+1                                                           
      JJ=JJ+1                                                           
      X(J)=W(K)                                                         
      W(JJ)=W(K)-W(KK)                                                  
   71 CONTINUE                                                          
      DO 72 J=1,M                                                       
      K=K+1                                                             
      KK=KK+1                                                           
      JJ=JJ+1                                                           
      F(J)=W(K)                                                         
      W(JJ)=W(K)-W(KK)                                                  
      SUM=SUM+W(JJ)*W(JJ)                                               
      DM=DM+F(J)*W(JJ)                                                  
   72 CONTINUE                                                          
      GO TO (73,10),ISS                                                 
   73 IF (ILESS.EQ.0) GO TO 83                                          
      J=KINV                                                            
      KK=NPLUS-KL                                                       
      DO 76 I=1,KL                                                      
      K=J+KL-I                                                          
      J=K+KK                                                            
      W(I)=W(K)                                                         
      W(K)=W(J-1)                                                       
   76 CONTINUE                                                          
      IF (KL-N) 77,78,78                                                
   77 KL=KL+1                                                           
      JJ=K                                                              
      DO 79 I=KL,N                                                      
      K=K+1                                                             
      J=J+NPLUS-I                                                       
      W(I)=W(K)                                                         
      W(K)=W(J-1)                                                       
   79 CONTINUE                                                          
      W(JJ)=W(K)                                                        
      B=1.0D0/W(KL-1)                                                      
      W(KL-1)=W(N)                                                      
      GO TO 88                                                          
   78 B=1.0D0/W(N)                                                         
   88 K=KINV                                                            
      DO 80 I=1,ILESS                                                   
      BB=B*W(I)                                                         
      DO 81 J=I,ILESS                                                   
      W(K)=W(K)-BB*W(J)                                                 
      K=K+1                                                             
   81 CONTINUE                                                          
      K=K+1                                                             
   80 CONTINUE                                                          
   83 IF (FMIN.LT.FF) GO TO 82                                          
      CHANGE=0.0D0                                                         
      GO TO 84                                                          
   82 FF=FMIN                                                           
      CHANGE=DABS(XC)*CHANGE                                            
   84 XL=-DM/FMIN                                                       
      SUM=1.0D0/DSQRT(DMAX1(1.0D-16,SUM+DM*XL))                            
      K=KSTORE                                                          
      DO 85 I=1,N                                                       
      K=K+1                                                             
      W(K)=SUM*W(K)                                                     
      W(I)=0.0D0                                                           
   85 CONTINUE                                                          
      DO 86 I=1,M                                                       
      K=K+1                                                             
      W(K)=SUM*(W(K)+XL*F(I))                                           
      KK=NN+I                                                           
      DO 87 J=1,N                                                       
      KK=KK+MPLUSN                                                      
      W(J)=W(J)+W(KK)*W(K)                                              
   87 CONTINUE                                                          
   86 CONTINUE                                                          
      GO TO 12                                                          
      RETURN
      END                                                               
C
      SUBROUTINE VD01A (ITEST,X,F,MAXFUN,ABSACC,RELACC,XSTEP)           
C                                                                       
C     THIS ROUTINE IS USED BY VA02A                                     
C                                                                       
C     11/6/75 - STOLLEN FROM AMDLIB                                     
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
      GO TO (1,2,2),ITEST                                               
    2 IS=6-ITEST                                                        
      ITEST=1                                                           
      IINC=1                                                            
      XINC=XSTEP+XSTEP                                                  
      MC=IS-3                                                           
      IF (MC) 4,4,15                                                    
    3 MC=MC+1                                                           
      IF (MAXFUN-MC) 12,15,15                                           
   12 ITEST=4                                                           
   43 X=DB                                                              
      F=FB                                                              
      IF (FB-FC) 15,15,44                                               
   44 X=DC                                                              
      F=FC                                                              
   15 RETURN                                                            
    1 GO TO (5,6,7,8),IS                                                
    8 IS=3                                                              
    4 DC=X                                                              
      FC=F                                                              
      X=X+XSTEP                                                         
      GO TO 3                                                           
    7 IF (FC-F) 9,10,11                                                 
   10 X=X+XINC                                                          
      XINC=XINC+XINC                                                    
      GO TO 3                                                           
    9 DB=X                                                              
      FB=F                                                              
      XINC=-XINC                                                        
      GO TO 13                                                          
   11 DB=DC                                                             
      FB=FC                                                             
      DC=X                                                              
      FC=F                                                              
   13 X=DC+DC-DB                                                        
      IS=2                                                              
      GO TO 3                                                           
    6 DA=DB                                                             
      DB=DC                                                             
      FA=FB                                                             
      FB=FC                                                             
   32 DC=X                                                              
      FC=F                                                              
      GO TO 14                                                          
    5 IF (FB-FC) 16,17,17                                               
   17 IF (F-FB) 18,32,32                                                
   18 FA=FB                                                             
      DA=DB                                                             
   19 FB=F                                                              
      DB=X                                                              
      GO TO 14                                                          
   16 IF (FA-FC) 21,21,20                                               
   20 XINC=FA                                                           
      FA=FC                                                             
      FC=XINC                                                           
      XINC=DA                                                           
      DA=DC                                                             
      DC=XINC                                                           
   21 XINC=DC                                                           
      IF ((D-DB)*(D-DC)) 32,22,22                                       
   22 IF (F-FA) 23,24,24                                                
   23 FC=FB                                                             
      DC=DB                                                             
      GO TO 19                                                          
   24 FA=F                                                              
      DA=X                                                              
   14 IF (FB-FC) 25,25,29                                               
   25 IINC=2                                                            
      XINC=DC                                                           
      IF (FB-FC) 29,45,29                                               
   29 D=(FA-FB)/(DA-DB)-(FA-FC)/(DA-DC)                                 
      IF (D*(DB-DC)) 33,37,37                                           
   37 D=0.5D0*(DB+DC-(FB-FC)/D)                                           
      IF (DABS(D-X)-DABS(ABSACC)) 34,34,35                              
   35 IF (DABS(D-X)-DABS(D*RELACC)) 34,34,36                            
   34 ITEST=2                                                           
      GO TO 43                                                          
   36 IS=1                                                              
      X=D                                                               
      IF ((DA-DC)*(DC-D)) 3,26,38                                       
   38 IS=2                                                              
      GO TO (39,40),IINC                                                
   39 IF (DABS(XINC)-DABS(DC-D)) 41,3,3                                 
   33 IS=2                                                              
      GO TO (41,42),IINC                                                
   41 X=DC                                                              
      GO TO 10                                                          
   40 IF (DABS(XINC-X)-DABS(X-DC)) 42,42,3                              
   42 X=0.5D0*(XINC+DC)                                                   
      IF ((XINC-X)*(X-DC)) 26,26,3                                      
   45 X=0.5D0*(DB+DC)                                                     
      IF ((DB-X)*(X-DC)) 26,26,3                                        
   26 ITEST=3                                                           
      GO TO 43                                                          
      RETURN
      END                                                               
C
      SUBROUTINE VMMIN(DIAG,EPS,X,XP,G,GP,S,T,GB,Q1,Q2,H,FS,FJACOB,NA,
     * N,IPRINT,NRAN,MAXFUN,ISIG)                                        
C                                                                       
C     AMD-DAVIDON VARIABLE-METRIC MINIMIZATION ADAPTED FOR LEAST-SQUARES
C                                                                       
C     THIS IS A VERSION OF THE MAIN PROGRAM "MATMIN" TO BE FOUND IN     
C     ARGONNE'S AMDLIB.  IT IS THE ORIGINAL DAVIDON VARIABLE METRIC     
C     FITTER CONVERTED TO SUBROUTINE FORM AND PROVIDED WITH THE         
C     ABILITY TO BASE THE INITIAL METRIC ON THE (JACOBIAN TRANSPOSE)*   
C     JACOBIAN APPROXIMATION TO THE SECOND DERIVATIVE MATRIX THAT IS    
C     APPROPRIATE TO SUM-OF-SQUARES FUNCTIONS.  THE STOPPING CRITERION  
C     HAS ALSO BEEN CHANGED, AND THE ABILITY TO IMPOSE LINEAR           
C     CONSTRAINTS ON THE PARAMETERS HAS BEEN REMOVED.  OTHERWISE        
C     THE LOGIC OF THE PROGRAM IS UNCHANGED FROM                        
C     THAT OF DAVIDON -- THE SUM-OF-SQUARES PROPERTY IS USED ONLY FOR   
C     THE INITIAL METRIC.                                               
C                                                                       
C     IT IS POSSIBLE TO USE THIS SUBROUTINE TO FIND THE MINIMA          
C     OF GENERAL FUNCTIONS BUT IN SUCH CASES DIAG MUST BE SET           
C     TO A NON-ZERO NUMBER.                                             
C                                                                       
C     PARAMETER:                                                        
C                                                                       
C     DIAG - CONTROLS THE INITIAL CHOICE OF METRIC:                     
C        0 - METRIC = ( 2 JTJ )**(-1)  WHERE J IS THE JACOBIAN.         
C       >0 - METRIC =  DIAG * UNIT MATRIX.                              
C       <0 - METRIC =  \DIAG\ * DELTA(I,J) * X(I)**2; I.E. A SCALED     
C                      DIAGONAL MATRIX.                                 
C                                                                       
C     EPS - REQUIRED RELATIVE ACCURACY.  THE FITTER WILL STOP IF EITHER 
C        OF THE FOLLOWING TWO CONDITIONS OBTAINS:                       
C        A)  \ X(I-1) - X(I) \ < EPS*\X(I)\  FOR ALL PARAMETERS         
C        OR                                                             
C        B)  FUNC(I-2) - FUNC(I) < EPS*FUNC(I)                          
C        HERE I LABELS THE ITERATION.                                   
C                                                                       
C     X - INITIAL AND FINAL PARAMETER VECTOR.  THE FIT STARTS AT THE    
C        INPUT VALUE OF X AND X CONTAINS THE LOCATION OF THE MINIMUM    
C        UPON RETURN.                                                   
C                                                                       
C     XP, G, GP, S, T, GB, Q1, Q2 - WORK ARRAYS OF DIMENSION N.         
C                                                                       
C     H - AN N X N WORK ARRAY THAT WILL CONTAIN THE VARIANCE MATRIX     
C            ( .5  D**2 F / DX DX ) ** (-1)                             
C         UPON RETURN.                                                  
C                                                                       
C     FS - WORK ARRAY OF DIMENSION NA                                   
C     FJACOB - WORK ARRAY OF DIMENSION  NA X N                          
C                                                                       
C     NA - NUMBER OF TERMS IN THE SUM-OF-SQUARES                        
C     N - NUMBER OF PARAMETERS THAT ARE TO BE VARIED IN THE FIT.        
C                                                                       
C     IPRINT - CONTROLS PRINTING DURING THE FIT:                        
C        0 - ONLY PRINT ERROR MESSAGES                                  
C        1 - CALL MINMON AT EACH ITERATION TO ALLOW IT TO PRINT         
C                                                                       
C     NRAN - NUMBER OF RANDOM STEPS TO USE TO CONFIRM A MINIMUM         
C                                                                       
C     MAXFUN - MAXIMUM ALLOWED NUMBER OF FUNCTION REFERENCES.           
C                                                                       
C     ISIG - A RETURN CODE IS PLACED HERE:                              
C        0 - ALL WENT O.K.                                              
C        1 - EXCEEDED MAXFUN, FIT HAS STOPPED WITH BEST VALUE FOR THE   
C            MINIMUM IN X.                                              
C                                                                       
C                                                                       
C     SUBROUTINES THAT MUST BE SUPPLIED BY THE USER:                    
C                                                                       
C                                                                       
C     F = CALGRA ( NA, N, FS, X, FJACOB, G )                            
C                                                                       
C     THIS IS CALLED EACH TIME A FUNCTION AND ITS GRADIENT IS DESIRED.  
C                                                                       
C     NA AND N ARE AS ABOVE.                                            
C     FS IS AS ABOVE AND SHOULD BE DEFINED TO HAVE THE INDIVIDUAL       
C        TERMS IN THE SUM-OF-SQUARES.                                   
C     F (THE FUNCTIONAL VALUE) IS THE SUM OF THE SQUARES OF THE FS.     
C     X IS AS ABOVE AND IS THE LOCATION AT WHICH F, FS, FJACOB AND G    
C       ARE TO BE EVALUATED.                                            
C     FJACOB IS AS ABOVE AND MUST BE SET TO THE JACOBIAN:               
C          FJACOB(I,J) = D(FS(I))/D(X(J))                               
C     G IS AS ABOVE AND MUST BE SET TO THE GRADIENT OF F:               
C          G(J) = SUM(I) 2*FS(I)*FJACOB(I,J)                            
C                                                                       
C     NOTE:  IT IS NECESSARY TO COMPUTE FS AND FJACOB ONLY ON THE       
C            FIRST FUNCTION CALL.  IF DIAG IS NON-ZERO IT IS NEVER      
C            NECESSARY TO COMPUTE FS AND FJACOB.  IN THE LATTER CASE    
C            IT IS NOT NECESSARY THAT THE FUNCTION BE A SUM-OF-SQUARES  
C            AND THE ARRAYS FS AND FJACOB AND THE PARAMETER NA WILL     
C            NOT BE USED.                                               
C                                                                       
C                                                                       
C     CALL MINMON ( IT, NCALL, NCALL, N, NA, X, F, FJACOB )             
C                                                                       
C     IF IPRINT = 1, THIS CALL IS MADE AT EACH ITERATION TO ALLOW       
C     PRINTING TO SHOW THE PROGRESS OF THE FIT.  THE PARAMETERS ARE     
C     AS ABOVE; THE FIRST THREE PARAMETERS ARE:                         
C                                                                       
C     IT - THE ITERATION COUNTER - INCREMENTED BY 1 FOR EACH ITERATION. 
C     NCALL - THE NUMBER OF FUNCITON REFERENCES THAT HAVE BEEN MADE.    
C                                                                       
C                                                                       
C     1/23/76 - CONVERTED FROM MATMIN - M. H. MACFARLANE                
C     2/19/76 - ADD COMMENTS, DIAG < 0 FEATURE - S. C. PIEPER           
C                                                                       
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,T,GB,Q1,Q2,H,FS,FJACOB                  
C                                                                       
      COMMON/VMM/F,GS,EL,SL,FP,GSP,T0,Z,Q,A,GSS,F0,GTP,FB,GTT,GSB,F1,F2
     *,F3,FFIRST,FAC,DELTA,E,P,NCALL,K,NC,NSSW1,NSSW2,MS,IT,M1,L        
      COMMON/VMM2/IY                                                    
C                                                                       
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),T(N),GB(N),Q1(N),Q2(N),
     * H(N,N),FS(NA),FJACOB(NA,N)                                             
C                                                                       
      ISIG=0                                                            
      FAC= DABS(DIAG)                                                   
      E=EPS                                                             
      NSSW1=IPRINT                                                      
      NSSW2=0                                                           
      K=IABS(NRAN)                                                      
      NC=0                                                              
      DELTA=1                                                           
      NCALL=0                                                           
      IY=1                                                              
C                                                                       
      IF(FAC.EQ.0.0D0) GO TO 1070                                           
C                                                                       
C     HERE WE USE A DIAGONAL MATRIX AS THE INITIAL METRIC               
C                                                                       
      DO 35 II=1,N                                                      
      DO 35 JJ=1,N                                                      
      H(II,JJ)=0.0D0                                                        
      IF (II-JJ) 35,30,35                                               
   30 H(II,JJ)= FAC                                                     
      IF ( DIAG .LT.0.0D0 ) H(II,JJ) = FAC*X(II)**2                        
   35 CONTINUE                                                          
 1070 M1=4                                                              
      GO TO 840                                                         
   65 M1=1                                                              
      IF(FAC.NE.0.0D0) GO TO 810                                            
   70 MS=0                                                              
C                                                                       
  100 IT=0                                                              
      F=CALGRA(NA,N,FS,X,FJACOB,G)                                      
      NCALL=NCALL+1                                                     
      F1=F                                                              
      F2=0.0D0                                                              
      F3=0.0D0                                                              
      FFIRST=F                                                          
      IF(NSSW1.EQ.1) CALL MINMON(IT,NCALL,NCALL,N,NA,X,F,FJACOB)        
      IF(FAC.NE.0.0D0) GO TO 840                                            
      DO 1050 II=1,N                                                    
      DO 1050 JJ=II,N                                                   
      H(II,JJ)=0.0D0                                                        
      DO 1040 IA=1,NA                                                   
 1040 H(II,JJ)=H(II,JJ)+FJACOB(IA,II)*FJACOB(IA,JJ)                     
      H(II,JJ)=2.0D0*H(II,JJ)                                               
 1050 H(JJ,II)=H(II,JJ)                                                 
      CALL SYMINV(H,N,Q1,0,DET,N)                                       
      DO 1060 II=2,N                                                    
      JTOP=II-1                                                         
      DO 1060 JJ=1,JTOP                                                 
 1060 H(II,JJ)=H(JJ,II)                                                 
      GO TO 840                                                         
C                                                                       
  120 CONTINUE                                                          
      M1=2                                                              
      IF(NCALL.LE.MAXFUN) GO TO 250                                     
      ISIG=1                                                            
      GO TO 800                                                         
  250 CALL READY(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                     
      IF(L-2)570,300,500                                                
  300 CALL AIM(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                       
      IF(L.NE.1)GO TO 500                                               
      CALL FIRE(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                      
      IF(L.GT.2)GO TO 300                                               
  500 CALL DRESS(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                     
      GO TO 120                                                         
  570 IF (NSSW1) 600,580,600                                            
  580 M1=5                                                              
      GO TO 840                                                         
  590 M1=2                                                              
  600 CALL STUFF(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                     
      IF(L.EQ.1)GO TO 100                                               
C                                                                       
 800  F = CALGRA( NA, N, FS, X, FJACOB, G )                             
      NCALL=NCALL+1                                                     
      M1=3                                                              
 810  IF ( M1 .EQ. 1 ) GO TO 70                                         
  840 GO TO (120, 120, 880, 65, 590), M1                                
  880 CONTINUE                                                          
C                                                                       
      DO 1200 I=1,N                                                     
      DO 1200 J=1,N                                                     
 1200 H(I,J)=2*H(I,J)                                                   
      RETURN                                                            
C                                                                       
      END                                                               
C                                                                       
      SUBROUTINE STUFF(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)               
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,T,GB,FS,H,FJACOB                        
      COMMON/VMM/F,GS,EL,SL,FP,GSP,T0,Z,Q,A,GSS,F0,GTP,FB,GTT,GSB,F1,F2
     *,F3,FFIRST,FAC,DELTA,E,P,NCALL,K,NC,NSSW1,NSSW2,MS,IT,M1,L        
      COMMON/VMM2/IY                                                    
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),T(N),GB(N), H(N,N),FS(NA),
     * FJACOB(NA,N)                                                         
C                                                                       
      L=2                                                               
      K=K-1                                                             
      IF (K) 640,610,610                                                
C                                                                       
  610 L=1                                                               
      PSQ=2.0D0*DSQRT(FP)                                                   
      P=DSQRT(PSQ)                                                      
      MS=MS+1                                                           
      WRITE(6,1)MS,DELTA,GS                                               
      DO 620 I=1,N                                                      
      CALL RRANDU (IY,YFL)
  620 T(I)=YFL-0.5D0                                                       
      CALL MATMPY(N,N,H,T,S)                                            
      CALL MATMPY(1,N,S,T,EL)                                           
      EL=P/DSQRT(EL)                                                    
      DO 630 I=1,N                                                      
  630 X(I)=X(I)+EL*S(I)                                                 
C                                                                       
  640 RETURN                                                            
  1   FORMAT ( ' RANDOM STEP', I4,'  DELTA =', D14.5,'  GS =', D14.5 )                                                                 
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE RRANDU(IX,YFL)
      IMPLICIT REAL*8 (A-H,O-Z)                                         
      REAL*8 YFL                                                        
      IY=IX*65539                                                       
      IF (IY)5,6,6                                                      
    5 IY=IY+2147483647+1                                                
    6 YFL=IY                                                            
      YFL=YFL*0.4656613D-9                                               
      IX=IY                                                             
      RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE READY(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)               
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,T,GB,FS,H,FJACOB                        
      COMMON/VMM/F,GS,EL,SL,FP,GSP,T0,Z,Q,A,GSS,F0,GTP,FB,GTT,GSB,F1,F2
     *,F3,FFIRST,FAC,DELTA,E,P,NCALL,K,NC,NSSW1,NSSW2,MS,IT,M1,L        
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),T(N),GB(N), H(N,N),FS(NA),
     * FJACOB(NA,N)                                                         
C                                                                       
      L=1                                                               
      CALL MATMPY(N,N,H,G,S)                                            
      DO 205 I=1,N                                                      
  205 S(I)=-S(I)                                                        
      CALL MATMPY(1,N,S,G,GS)                                           
      DO 250 I=1,N                                                      
      IF(DABS(S(I)).GT.E*DABS(X(I))) GO TO 210                          
  250 CONTINUE                                                          
      GO TO 240                                                         
C                                                                       
  210 L=2                                                               
      EL=2.0D0                                                            
      T0=EL*F/GS                                                        
      IF (T0+EL) 213, 213, 212                                          
  212 EL=-T0                                                            
  213 SL=-GS                                                            
      DO 215 I=1,N                                                      
  215 XP(I)=X(I)+EL*S(I)                                                
      FP=CALGRA(NA,N,FS,XP,FJACOB,GP)                                   
      NCALL=NCALL+1                                                     
      IF(FP.GT.100.0D0*FFIRST) FP=10.0D0*FFIRST                                 
      F3=F2                                                             
      F2=F1                                                             
      F1=FP                                                             
      IF((DABS(F1-F3).GT.E*DABS(F1)).OR.(FP.EQ.10.0D0*FFIRST))GO TO 260    
      L=1                                                               
      GO TO 240                                                         
  260 CALL MATMPY(1,N,S,GP,GSP)                                         
      IF (-GSP) 240,240,220                                             
  220 IF (F-FP) 240,240,225                                             
C                                                                       
  225 L=3                                                               
      FB=FP                                                             
      DO 230 I=1,N                                                      
      GB(I)=GP(I)                                                       
      T(I)=XP(I)                                                        
  230 CONTINUE                                                          
      IF (EL-2.0D0) 240,235,240                                           
C                                                                       
  235 L=4                                                               
      DELTA=DELTA+DELTA                                                 
      T0=1.0D0/SL                                                         
C                                                                       
  240 RETURN                                                            
      END                                                               
C                                                                       
C                                                                       
      SUBROUTINE AIM(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                 
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,T,GB,FS,H,FJACOB                        
      COMMON/VMM/F,GS,EL,SL,FP,GSP,T0,Z,Q,A,GSS,F0,GTP,FB,GTT,GSB,F1,F2
     *,F3,FFIRST,FAC,DELTA,E,P,NCALL,K,NC,NSSW1,NSSW2,MS,IT,M1,L        
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),T(N),GB(N), H(N,N),FS(NA),
     * FJACOB(NA,N)                                                         
C                                                                       
      L=1                                                               
      Z=3.0D0/EL*(F-FP)+GS+GSP                                            
      Q=DABS(Z*DSQRT(1.0D0-(GS/Z)*(GSP/Z)))                               
      A=(Q-Z+GSP)/(Q+Q-GS+GSP)                                          
      T0=EL/3.0D0*(Q+Q+Z+GSP)*A*A                                         
      F0=FP-T0                                                          
      CALL MATMPY(N,N,H,GP,T)                                           
      DO 305 I=1,N                                                      
  305 T(I)=(GSP/SL)*S(I)-T(I)                                           
      CALL MATMPY(1,N,T,GP,GTP)                                         
      IF (T0+T0+GTP) 315,310,310                                        
  310 DO 312 I=1,N                                                      
  312 T(I)=XP(I)+A*(X(I)-XP(I))                                         
      GO TO 340                                                         
  315 IF (F+F+GTP) 310,320,320                                          
  320 DO 322 I=1,N                                                      
  322 T(I)=T(I)+XP(I)                                                   
      FB=CALGRA(NA,N,FS,T,FJACOB,GB)                                    
      NCALL=NCALL+1                                                     
      IF(FB.GT.100.0D0*FFIRST) FB=10.0D0*FFIRST                                 
      IF (F0-FB) 310,325,325                                            
C                                                                       
  325 L=3                                                               
      DO 327 I=1,N                                                      
  327 S(I)=T(I)-XP(I)                                                   
      CALL MATMPY(1,N,S,GB,GTT)                                         
      GTT=GTT-GTP                                                       
      IF (GTT) 340,330,330                                              
C                                                                       
  330 L=2                                                               
      GSS=GTT                                                           
      SL=-GTP                                                           
      EL=1.0D0                                                            
C                                                                       
  340 RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE FIRE(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)                
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,T,GB,FS,H,FJACOB                        
      COMMON/VMM/F,GS,EL,SL,FP,GSP,T0,Z,Q,A,GSS,F0,GTP,FB,GTT,GSB,F1,F2
     *,F3,FFIRST,FAC,DELTA,E,P,NCALL,K,NC,NSSW1,NSSW2,MS,IT,M1,L        
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),T(N),GB(N), H(N,N),FS(NA),
     * FJACOB(NA,N)                                                         
      EQUIVALENCE (TEMP,GTT)                                            
C                                                                       
      L=1                                                               
      TEMP=A/(1.0D0-A)                                                    
      FB=CALGRA(NA,N,FS,T,FJACOB,GB)                                    
      NCALL=NCALL+1                                                     
      IF(FB.GT.100.0D0*FFIRST) FB=10.0D0*FFIRST                                 
      CALL MATMPY(1,N,S,GB,GSB)                                         
      T0=F                                                              
      IF (T0-FP) 403, 403, 402                                          
  402 T0=FP                                                             
  403 IF (T0-FB+E*DABS(FB)) 415, 405, 405                               
  405 GSS=Q+Q                                                           
      T0=GSB*(TEMP-1.0D0/TEMP)                                            
      IF (DABS(T0)-Q) 430,410,410                                       
C                                                                       
  410 L=2                                                               
      GO TO 440                                                         
C                                                                       
  415 L=3                                                               
      IF (FP-F) 425,420,420                                             
C                                                                       
 420  EL=(1.0D0-A)*EL                                                     
      FP=FB                                                             
      GSP=GSB                                                           
      DO 422 I=1,N                                                      
      XP(I)=T(I)                                                        
      GP(I)=GB(I)                                                       
  422 CONTINUE                                                          
      GO TO 440                                                         
C                                                                       
 425  EL=EL*A                                                           
      F=FB                                                              
      GS=GSB                                                            
      DO 427 I=1,N                                                      
      X(I)=T(I)                                                         
      G(I)=GB(I)                                                        
  427 CONTINUE                                                          
      GO TO 440                                                         
C                                                                       
  430 GSS=GSS+T0                                                        
      DO 435 I=1,N                                                      
  435 G(I)=(GB(I)-G(I))*TEMP+(GP(I)-GB(I))/TEMP                         
C                                                                       
  440 RETURN                                                            
      END                                                               
C                                                                       
      SUBROUTINE DRESS(NA,N,X,XP,G,GP,S,T,GB,FS,H,FJACOB)               
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL X,XP,G,GP,S,T,GB,H,FS,FJACOB                        
      COMMON/VMM/F,GS,EL,SL,FP,GSP,T0,Z,Q,A,GSS,F0,GTP,FB,GTT,GSB,F1,F2
     *,F3,FFIRST,FAC,DELTA,E,P,NCALL,K,NC,NSSW1,NSSW2,MS,IT,M1,L        
      DIMENSION X(N),XP(N),G(N),GP(N),S(N),T(N),GB(N), H(N,N),FS(NA),
     * FJACOB(NA,N)                                                         
C                                                                       
      GO TO (505,520,530,525), L                                        
C                                                                       
  505 CALL MATMPY(N,N,H,G,X)                                            
      CALL MATMPY(1,N,X,G,T0)                                           
      IF (T0-GSS**2/SL-E*DABS(T0)) 515,510,510                          
  510 DO 512 II=1,N                                                     
      DO 512 JJ=1,N                                                     
  512 H(II,JJ)=H(II,JJ)-X(II)*X(JJ)/T0                                  
      DELTA=DELTA*(EL*GSS/T0)                                           
      T0=EL/GSS                                                         
      GO TO 525                                                         
C                                                                       
  515 CONTINUE                                                          
C                                                                       
  520 DELTA=DELTA*(EL*SL/GSS)                                           
      T0=EL/GSS-1.0D0/SL                                                  
C                                                                       
  525 DO 527 II=1,N                                                     
      DO 527 JJ=1,N                                                     
  527 H(II,JJ)=H(II,JJ)+T0*S(II)*S(JJ)                                  
C                                                                       
  530 IT=IT+1                                                           
      F=FB                                                              
      DO 532 I=1,N                                                      
      G(I)=GB(I)                                                        
      X(I)=T(I)                                                         
  532 CONTINUE                                                          
      IF(NSSW1.EQ.1) CALL MINMON(IT,NCALL,NCALL,N,NA,X,F,FJACOB)        
      RETURN                                                            
C                                                                       
      END                                                               
C                                                                       
      SUBROUTINE MATMPY(M,N,H,G,S)                                      
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL H,G,S                                               
      REAL*8 H,G,S                                                      
      DIMENSION H(N,N),G(N),S(N)                                        
C                                                                       
      DO 720 II=1,M                                                     
      S(II)=0.0D0                                                         
      DO 720 JJ=1,N                                                     
  720 S(II)=S(II)+H(JJ,II)*G(JJ)                                        
C                                                                       
      RETURN                                                            
      END                                                               

      FUNCTION CALFUN ( NPTST, MPARA, FS, XS )                          
CDC      FUNCTION  CALFN2 ( NPTST, MPARA, FS, XS, FGRAD, TGRAD )        
C                                                                       
C  CALCULATE CHI SQUARE SUM WITH GIVEN OPTICAL PARAMETER SET            
C                                                                       
C     ENTRIES:                                                          
C       CALFUN - COMPUTES THE CHI-SQUARED SUM AND THE ARRAY FS          
C                OF INDIVIDUAL ERRORS                                   
C         CHI-SUM  =  SUM(I) FS(I)**2                                   
C                                                                       
C       CALGRA - COMPUTES THE ABOVE AND THE JACOBIAN ARRAY FGRAD        
C           FGRAD(I,M) = D( FS(I) ) / D( PARAMETER(M) )                 
C                 ALSO FINDS THE TOTAL GRADIEND TGRAD                   
C           TGRAD(M) = D( CHI-SUM ) / D( PARAMETER(M) )                 
C                                                                       
C     IF  MOD ( IPRINT, 100 )  IS NON-ZERO, WE PRINT A LINE EACH        
C     TIME CALFUN OR CALGRA IS CALLED.                                  
C                                                                       
C  THIS ROUTINE IS CALLED FROM VA02A                                    
C  AND FROM FITEL WHEN ROCORD SPECIFIED                                 
C     AND VIA INTERFACES FOR THE OTHER FITTERS.                         
C                                                                       
C *************************************************************         
C *   THIS ROUTINE WILL NEED ADDITIONAL WORK FOR S.O. FORCES  *         
C *************************************************************         
C                                                                       
C  THIS FIRST VERSION 8/15/75  D.G.                                     
C     10/2/75 - MERGE INTO SEPTEMBER PTOLEMY - S.P.                     
C     10/7/75 - USE MAKPOT FOR POTENTIALS - S.P.                        
C     10/11/75 - TABLE OF SIN/COS(SIGMA)                                
C     10/15/75 - ELAB-DEPENDENT LMIN/LMAX AND STEPSIZE                  
C     10/21/75 - ESQ PARAMETERS                                         
C     10/23/75 - PRINT POTENTIALS; PRINT EQUALS                         
C     10/28/75 - COMPUTE GRADIENTS                                      
C     11/4/75 - OUT-OF-BOUNDS STEPS; SEPARATE PRTCHI.                   
C     11/6/75 - IMPROVED TIMING; OTHER SMALL FIXES                      
C     11/12/75 - FIX RENORMALIZATION AND WEIGHT - S.P.                  
C     11/20/75 - LAB-ANGLE SHIFT - S. P.                                
C     12/28/75 - PRINT TRANSFORMED DATA                                 
C     1/15/76 - MADE INTO A FUNCTION; USE IPRINT                        
C     1/31/76 - SMALL PRINTOUT FIXES                                    
C     7/31/76 - NEW DATA STRUCTURE; LAB ANGLES; CDC/IBM                 
C     4/15/77 - MULTIPLE CHANNELS; DATA IN STANDARD FORM; CHANGE PRINTIN
C     4/21/77 - RENORM AND ANGLE SHIFT CAN BE FIT                       
C     4/26/77 - IDENTICAL PARTICLES; DATA HAS LENGTH INDICATORS         
C     11/22/77 - NEW SETFIT, PRTCHI - S.P.                              
C     12/1/77 - RESTORE LINKULE STUFF - S.P.                            
C     2/14/78 - NUMERICAL DERIVATIVES                                   
C     4/15/78 - CALGR2 FOR CDC OVERLAY; IESIZE INSTEAD OF 9 - S.P.      
C     4/17/78 - PRINT DIAGNOSTIC MESSAGE AT ENTRY - S.P.                
C     4/24/78  - REVERSE LOMOST FOR LIMOST - S.P.                       
C         5/22/78 - SEPARATE OUT CALANG - S.P.                          
C     5/4/79 - RATMAS FROM DATA - S.P.                                  
C     12/26/79 - STANDARD MORTRAN - RPG                                 
C                                                                       
      IMPLICIT REAL*8 ( A-H, O-Z )                                      
CDC      ALLOCLEVEL FS,XS,FGRAD,TGRAD                                   
C                                                                       
      DIMENSION FS(NPTST),XS(MPARA)                                     
      DIMENSION FGRAD(NPTST, MPARA), TGRAD(MPARA)                       
C                                                                       
      INTEGER*4 Z                                                       
      COMMON /LOCKT/ Z(100)                                               
      INTEGER FACFR4, FACFR2, FACFR1                                    
      COMMON /ALLOCS/ FACFR4, FACFR2, FACFR1                            
      DIMENSION ALLOC(1)                                                
      LOGICAL*1 MLLOC(1)                                                
      REAL*4 ALLOC4(1)                                                  
      INTEGER ILLOC(1)                                                  
      INTEGER*2 ILLOC2(1)                                               
      EQUIVALENCE ( FACFR4,ALLOC(1),ILLOC(1),ILLOC2(1),ALLOC4(1),
     * MLLOC(1) )                                                          
C                                                                       
      COMMON /FLOAT/ A, ACCURA, AI, ANGSTP, ANGMIN, ANGMAX, ASO, ASOI, A
     *SYMPT, DELTVK, DWCUT, E, ECM, ELAB, EXS(5), GAMDIF, GAMSUM, TAU, T
     *AUI, AM, AMA, AMB, AMBIGA, AMBIGB, AMX, AMXCS(5), AMXCGS(5), AMP, 
     *AMT, Q, R, R0, RC, RC0, RI, RI0, RSO, RSO0, RSOI, RSOI0, SPAM, SPA
     *MP, SPAMT, STEPSZ, STEP1R, STEP1I, SUMMAX, SUMMID, SUMMIN, V, VI, 
     *VSO, VSOI, FITMUL, WRITES, ASI, RSI, RSI0, VSI, AMDMLT, STPSPR, R0
     *E, RI0E, AE, AIE, VE, VIE, FITACC, ALMNMT, SPFACP, SPFACT, R0ESQ, 
     *RI0ESQ, AESQ, AIESQ, VESQ, VIESQ, DERIVS, FITRAT, RCP, RCT, RC0P, 
     *RC0T, SUMPTS, ACCINE, COULML, ALMXMT, EXSPT(2), POWRL, POWIM, BNDA
     *SY, SCTASY, VTR,VTRI,VTL,VTLI,VTP,VTPI, PARAM1, PARAM2, PARAM3, PA
     *RAM4, PARAM5, PAR620(15), RTR, RTRI, RTL, RTLI, RTP, RTPI, RTR0, R
     *TRI0, RTL0, RTLI0, RTP0, RTPI0, ATR, ATRI, ATL, ATLI, ATP, ATPI, A
     *MXGPT(2), PHIMID, GAMPHI                                          
      DIMENSION PARAMS(20)                                              
      EQUIVALENCE ( PARAMS(1), PARAM1 )                                 
      DIMENSION AMS(5)                                                  
      EQUIVALENCE ( AMS(1), AMA )                                       
      DIMENSION VTEN(6), RTEN(6), RTEN0(6), ATEN(6)                     
      EQUIVALENCE ( VTEN(1), VTR ), (RTEN(1), RTR ), (RTEN0(1), RTR0 ), 
     *(ATEN(1), ATR )                                                   
      INTEGER PARITS, PARIPT, PARITY                                    
      COMMON /INTGER/ L, LMAX, LMIN, LOOKST, LX, MAPSUM, MAPDIF, MAXITR,
     * NAITKN, NBACK, NCOSIN, NODES, NPDIF, NPHIAD, NPPHI, NPSUM, IPRINT
     *, ISIZE, IZS(5), IZP, IZT, LSTEP, LBACK, MAXLEX, LMAXAD, MAXFUN, I
     *REINI, MODEFT, NUMRAN, NVPOLY, NPCOUL, MXCOUL, LMINSB, PARITS(5), 
     *PARIPT(2), PARITY, IPARM1, IPARM2, IPARM3, IPARM4, IPARM5, NBLKSZ,
     * MAPPHI, NFIROF, LEBACK                                           
      COMMON /CNSTNT/ PI, RT4PI, PIINV, RADIAN, DEGREE, HBARC, AMUMEV, A
     *FINE, BIGEST, SMLEST, BIGNUM, SMLNUM, BIGLOG                      
      LOGICAL*1 WASSET                                                  
      COMMON /INTRNL/ UNDEF, NOTDEF, ICHANB, ICHANW, ICNTT, EBNDS(2), ID
     *ONE, ISTRIP, IHSAVE, NEXTSV, IPOT, IWFN, LIMOST, IEXCIT, R0MASS, W
     *ASSET(40), RATMAS, LSPECS(4), NODESP(4), LSPCPT(2), NODEPT(2),ICD3
     *I                                                                 
      COMMON /LOCFIT/ IDATA, LPL, LSIGMA, LSMAT, IESETS, IMASK, IEQUAL, 
     *NPARA, IPARA, ITEMP, LSCSIG, LEPARM, LFS, NPTS, LEQNAM, MAXEQL, MA
     *XVS, LVDERV, LSDELS, LDVTYP, LDSDP, LDFDX, LFGRAD, LTGRAD, LMASK, 
     *LEQUAL, MNDRIV, MXDRIV, IESIZE, LENPRE, LENKEY, LENHED, ILNKBL    
      LOGICAL PWAVSW, STANSW, SOSWS, PWBGSW, TCSWS, SOSWT               
      COMMON /WAVCOM/ HS(2), RSTEPS(2), NSTPSS(2), NFS(2), NGS(2), NF1S(
     *2), NG1S(2), STEPR, STEPI, IWAVR, IWAVI, NUMAIT, NUMFIT, IOFFIT, P
     *WAVSW, NBAKCM, STANSW, SOSWS(2), NUMJS(2), RIS(2), AIS(2), RSOS(2)
     *, ASOS(2), RSOIS(2), ASOIS(2), IRLVS(2), IIMVS(2), IZPS(2), IZTS(2
     *), ASYMPS(2), V0RS(2), V0IS(2), V0SORS(2), V0SOIS(2), ES(2), JSPS(
     *2), RSIS(2), ASIS(2), V0SIS(2), ISORS(2), ISOIS(2), ALPHAR, ALPHAI
     *, NFIRST, PWBGSW, LASTL(2), LASTNF(2), LSKIPS(2), ISTATS(2), POWRL
     *S(2), POWIMS(2), NSTP2S(2), XFACS(3,2), IVWORK, ICENTR(2), TCSWS(2
     *), NWAVEF, VTENS(2,6), RTENS(2,6), ATENS(2,6), ITENS(2,6), IWAVRS(
     *2), IWAVIS(2), LASTZR, SOSWT, IINDXE(2), ITOCE(2)                 
      COMMON /KANDM/ AKI, AKO, REDMI, REDMO, RIMAX, ROMAX, ISIGS(2), RSC
     *TS(2), ASCTS(2), RCSCTS(2) ,ISMATS(2), LOMOST, LCRIT, ISCTMN, ISCT
     *CR, ETAS(2), RCSCTP(2), RCSCTT(2), LCRITS(2), RATMSS(2)           
      DIMENSION AKS(2), REDMS(2), RMAXS(2)                              
      EQUIVALENCE (AKI,AKS(1)),(REDMI, REDMS(1) ),(RIMAX,RMAXS(1)) 
      COMMON /LNKBLK/ NUMLNK, IUNIQU, LNKADR(6,12)                      
      REAL*8 LNKNAM(3,12)                                               
      EQUIVALENCE ( LNKNAM(1,1), LNKADR(1,1) )                          
C                                                                       
      DIMENSION DUMMM(1)
C                                                                       
      LOGICAL GRADSW, BACKSW, PRNTSW                                    
      CHARACTER*1 CHAR, BLANK, STAR                                     
      DATA STAR/'*'/                                                     
      DATA BLANK/' '/                                                    
C                                                                       
      GRADSW = .FALSE.                                                  
      NUMCAL = NUMCAL + 1                                               
      GO TO 1                                                           
C                                                                       
C     FOLLOWING ENTRY IS FOR GRADIANT COMPUTATIONS                      
C                                                                       
      ENTRY CALGRA ( NPTST, MPARA, FS, XS, FGRAD, TGRAD )               
CDC      ENTRY CALGR2                                                   
C
      GRADSW = .TRUE.                                                   
      NUMGRA = NUMGRA + 1                                               
C                                                                       
   1  IRET=0                                                            
      CHISQ=0.0D0                                                           
C                                                                       
      PRNTSW = MOD(IPRINT/100, 10) .GT. 0                       
      IF ( PRNTSW ) WRITE(6,3) XS                                         
 3    FORMAT (' CALFUN INPUT:', ( T16, 8G14.5 ) )                       
C                                                                       
      IF ( .NOT. GRADSW ) GO TO 10                                      
      DO 5 I = 1, MPARA                                                 
 5    TGRAD(I) = 0.0D0                                                      
 10   LEKEY = IDATA + LENPRE-1                                          
      ICNTPT=0                                                          
      LWAVR = Z(IWAVR)                                                  
      LWAVI = Z(IWAVI)                                                  
      LVREAL = Z(IRLVS(1))                                              
      LVIMAG = Z(IIMVS(1))                                              
      LCENTR = Z(ICENTR(1))                                             
C                                                                       
C     LOOP OVER THE DATA SETS WE ARE FITTING TO ( MAY BE ONLY ONE ).    
C                                                                       
      DO 789 IOFFIT = 1, IESETS                                         
C                                                                       
      LEHED = LEKEY + LENKEY+1                                          
      LLNKBL = FACFR4*Z(ILNKBL)-FACFR4 + 3*NUMLNK*(IOFFIT-1) - 3        
      LADATA = LEHED + LENHED                                           
      LLEPRM = LEPARM + IESIZE*(IOFFIT-1)                               
C                                                                       
C     GET ALL THE CHANNEL PARAMETERS                                    
C                                                                       
      ELAB = ALLOC4(LEHED+4)                                            
      IZP = ILLOC(LEHED+5)                                              
      IZT = ILLOC(LEHED+6)                                              
      AMP = ALLOC4(LEHED+7)                                             
      AMT = ALLOC4(LEHED+8)                                             
      LSKIP = ILLOC(LEHED+22)                                           
      RATMAS = ALLOC4(LEHED+24)                                         
      LLMIN = ILLOC(FACFR4*LLEPRM+1)                                    
      LLMAX = ILLOC(FACFR4*LLEPRM+2)                                    
      NSTPSS(1) = ILLOC(FACFR4*LLEPRM+3)                                
CDC      LLEPRM = LLEPRM + 2                                            
      RSTEPS(1) = ALLOC(LLEPRM+3)                                       
      ETAS(1) = ALLOC(LLEPRM+4)                                         
      AKS(1) = ALLOC(LLEPRM+5)         
      E = ALLOC(LLEPRM+6)                                               
      R0MASS = ALLOC(LLEPRM+7)                                          
      RCP = ALLOC(LLEPRM+8)                                             
      RCT = ALLOC(LLEPRM+9)                                             
      ES(1) = E                                                         
      HS(1) = AKS(1)*RSTEPS(1)                                          
C                                                                       
      DO 69 I = 1, NUMLNK                                               
      DO 69 II = 1, 3                                                   
      LNKADR(3+II,I) = ILLOC(LLNKBL+3*I+II)                             
 69   CONTINUE                                                          
C                                                                       
C     GENERATE THE POTENTIAL ARRAYS FOR WAVELJ                          
C                                                                       
      CALL SETFIT ( XS, IRET, 1 )                                       
C                                                                       
C     DID SETFIT DETECT AN OUT-OF-BOUNDS STEP                           
C                                                                       
      IF ( IRET .LE. 0 ) GO TO 900                                      
C                                                                       
C     GET DATA DESCRIPTORS AND OVERALL WEIGHT, RENORMALIZATION AND      
C     AND ANGLE SHIFT                                                   
C     RENORM - RENORMALIZATION FACTOR FOR DATA                          
C               ( DATA VALUE = RENORM*INPUT )                           
C     ASHIFT - LAB ANGLE SHIFT                                          
C     POLARIZATION DATA 0  ( NOT ACTUALLY USED YET )                    
C                                                                       
CCCC      IPOL = ILLOC(LEKEY+3)                                         
C                                                                       
      LENDAT = ILLOC(LEKEY+5)                                           
C                                                                       
      NPTSE=ILLOC(LEHED)                                                
      WEIGHT =ALLOC4(LEHED+1)                                           
      WT = DSQRT(WEIGHT)                                                
      RENORM=ALLOC4(LEHED+2)                                            
      ASHIFT = ALLOC4(LEHED+3)                                          
C                                                                       
      CHAR = BLANK                                                      
      IF ( ASHIFT.NE.0.0D0.OR.RENORM.NE.1.0D0) CHAR = STAR               
C                                                                       
      LSIG = LSIGMA + (IOFFIT-1)*(LOMOST+1)                             
      LSC = LSCSIG + (IOFFIT-1)*(LOMOST+1)*2                            
C                                                                       
C     FOR INITIAL AND FINAL FUNCITON REFS, PRINT STUFF                  
C                                                                       
      IF ( ITEMP .NE. 2 ) CALL PRTCHI ( WEIGHT, ASHIFT, ILLOC(LEHED+9),
     * RENORM )                                                         
C                                                                       
C     SETUP POTENTIAL GRADIENT ARRAY                                    
C                                                                       
      IF (GRADSW) CALL MAKDER ( NSTPSS(1)+1, NUMVS, MAXVS, MPARA,  XS, D
     *ERIVS, ALLOC(LVDERV), ALLOC(LDVTYP), ALLOC(Z(IRLVS(1))),  ALLOC(Z(
     *IIMVS(1))), ALLOC(LWAVR), ALLOC(LWAVI) )                          
C                                                                       
C     COMPUTE ALL S-MATRICES AND THEIR GRADIENTS                        
C                                                                       
      DO 100 LL = LLMIN, LLMAX, LSKIP                                   
C                                                                       
C  FOR NOW NO SPIN ORBIT -- ONLY ONE JVALUE                             
C                                                                       
      CALL WAVELJ(LL, 0, ICHANW, 0, DUMMM, DUMMM, DUMMM,  ALLOC(LWAVR), 
     *ALLOC(LWAVI), ALLOC(LVREAL), ALLOC(LVIMAG),  ALLOC(LCENTR) )      
C                                                                       
      IF (GRADSW) CALL SDERIV ( NFIRST, NSTPSS(1)+1, NUMVS, MAXVS,  MPAR
     *A, ALLOC(LWAVR), ALLOC(LWAVI), ALPHAR, ALPHAI,  ALLOC(LVDERV), ALL
     *OC(LDVTYP), ALLOC(LSDELS),  ALLOC(LDSDP+(2*MPARA)*LL) )           
C                                                                       
100   CONTINUE                                                          
C                                                                       
C     SDERIV HAS FOUND ONLY THE ANALYTIC DERIVATIVES                    
C     IF THERE ARE SOME NUMERIC ONES TO BE DONE; WE PERTURB             
C     THINGS AND DO THEM NOW.                                           
C                                                                       
      IF ( .NOT. GRADSW ) GO TO 150                                     
      II = 0                                                            
      DO 149 NX = 1, MPARA                                              
      IF ( ILLOC(FACFR4*LDVTYP-FACFR4+NX) .LT. 4 ) GO TO 149            
      II = 1                                                            
      PSAVE = XS(NX)                                                    
      DEL = FITACC*PSAVE                                                
      IF ( PSAVE .EQ.0.0D0) DEL = FITACC                                  
      XS(NX) = PSAVE+DEL                                                
      CALL SETFIT ( XS, IRET, 1 )                                       
      DO 139 LL = LLMIN, LLMAX, LSKIP                                   
      S1 = ALLOC(LSMAT+2*LL)                                            
      S2 = ALLOC(LSMAT+2*LL+1)                                          
      CALL WAVELJ ( LL, 0, ICHANW, 0, DUMMM,DUMMM,DUMMM,ALLOC(LWAVR),
     * ALLOC(LWAVI), ALLOC(LVREAL), ALLOC(LVIMAG),  ALLOC(LCENTR) )    
      ALLOC(LDSDP+2*NX+(2*MPARA)*LL-2)=(ALLOC(LSMAT+2*LL)-S1)/DEL  
      ALLOC(LDSDP+2*NX+(2*MPARA)*LL-1)=(ALLOC(LSMAT+2*LL+1)-S2)/DEL                                                                  
      ALLOC(LSMAT+2*LL) = S1                                            
      ALLOC(LSMAT+2*LL+1) = S2                                          
 139  CONTINUE                                                          
      XS(NX) = PSAVE                                                    
 149  CONTINUE                                                          
      IF ( II .EQ. 1 ) CALL SETFIT ( XS, IRET, 0 )                      
C                                                                       
 150  CONTINUE
C                                                                       
C  GET CROSS SECTIONS                                                   
C                                                                       
      IF ( ITEMP .NE. 2 ) WRITE(6,153) CHAR, CHAR                         
 153  FORMAT( 9X, 'ANGLE', A1, T25,  'EXPERIMENT', A1, T44, 'FIT', 
     * T61, 'FIT', T75, 'EXP ERROR', T89, '(FIT-EXP)/ERROR', T109,  '%/L
     * AT LMAX' /  9X, 'C.M.', T25, 'SIGMA/RUTH', T41,  'SIGMA/RUTH', T6
     *0, '(M.B.)', T79, '%', T91, '(UNWEIGHTED)' / )                    
      CHISQR=0.0D0                                                          
      SUM1=0.0D0                                                            
      SUM2=0.0D0                                                            
      SUM3=0.0D0                                                            
C                                                                       
C     DO ALL ANGLES AT THIS ELAB                                        
C                                                                       
      DO 589 IANG=1,NPTSE                                               
C                                                                       
C     GET THE DATA.  IT HAS BEEN STORED IN THE STANDARD FORM OF         
C       C.M. ANGLE,  C.M. SIGMA (MB),  % ERROR                          
C     ANGL - C.M. ANGLE IN DEGREES                                      
C     ANGLE - C.M. ANGLE IN RADIANS                                     
C             (BOTH OF ABOVE ARE MODIFIED BY LAB ANGLE SHIFT)           
C     ERR - ERROR IN PERCENT                                            
C     ERROR - WILL BE SET TO ERROR IN M.B. (INCLUDES EFFECT OF BOTH     
C             LAB ANGLE SHIFT AND RENORM.)                              
C     SIGMA - EXP. CROSS SEC IN M.B. (WILL BE MODIFIED BY LAB ANGLE     
C             SHIFT BUT NOT BY RENORM.)                                 
C     RSIGMA - RENORMALIZED EXP. CROSS SEC. IN M.B. (WILL INCLUDE       
C             BOTH LAB ANGLE SHIFT AND RENORM.)                         
C     CROSS - FITTED VALUE OF SIGMA (M.B.).                             
C     RUTH - RUTHERFORD CROSS SECTION IN M.B.                           
C                                                                       
C                                                                       
      ANGL=ALLOC4(LADATA)                                               
      ANGLE=ANGL*DEGREE                                                 
      ERR=ALLOC4(LADATA+1)                                              
      SIGMA=ALLOC4(LADATA+2)                                            
C                                                                       
C     IF WE ARE MAKING AN ANGLE SHIFT IN THE LAB, DO IT NOW             
C                                                                       
      IF (ASHIFT.EQ.0.0D0) GO TO 250                                    
C                                                                       
C     FIRST GET OLD JACOBIAN AND LAB ANGLE                              
C                                                                       
      TEMP =1.0D0 + RATMAS*( RATMAS + 2.0D0*DCOS(ANGLE) )                      
      RATJAC=(1.0D0+RATMAS*DCOS(ANGLE))/(TEMP*DSQRT(TEMP))              
      BACKSW = DCOS(ANGLE) .LT. -1.0D0/RATMAS                               
      ANGLE = DATAN2(DSIN(ANGLE),DCOS(ANGLE)+RATMAS )                 
C                                                                       
C     NOW SHIFT IN LAB AND CONVERT BACK TO C.M.  ALSO GET NEW JACOBIAN  
C                                                                       
      ANGLE = ANGLE + ASHIFT*DEGREE                                     
      TEMP = 1.0D0 - (RATMAS*DSIN(ANGLE))**2                                
      IF ( TEMP .LT.0.0D0 ) TEMP = 0.0D0                                       
      TEMP =DSQRT(TEMP)                                                 
      IF ( BACKSW ) TEMP = -TEMP                                        
      ANGL = DCOS(ANGLE)*TEMP - RATMAS*DSIN(ANGLE)**2                   
      ANGLE = ACOS(ANGL)                                                
      TEMP = 1.0D0 + RATMAS*(RATMAS+2.0D0*ANGL)                                 
      RATJAC = (1.0D0+RATMAS*ANGL)/(RATJAC * TEMP*DSQRT(TEMP))              
      ANGL = ANGLE*RADIAN                                               
C                                                                       
C     COMPUTE TRANSFORMED CROSS SECTION                                 
C                                                                       
      SIGMA = SIGMA*RATJAC                                              
C                                                                       
C     GET PREDICTED CROSS SECTION AND GRADIENT                          
C                                                                       
 250  CALL CALANG ( GRADSW, LLMIN, LLMAX, LSKIP, MPARA,  ANGLE,  ALLOC(L
     *PL), ALLOC(LSIG), ALLOC(LDFDX),  ALLOC(LSC), ALLOC(LSMAT), ALLOC(L
     *DSDP),  FR, FI, RUTH, HIGHR, HIGHI )                              
C                                                                       
C                                                                       
C  CROSS SECTIONS IN MILLI BARNS                                        
C                                                                       
      CROSS= 10.0D0 * (FR*FR + FI*FI)                                   
      RSIGMA = RENORM*SIGMA                                             
C                                                                       
C     CONVERT % ERRORS TO MB                                            
C                                                                       
      ERROR = .010D0*ERR*RSIGMA                                         
C                                                                       
C     COMPUTE DEVIATION                                                 
C                                                                       
      RATJAC = (CROSS-RSIGMA)/ERROR                                     
      FS(ICNTPT+IANG) = WT*RATJAC                                       
      CHISQ = CHISQ + (WT*RATJAC)**2                                    
C                                                                       
C     COMPUTE GRADIANT IF NEEDED                                        
C                                                                       
      IF ( .NOT. GRADSW ) GO TO 450                                     
      RAT = 20.0D0*WT/ERROR                                                 
      DO 439 I = 1, MPARA                                               
      KEY = ILLOC(LMASK+I-1)                                            
      IF ( KEY .LT. 0 ) GO TO 410                                       
      TEMP = RAT*( FR*ALLOC(LDFDX+2*I-2)  + FI*ALLOC(LDFDX+2*I-1) )     
      GO TO 435                                                         
C                                                                       
C     THE PARAMETER IS A RENORMALIZATION OR ANGLE SHIFT -               
C     IT MAY NOT EFFECT THIS DATA SET.                                  
C                                                                       
 410  TEMP = 0.0D0                                                          
      ISET = MOD( -KEY/10000, 100 )                                     
      ITYPE = -KEY/100 0000                                             
      IF ( ISET .EQ. IOFFIT ) GO TO 425                                 
      INDEX = LEQUAL                                                    
 415  IMANY = ILLOC(INDEX)                                              
      IF ( IMANY .EQ. 0 ) GO TO 435                                     
      IF ( ILLOC(INDEX+1) .EQ. KEY ) GO TO 420                          
      INDEX = INDEX + IMANY+2                                           
      GO TO 415                                                         
 420  DO 424 IM = 1, IMANY                                              
      IF(MOD( -ILLOC(INDEX+1+IM)/10000, 100 ).EQ.IOFFIT ) GO TO 425
 424  CONTINUE                                                          
      GO TO 435                                                         
C                                                                       
C     IT DOES EFFECT THIS GUY.                                          
C    WE CAN COMPUTE THE GRADIENT ONLY FOR RENORMALIZATION               
C                                                                       
 425  IF ( ITYPE .EQ. 0 ) TEMP = -WT*CROSS/(ERROR*RENORM)               
C                                                                       
C     STORE GRADIENT OF F AND COMPUTE GRADIENT OF CHI SQUARE            
C                                                                       
 435  FGRAD(ICNTPT+IANG,I) = TEMP                                       
      TGRAD(I) = TGRAD(I) + 2*TEMP*(WT*RATJAC)                          
  439  CONTINUE                                                          
C                                                                       
 450  IF(ITEMP.EQ.2)GO TO 580                                           
C                                                                       
C     WE PRINT SIGMA/RUTH                                               
C                                                                       
      RSIGMA = RSIGMA/RUTH                                              
      RCROSS = CROSS/RUTH                                               
      PERCEN=(1.0D0-(HIGHR*HIGHR+HIGHI*HIGHI)*10.0D0/CROSS) * 50        
      WRITE(6,501)ANGL,RSIGMA,RCROSS,CROSS,ERR,RATJAC,PERCEN            
501   FORMAT ( F15.4, 2X, 5G17.4, F14.2 )                               
      CHISQR=CHISQR+RATJAC**2                                           
      ERR=ERROR*ERROR                                                   
      SUM1=SUM1+CROSS*SIGMA/ERR                                         
      SUM2=SUM2+SIGMA*SIGMA/ERR                                         
      SUM3=SUM3+CROSS*CROSS/ERR                                         
C                                                                       
C     ADVANCE TO NEXT ANGLE                                             
C                                                                       
 580  LADATA = LADATA + LENDAT                                          
C                                                                       
589   CONTINUE                                                          
C                                                                       
C     END OF THE LOOP OVER ANGLES FOR THIS ELAB                         
C                                                                       
      IF(ITEMP.EQ.2)GO TO 780                                           
C                                                                       
C     PRINT A SUMMARY OF THIS DATA SET                                  
C                                                                       
      IF ( CHAR .EQ. STAR ) WRITE(6,623) ASHIFT, RENORM                   
 623  FORMAT ( ' * INCLUDES LAB-ANGLE SHIFT OF', F7.3,  ' DEGREES AND RE
     *NORMALIZATION OF', F7.4 )                                         
      THISCN = WEIGHT * CHISQR/NPTST                                    
      CHISQR=CHISQR/NPTSE                                               
      WCHISQ = WEIGHT*CHISQR                                            
C                                                              
C     ROPT IS A FACTOR TO SCALE THE DATA TO THE FIT:                    
C       BEST DATA = ROPT * DATA                                         
C       BEST FIT  = FIT / ROPT                                          
C                                                                       
      ROPT=SUM1/SUM2                                                    
      FOPT = (ROPT*ROPT*SUM2-2.0D0*ROPT*SUM1+SUM3)/NPTSE                    
      WRITE(6,633) NPTSE, RENORM, CHISQR, WCHISQ, ROPT,FOPT             
 633  FORMAT ( X,I5, ' POINTS' /  ' RENORMALIZATION =',F8.4,'  GIVES 
     *UNWEIGHTED CHI.SQR./PT. =', G12.4, 5X, 'WEIGHTED CHI.SQR./PT. ='
     *,G12.4 /  ' OPTIMUM RENORMALIZATION =', G12.4,  5X, 'OPTIMUM 
     *UNWEIGHTED CHI.SQR./PT. =', G12.4 )                                   
      IF ( IESETS .GT. 1 ) WRITE(6,643) THISCN                            
 643  FORMAT (' CONTRIBUTION OF THIS DATASET TO TOTAL CHI SQ/PT =', 
     *G12.4 )                                                             
C                                                                       
C     END OF THIS ELAB, GO TO THE NEXT ONE                              
C                                                                       
 780  ICNTPT=ICNTPT+NPTSE                                               
      LEKEY = LADATA-1                                                  
 789  CONTINUE                                                          
C                                                                       
C     THE FUNCTION VALUE IS THE CHI-SQ SUM.                             
C                                                                       
      CALFUN = CHISQ                                                    
      CHISQ = CHISQ/NPTST                                               
C                                                                       
      IF ( PRNTSW ) WRITE(6,803) CHISQ                                    
 803  FORMAT ( ' CALFUN OUTPUT =', G20.10 )                             
C                                                                       
      IF(ITEMP.EQ.2 )GO TO 950                                          
C                                                                       
      CHISQR = CALFUN/(NPTST-MPARA)                                     
      WRITE(6,813) MPARA, NPTST, CHISQ, CHISQR                          
 813  FORMAT ( X, I3, ' PARAMETERS,',  I5, ' POINTS,    TOTAL CHI.SQR.
     */PT. =', G12.4,  5X, 'TOTAL CHI.SQR./DEG.FREEDOM =', G12.4 // )   
      GO TO 950                                                         
C                                                                       
C     ATTEMPTED TO STEP OUT OF BOUNDS                                   
C                                                                       
 900  CALFUN = 0.0D0                                                        
      DO 959 I = 1, NPTST                                               
      FS(NPTST) = 1.D+10                                                
      CALFUN = CALFUN + 1.D+20                                          
 959  CONTINUE                                                          
 950  CONTINUE
CDC      CALFN2 = CALFUN                                                
      RETURN                                                            
C                                                                       
CDC      END                                                            
CDC      FUNCTION CALFUN ( NPTST, MPARA, FS, XS )                       
CDC      IMPLICIT                                                       
CDC      ALLOCLEVEL FS,XS                                               
CDC      DIMENSION FS(1), XS(1)                                         
CDC      CALFUN = CALFN2 ( NPTST, MPARA, FS, XS, A, B )                 
CDC      RETURN                                                         
CDC      END                                                            
CDC      FUNCTION CALGRA ( NPTST, MPARA, FS, XS, FT, TG )               
CDC      IMPLICIT                                                       
CDC      DIMENSION FS(1), XS(1), FT(1), TG(1)                           
CDC      ALLOCLEVEL FS, XS, FT, TG                                      
CDC      CALGRA = CALGR2( NPTST, MPARA, FS, XS, FT, TG )                
CDC      RETURN                                                         
C                                                                       
      END                                                              

      SUBROUTINE GETIT(IRET)
      RETURN
      END
      SUBROUTINE DUMMY2
      RETURN
      END
      SUBROUTINE CKSET
      RETURN
      END

      SUBROUTINE EPSLON(X,N,FIN,DEPS,DERR,IER)
      IMPLICIT REAL*8(A-H,O-Z)
      COMPLEX*16 X(N),FIN
      COMPLEX*16 T,W1,W2,W3,W4,W5,W6,W7
      COMPLEX*16 Z
      LOGICAL ISW1,ISW2
      COMMON/CNSTNT/ PI,RT4PI,PIINV,RADIAN,DEGREE,HBARC,
     & AMUMEV,AFINE,BIGEST,SMLEST,BIGNUM,SMLNUM,BIGLOG
      LOGICAL NONZER,ZERO
      APXABS(Z)=DABS(DREAL(Z))+DABS(DIMAG(Z))
      NONZER(Z)=APXABS(Z).NE.0.0D0
      ZERO(Z)=APXABS(Z).EQ.0.0D0
      IER=-1
      IF(N.LT.5)RETURN
      BIG=1.D+38
      ACC=DMAX1(1.D-8,DEPS**2)
 41   DERR=0.0D0
      FIN=X(N)
      DO 43 I=1,N
      DERR=DMAX1(DERR,APXABS(FIN-X(I)))
 43   CONTINUE
      DERR=DERR/(APXABS(FIN)+SMLNUM)
      IER=0
      IF(DERR.LE.ACC)RETURN
      ACC=DEPS
      IF(N.GE.6)GO TO 2
      IER=1
      RETURN
2     ISW1=.FALSE.
      ISW2=.FALSE.
      W1=BIG
      W7=X(4)-X(3)
      IF(NONZER(W7))W1=1.0D0/W7
      W5=BIG
      W7=X(2)-X(1)
      IF(NONZER(W7))W5=1.0D0/W7
      W4=X(3)-X(2)
      IF(NONZER(W4))GO TO 9
      W4=BIG
      T=X(2)
      W2=X(3)
 8    W3=BIG
      GO TO 17
 9    W4=1.0D0/W4
      T=BIG
      W7=W4-W5
      IF(NONZER(W7))T=X(2)+1.0D0/W7
      W2=W1-W4
      IF(NONZER(W2))GO TO 15
      W2=BIG
      ISW2=DREAL(T).NE.BIG
      W3=W4
      GO TO 17
 15   W2=X(3)+1.0D0/W2
      W7=W2-T
      IF(ZERO(W7))GO TO 8
      W3=W4+1.0D0/W7
 17   ISW1=ISW2
      ISW2=.FALSE.
      IMIN=4
      DO 40 I=5,N
      IAUS=I-IMIN
      W4=BIG
      W5=X(I-1)
      W7=X(I)-X(I-1)
      IF(ZERO(W7))GO TO 24
      W4=1.0D0/W7
      IF(DREAL(W1).EQ.BIG)GO TO 25
      W6=W4-W1
      IF(APXABS(W6).GT.1.D-12*APXABS(W4))GO TO 22
      ISW2=.TRUE.
      IF(NONZER(W6))GO TO 22
      W5=BIG
      W6=W1
      IF(DREAL(W2).NE.BIG)GO TO 28
      GO TO 26
  22  W5=X(I-1)+1.0D0/W6
      IF(APXABS(W5).GE.1.D-10*APXABS(X(I-1)))GO TO 24
      IF(NONZER(W5))GO TO 36
 24   W7=W5-W2
      IF(NONZER(W7))GO TO 27
 25   W6=BIG
 26   ISW2=.FALSE.
      X(IAUS)=W2
      GO TO 37
 27   W6=W1+1.0D0/W7
 28   IF(.NOT.ISW1)GO TO 33
      IF(DREAL(W2).EQ.BIG)GO TO 32
      W7=W5/(W2-W5)+T/(W2-T)+X(I-2)/(X(I-2)-W2)
      IF(ZERO(W7+1.0D0))GO TO 38
      X(IAUS)=W7*W2/(1.0D0+W7)
      GO TO 39
 32   X(IAUS)=W5+T-X(I-2)
      GO TO 39
 33   W7=W6-W3
      IF(ZERO(W7))GO TO 38
      X(IAUS)=W2+1.0D0/W7
      IF(APXABS(X(IAUS)).GE.1.D-10*APXABS(W2))GO TO 37
      IF(NONZER(X(IAUS)))GOTO 36
37    IF(DREAL(W2).NE.BIG)GO TO 39
 38   X(IAUS)=BIG
      IMIN=I
 39   W1=W4
      T=W2
      W2=W5
      W3=W6
      ISW1=ISW2
40    ISW2=.FALSE.
      N=N-IMIN
      GO TO 41
 36   N=IAUS-1
      GO TO 41
      RETURN
      END
