C @(#)dlsq.for	17.1.1.1 (ES0-DMD) 01/25/02 17:17:16
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C  @(#)dlsq.for	17.1.1.1  (ESO-IPG)  01/25/02  17:17:16
      SUBROUTINE DLSQ(YP,NPTS,NPARMS,NBLOX,MBLOK,NIX,NARGS,NKARGS,TEST)
C
C  14 OCT 1985 VERSION.    **********     USES PARTITIONING.
C
C       Copyright (C) Andrew T. Young, 1990, 1991
C
C  Fits data to nonlinear models by Damped Least Squares.
C  (This version uses robust estimation.)
C
C
C   ARGUMENTS:
C
C  YP      is the subroutine called to calculate YT and PARTial derivs.
C             --> MUST be declared EXTERNAL in calling pgm!
C  NPTS    is the number of Data POINTS in the problem.
C             --> must be .LE. MXPT.
C  NPARMS  is the Number of PARAMETERS actually adjusted.
C             --> must be .LE. MXP.
C  NBLOX   is the Number of BLOCKS in the block-diagonal (first)
C             part of the normal eqns.
C  MBLOK   is the Number of PARAMETERS in a Block.
C             --> (NBLOX*MBLOK) must be .LE. NPARMS.
C  NIX     is the Number of Parameters that are held FIXed.
C             --> must be .LE. NPARMS.
C  NARGS   is the Number of REAL independent variables actually used.
C             --> must be .LE. MXX.
C  NKARGS  is the Number of INTEGER independent variables actually used.
C             --> must be .LE. MKX.
C  TEST    When the relative change in every parameter is less than
C             TEST, we are done.
C
C
        IMPLICIT DOUBLE PRECISION (A-H,O-Z)
        IMPLICIT INTEGER (I-N)
C
        EXTERNAL YP
        REAL TEST
C
      INCLUDE 'MID_REL_INCL:kingfit.inc'
C     PARAMETER (MXPT=xxx, MXX=1, MKX=x, MXP=xxx)
C
C       MXPT = MAXimum number of data POINTS.
C
C       MXP  = MAXimum number of Parameters.
C
C       MXX  = MAX.no.of real independent variables.
C
C       MKX  = MAX.no.of integer independent variables.
C
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C       Defines common/LEASTC/; Equiv's (KDUM,MAXK etc.)
C
C
C  The above COMMON & PARAMETER statements MUST be in both the subroutine
C  YP, which the user supplies, and the main program that uses these
C  subroutines.  ****  The ORDERING of  /LEASTC/  MUST be identical.
C
C
C  Here are the variables in /LEASTC/:
C
C   A    is the matrix of normal equations, or (on exit) its inverse.
C  DUM   is a DUMmy set of Real numbers carried throughout the problem.
C        These may be used to transmit information to YP directly
C        from the calling program.  May be removed or redimensioned.
C  KDUM  is a DUMmy set of integers (same purpose).
C  IDF   is the number of Degrees of Freedom, not counting points
C           with zero weight.
C  KIP   is the number of active Partials for each datum.  (set in YP.)
C  IP    Indexes the active Partials for each datum.  (set in YP.)
C  IPR   is the print flag.  IPR > 0 produces output on the progress of
C             the iterations.
C          IPR = 1 gives WVAR, DAMP, & MODE; IPR = 2 gives P's as well;
C          IPR = 3 adds internal state & SCALE; IPR=7 checks derivatives.
C  IT    is the ITeration number.  (limited to MAXIT.)
C  IX    are the 'serial numbers' of parameters that are held FIXed
C           for a particular problem.  Only the first NIX of the
C           IX's are used  (see example below).
C  MAXP  is the MAX.no.of P's being used (=NPARMS, to YP via /LEASTC/).
C  NEQS  is the number of equations being used ( = NPARMS-NIX ).
C  NONEQ tells subroutine YP whether partial derivatives are needed.
C           NONEQ = 0 if needed,   NONEQ =1 if "NO Normal EQuations".
C  NPT   is the Number of the datum Y() being processed in YP.
C  PG    are Preliminary Guesses (starting values).
C  P     are the values of the PARAMETERS being solved for.
C  PART  are the PARTial derivatives of YT with respect to the P's.
C  SP    are the Standard deviations of the Parameters.
C  W     are the Weights associated with the  Y  data.
C  WVAR  is the Weighted VARiance.  STD.ERR.(unit wt) = sqrt(WVAR).
C  X     are the Real  INDEPENDENT  variables.  The I th value of the
C           J  th indepedent variable is  X(J,I).
C  KX    are the Integer  INDEPENDENT  variables (similar indexing).
C  Y     are the  DEPENDENT  variables. (data being fit)
C  YC    are the values of Y Calculated for the
C           current values of the parameters, P.
C  YT    is the value of the fitted fcn.calc.by YP for arguments Z, KZ.
C  Z(J)  = X(J,NPT).
C KZ(J)  = KX(J,NPT).
C
C  Example of use of IX:  to fix the 2nd,5th, and 6th parameters,
C              set NIX=3, IX(1)=2, IX(2)=5, and IX(3)=6.
C  NOTE: Order is unimportant.  If all P are held fixed, DLSQ gives
C        residuals.
C
C
C   Use of KIP, IP, etc. and the YP subroutine:
C
C   YP generates the value of YT for each equation of condition.
C      It also generates the values of KIP partial derivatives, PART.
C      In block-diagonal problems, these are only a small subset of the
C      total number of partials.  The ones that are generated are listed
C      in the IP array; e.g., if the 4 partials PART(1), PART(2), PART(5),
C      and PART(9) are non-zero for a given datum, then set KIP=4 and
C      IP(1)=1, IP(2)=2, IP(3)=5, and IP(4)=9 in YP.  The order is not
C      important -- it's just a list to tell this subroutine which PARTs
C      are non-zero.  This saves a lot of time in building the normal
C      equations in sparse, block-diagonal problems.  YP may generate
C      other, unused partials, but they should not be returned to this
C      program in the IP list.
C
C   Note that MAXP (= NPARMS) and NEQS are both available in /LEASTC/
C   so that YP can decide which PARTs to index in IP.  Note also that
C   the block to which a given datum contributes (in block-diagonal
C   problems) is usually inferred from the KZ (KX in calling pgm.)
C   integer variables.
C
C
C   There is a conflict between this and earlier versions of DLSQ.
C   To preserve the block-diagonal format of the normal equations, we
C   have to carry ALL the P's through the solution, even if some are
C   kept fixed.  This is handled here by NOT indexing the normal eqns.
C   The equations corresponding to a non-zero IX(n) are simply replaced
C   by a 1 on the main diagonal, and zeros on RHS (i.e., in B).
C
C   This makes this version less efficient than the older one for simple
C   problems.  Also, we lose more efficiency here in problems with NBLOX=1
C   because we explicitly invert the matrix.  But for little problems, who
C   cares?
C
C   Note that this version stores the input data (X, Y, W, and PG) as REAL,
C   but does the matrix work (A, B, P, YT, YC, WVAR) in DOUBLE PRECISION.
C
C
C   LOCAL variables:
C
      DIMENSION B(MXP),DP(MXP),PC(MXP),PLOW(MXP)
      DIMENSION WT(MXPT)
      REAL S1(MXPT)
C
C   B    is the RHS vector (in the normal equations: AX = B ).
C   DP   is the change to P(I) on each iteration.
C   PC   = P(I) + DP(I)  is the Corrected value of P.
C   PLOW is the set of P's that gave the LOWest value of WVAR in the past.
C   USED is .TRUE. if the corresponding P is *not* in the IX list.
C   USED is .FALSE. if the corresponding P *is* in the IX list.
C   PAGE is used to store the printed version of matrix rows.
C   S1   is used to save residuals for sorting.
C
C
      LOGICAL PARTSW, USED(MXP), HUBER
      CHARACTER*79 CARD,PAGE(MXP/4+1)
C
      DATA PARTSW/.FALSE./,MAXIT1/99/
C
C  ***********************  INITIALIZE EVERYTHING  ******************
C
      MODE=0
      NZD=0
      K=0
      IT=0
      IRESET=0
      LAST4=0
      NONEQ=0
C  Compute variance but not normal equations if NONEQ.NE.0.
      WVAR=0.
      RESETS=0.
      RESTART=0.
      SCALE=0.
      IF (MAXIT.EQ.0.) MAXIT=MAXIT1
      OLDVAR=-1.
      DAMP=1.E-4
      FACTOR=100.
      HUBER=.TRUE.
      RTTEST=SQRT(TEST)
      MAXP=NPARMS
      NEQS = NPARMS - NIX
C   NEQS=No. of parameters actually solved for.
      FNEQS=NEQS
C
C     Dummy write to space o/p.
      IF (IPR.GT.0) CALL SPACE2
      DO 80 I=1,NPARMS
      DP(I)=0.0
      PART(I)=0.0
      P(I)=PG(I)
      PC(I)=PG(I)
      PLOW(I)=PG(I)
      DO 75 J=1,NIX
         IF (IX(J).LE.0 .OR. IX(J).GT.NPARMS) THEN
            WRITE(CARD,'(A,I4,A,I4)')'Invalid entry in IX(',J,')=',IX(J)
            CALL TV(CARD)
            GO TO 80
         else IF(I.EQ.IX(J)) THEN
            USED(I)=.FALSE.
            DO 70 KK=J+1,NIX
C              make sure there are no dupes...
               IF (I.EQ.IX(KK)) THEN
                  WRITE(CARD,'(I4,A,I4,A)')i,' Duplicated in IX(',J,')'
                  CALL TV(CARD)
                  GO TO 80
               END IF
   70       CONTINUE
C           OK, I is a legal IX and is not repeated in rest of IX list.
            K=K+1
            GO TO 80
         END IF
   75 CONTINUE
      USED(I)=.TRUE.
   80 CONTINUE
      IF (NPTS.LT.NEQS) THEN
          CALL TV('More unknown parameters than data!')
          RETURN
      END IF
C
C  NOTE:
C    --> YP *MUST* return KIP and IP, as well as YT and PART's.
C
C  K should equal NIX, but is less if IX contains duplicates.
C
      IF(K.LT.NIX) THEN
          WRITE(PAGE,85)NIX,K,NPARMS,IX
   85 FORMAT('WARNING: The first',I4,' entries of IX contain only',I4,
     1' unique values .LE.',I4/(10I4/))
          CALL TV(PAGE(1))
          DO 86 I=2,NPARMS/4+1
   86     CALL TVN(PAGE(I))
          CALL STSEPI
      END IF
C
C  **********************************************   END SETUP  ******
C
      NONEQ = 1
C  NONEQ=1 suppresses partials in YP & normal eqns.
      IF (NEQS.EQ.0) GO TO 120
C
C  INITIALIZE FOR NORMAL EQNS.  **************BEGIN ITERATION *******
  100 IT=IT+1
C$$$
      IF(IT.GE.MAXIT)GOTO 905
C$$$
C  IT is the iteration counter.
C
  101 K=0
      NREJ=0
      NUSED=0
      IF (HUBER) THEN
         CTUNE=1.5
      ELSE
         CTUNE=6.
         SCLNUM=0.
         SCLDEN=0.
      END IF
      CS=CTUNE*SCALE
      DO 105 I=1,NPARMS
      B(I)=0.
      DO 105 J=1,I
      K=K+1
  105 A(K)=0.
C
  120 VAR=0.
C  NOW FORM NORMAL EQNS.  **********  BEGIN NORMAL EQUATIONS.********
C
C  Form EQ.OF COND. for NPT-th point.
      DO 129 NPT=1,NPTS
      DO 121 J = 1, NARGS
  121 Z(J) = X(J,NPT)
      DO 122 J=1,NKARGS
  122 KZ(J)=KX(J,NPT)
C
      CALL YP
C
C   note that all values are transmitted via COMMON /LEASTC/
C
      YC(NPT) = YT
      IF(W(NPT).EQ.0.) GO TO 129
C  Skip zero-weight points.   ******** ACCUMULATE VARIANCES. ********
      DY = Y(NPT) - YT
      IF (SCALE.NE.0.) THEN
C         normal iteration.
          U=DY/CS
C         CS valid for unit weight, so include weight here.
          USQ=U*U*W(NPT)
C
          IF(MODE.EQ.0) THEN
C             adjust weights.
              IF (HUBER) THEN
                  IF (USQ.GE.1.) THEN
                      WT(NPT)=W(NPT)/USQ
                      NREJ=NREJ+1
                  ELSE
                      WT(NPT)=W(NPT)
                  END IF
              ELSE
C                 Compute biweight:
                  IF (USQ.GE.1.) THEN
                      WT(NPT)=0.
                      GO TO 129
                  ELSE
                      WT(NPT)=W(NPT)*(1.-USQ)**2
                  END IF
              END IF
          END IF
C
          IF (WT(NPT).EQ.0.) GO TO 129
          VAR=VAR+WT(NPT)*DY*DY
C         prepare to reset scales.
          NUSED=NUSED+1
          IF (HUBER) THEN
C            use MAD:
             S1(NUSED)=W(NPT)*DY*DY
          ELSE
C            do sums for A-estimator:
             SCLNUM=SCLNUM + WT(NPT)*DY*DY*(1.-USQ)**2
             SCLDEN=SCLDEN + (1.-USQ)*(1.-5.*USQ)
          END IF
      ELSE
C         initial iteration.
          NUSED=NUSED+1
          S1(NUSED)=W(NPT)*DY*DY
      END IF
C  Skip normal equations if not needed.
      IF(NONEQ.NE.0) GO TO 129
      NULKIP=0
C       KIP and IP are set by YP.
      DO 123 J=1,KIP
          I1=IP(J)
          IF (.NOT.USED(I1) .OR. PART(I1).EQ.0.) THEN
C           shrink IP stack.
            NULKIP=NULKIP+1
          ELSE
            IF(NULKIP.NE.0) IP(J-NULKIP)=IP(J)
          END IF
  123 CONTINUE
C  Add CONTRIBUTION of NPT-th point to NORMAL EQUATIONS.
      DO 125 I=1,KIP-NULKIP
      I1=IP(I)
C ***  IF(kindex(I1).ne.0. .and. PART(I1).NE.0.)THEN
      TERM=PART(I1)*WT(NPT)
      B(I1)=B(I1)+TERM*DY
      DO 124 J=1,I
          J1=IP(J)
C *** if (kindex(J1).eq.0.) go to 124
          K= (MAX(J1,I1) -1) * MAX(J1,I1) /2
     1       +MIN(J1,I1)
          A(K)=A(K)+TERM*PART(J1)
  124 CONTINUE
C ***  END IF
  125 CONTINUE
C  The above is very efficient.  DON'T MESS!
  129 CONTINUE
C
C     Make sure VAR is OK:
      IF (VAR.GE.0. .AND. VAR.LT.3.E33) THEN
C        VAR is OK.
      ELSE
C        VAR is bad.  Force reset:
         IF (IPR.GT.0)
     1      WRITE (CARD,'(''VAR has unreasonable value: '',E8.2)')VAR
         CALL TV(CARD)
         VAR=CAP+1.
         GO TO 136
      END IF
C
      IF (SCALE.EQ.0.) THEN
C         set initial scale to 1.4*median of rough estimates.
          IF (NUSED.EQ.0) THEN
             CALL TV('NUSED=0.  Probably all data have weight 0.')
             CALL STSEPI
          END IF
          CALL SORT1(S1,NUSED)
          SCALE=SQRT(S1((NUSED+1)/2) + S1(NUSED/2+1))
          IF (NEQS.GT.0) NONEQ=0
          IF (SCALE.LE.0.) THEN
              SCALE=1.
              IF (IPR.GT.0) CALL TV('SCALE fell to zero at 129.')
          END IF
          GO TO 101
      END IF
C
      IF(IPR.LT.2 .OR. MODE.NE.5) GO TO 135
C
C  PRINT normal equations on final pass.....
      CALL SPACE2
      WRITE(CARD,131)
  131 FORMAT(3X,'I',30X,'A(I,J)',30X,'B(I)')
C       formats 131 - 133 changed to fit on 80-col.screen.
      CALL TV(CARD)
      CALL SPACE2
      K=0
      DO 134 I=1,NPARMS
         IF (.NOT.USED(I)) GO TO 134
         WRITE(PAGE,132)I,(A(J),J=K+1,K+I)
  132    FORMAT(I4,1P5E15.5/(E19.5,4E15.5))
         CALL PRPAGE(PAGE,(I+4)/5)
         WRITE(CARD,133) B(I)
  133    FORMAT(63X,1PE16.7)
         CALL TVN(CARD)
         CALL SPACE2
         K=K+I
  134 CONTINUE
C     Print a row of *'s:
      WRITE(CARD,555)
      CALL TV(CARD)
C
  135 IF(NEQS.EQ.0 .OR. MODE.EQ.5) GO TO 800
        IF(MODE.EQ.6) GO TO 498
C                          ****** NORMAL EQUATIONS COMPLETE.   ******
      IF(OLDVAR.EQ.-1.)THEN
         ABSMIN=VAR
         CAP=VAR+VAR
         GO TO 230
      END IF
C                                ********** TEST VARIANCES **********
C    ***
  136   IF(VAR.LT.CAP)THEN
C    ***
C  WELL-BEHAVED.
        IF (VAR.EQ.0.) GO TO 800
C   --
       IF(VAR.LE.OLDVAR .AND. NREJ*5.LT.NPTS)THEN
C   --
C   previous iteration went downhill, safely.
      MODE=2
      IF(DAMP.EQ.0.) NZD=NZD+1
C
        IF(VAR.LT.ABSMIN)THEN
C    we are at a new low.
          DO 139 I=1,NPARMS
  139     PLOW(I)=P(I)
          CAP=CAP+VAR
          ABSMIN=VAR
          IF(DAMP.NE.0.) FACTOR=FACTOR*1.2
          MODE=4
          LAST4=IT
        END IF
C
          IF (MODE.EQ.4) THEN
             DAMP=DAMP/FACTOR
          ELSE IF (MODE.EQ.2) THEN
             IF(DAMP.EQ.0. .AND. NZD.GE.5+IT/20)THEN
               NZD=0
C              Restore damping if undamped not converging.
               DAMP=TEST
            END IF
            DAMP=DAMP*(RESETS+1.)/(RESETS+2.)
          END IF
C   --
       ELSE
C   --
C    if many rejects, continue without adjustment.
         IF (VAR.LT.OLDVAR) THEN
           MODE=3
           GO TO 231
         END IF
C    previous iteration went uphill, but safe.
         IF(VAR.GT.SQRT(ABSMIN*CAP)) DAMP=(DAMP+1.E-9)*(FACTOR+FACTOR)
         FACTOR=SQRT(FACTOR)+1.
         MODE=1
         IF(LAST4.EQ.0) DAMP=DAMP*2.
         DAMP=DAMP*10.**(MAX(0,(IT-LAST4)+4/(IT+1-IRESET))/6) + 1.E-9
C   --
       END IF
C   --
C   prevent endless wandering above absmin.
      IF(MODE.NE.4) CAP=ABSMIN+CAP/2.
C    ***
        ELSE
C    ***
C  VAR.gt.CAP, or dangerous number of rejects.
C
      RESETS=RESETS+1.
      FACTOR=1.+SQRT(FACTOR)
          IF (CAP.LE.0. .OR. RESETS.GT.100.) THEN
             DO 205 NPT=1,NPTS
                IF (W(NPT).LT.0.) THEN
                   CALL TV('Negative weight -- PROGRAM ERROR !!')
                   CALL STSEPI
                END IF
  205        CONTINUE
             IF (IPR.GT.0) CALL TV('Weights checked -- OK.')
          END IF
                        IF(IPR.GT.0)THEN
                           WRITE(CARD,*)CAP,RESETS,NREJ,SCALE
                           CALL TV(CARD)
                           WRITE(CARD,*) FACTOR, DAMP
                           CALL TVN(CARD)
                        END IF
        IF(DAMP.EQ.0.)THEN
          DAMP=1.E-4*FACTOR
          IF (IPR.GT.0) CALL TV('DAMP reset to 1.E-4')
        ELSE IF(IRESET.EQ.IT) THEN
C         Reset loop.  Raise CAP.
               MODE=0
               OLDVAR=-1.
          CAP=VAR+2.*(VAR-ABSMIN)
          IF (CAP.LE.0.) CAP=1.+RESETS
                  CALL TV('Reset looping.')
          DAMP=FACTOR*(DAMP+1.E-4)
          IF (DAMP.GT.1./TEST) GO TO 905
        ELSE IF(RESETS.GT.1.) THEN
          DAMP=DAMP*FACTOR*RESETS
          CAP=(CAP+ABSMIN+ABSMIN)/2.
        ELSE
C         Big damp on 1st reset.
          DAMP=FACTOR
C         sets damp to 1.
        END IF
C
C  Restart from best previous position.          *****  RESET  *****
      OLDVAR=ABSMIN
      IF(IPR.GT.0) THEN
         WRITE(CARD,208)
  208    FORMAT(68X,'RESET *****')
         CALL TV(CARD)
      END IF
      IRESET=IT
      DO 210 I=1,NPARMS
  210 P(I)=PLOW(I)
      GO TO 101
C    ***
        END IF
C    ***
  230 OLDVAR=VAR
  231 CONTINUE
C
C  ****************************  NOW SOLVE NORMAL EQUATIONS. ********
      DO 490 I=1,NPARMS
      K=I*(I+1)/2
C       Diagonal element: ignore missing parms...
      IF (A(K).EQ.0.D0) A(K)=1.D0
  490 CONTINUE
       IF((1.+DAMP).EQ.1.)THEN
      DAMP=0.
      GO TO 500
       END IF
C                       Check derivatives if IPR=7:
                         IF (IPR.EQ.7) GO TO 905
      IF(DAMP.LE.1./TEST) GO TO 498
      WRITE(CARD,491)DAMP
  491 FORMAT('EXCESSIVE DAMPING',1PE10.1,
     1' possibly due to bad derivatives in YP,')
      CALL TV(CARD)
      CALL TVN('very bad starting values, or unsuitable model.')
      GO TO 905
C
  495 WRITE(CARD,496)DAMP,IT
  496 FORMAT('DAMP =',1PE10.2,' gave SINGULAR MATRIX on iteration',I3)
      CALL TV(CARD)
      IF (IPR.GT.0) THEN
         WRITE(CARD,*)'PARTIT was called with NEQS=',NEQS,
     1                     ' NBLOX=',NBLOX,' MBLOK=',MBLOK
         CALL TV(CARD)
      END IF
      IF(MODE.EQ.5) GO TO 905
      DAMP=(DAMP+1.E-9)/1.E-4
C   force reset.
      VAR=CAP+CAP+1.E-4
      GO TO 136
C
C       Apply damping...
  498 DO 499 I=1,NPARMS
      K=I*(I+1)/2
  499 A(K)=A(K)*(1.+DAMP)
        IF(MODE.EQ.6) GO TO 800
C
  500 CALL PARTIT(A,NPARMS,SP,NBLOX,MBLOK,*495)
C  use SP as SCRATCH SPACE to invert A, so X = (A-1)*B.
      K=0
C   FORM X IN SP.
      DO 520 I=1,NPARMS
      SP(I)=0.
C    FIRST PART.
      DO 510 J=1,I
  510 SP(I)=SP(I)+A(K+J)*B(J)
      K=K+I
C    SECOND PART.
      KOL=K+I
      DO 512 J=I+1,NPARMS
      SP(I)=SP(I)+A(KOL)*B(J)
  512 KOL=KOL+J
  520 CONTINUE
C
C     *****   NORMAL EQNS.NOW SOLVED, SOLN.IN SP.   REVISE PC'S. ******
C  UPDATE P'S.
      DO 540 I=1,NPARMS
      DP(I) =SP(I)
      PC(I) = P(I) + DP(I)
  540 Continue
C                                   **********  PC'S NOW SET. *******
      IF(IPR.EQ.0) GO TO 580
      WRITE(CARD,550) IT, VAR, SCALE, MODE
  550 FORMAT('IT =',I3,6X,'WVAR =',1PE22.15,6X,'SCALE =',E9.2,6X,
     1'MODE =',I2)
      CALL TV(CARD)
      IF (IPR.GT.1) THEN
         WRITE(CARD,551) FACTOR, CAP, DAMP
  551    FORMAT(13X,'FACTOR =',F6.2,5X,'CAP =',1PE9.2,5X,'DAMP =',E10.3)
         CALL TVN(CARD)
      END IF
      IF(IPR.EQ.1) GO TO 580
      WRITE(CARD,552)
  552 FORMAT('VAR.NO',6X,'PG',13X,'P',14X,'PC',13X,'DP')
      CALL TV(CARD)
      CALL SPACE
      DO 553 J=1,NPARMS
      IF (IPR.LT.3 .AND. .NOT.USED(J)) GO TO 553
      I=J
      WRITE(PAGE,132) I,PG(I),P(I),PC(I),DP(I)
      CALL PRPAGE(PAGE,1)
  553 CONTINUE
C
C  omit fixed P's unless this is last iteration.
       IF(MODE.EQ.5 .AND. NIX.NE.0)THEN
      WRITE(CARD,'('' Fixed Parameters:'')')
      CALL TVN(CARD)
      DO 554 J=1,NIX
      I=IX(J)
  554 WRITE(PAGE,132) I,PG(I),P(I)
      CALL PRPAGE(PAGE,1)
       END IF
C
      IF(IPR.GE.2)THEN
         WRITE(CARD,555)
  555 FORMAT(39(' *'))
         CALL TV(CARD)
      END IF
C
  580 DO 590 I=1,NPARMS
  590 P(I)=PC(I)
C     prevent changes to weights.
      IF(MODE.EQ.0) MODE=7
C            *************** END OF ITERATION ********************
      IF(IT.LT.MAXIT) GO TO 670
      WRITE(CARD,667)MAXIT
  667 FORMAT('DLSQ  DID NOT CONVERGE in',I3,' iterations.',16(' $'))
      CALL TV(CARD)
      GO TO 790
C
  670 CONTINUE
C     Recompute scale:
      IDF=NPTS-NEQS
      IF (IDF.LT.0) THEN
          CALL TV('You are trying to solve for more parameters than')
          CALL TVN('the number of available data!')
          CALL TV('QUITTING')
          CALL STSEPI
      ELSE IF (IDF.EQ.0) THEN
          CALL TV('No degrees of freedom left!')
          GO TO 790
      END IF
C
      IF(HUBER) THEN
C        check for excesssive rejected points (scale too small):
         IF ( (NREJ*5.GT.NPTS .AND. IT.GT.1) .OR.
     1          NPTS-NREJ.LT.NEQS) THEN 
                        IF (IPR.GT.0) CALL TVN('Restarted at 671')
            MODE=0
            RESTART=RESTART+1.
            SCALE=SCALE*3.
            OLDVAR=-1.
            DAMP=1.
            DO 671 I=1,NPARMS
  671       P(I)=PLOW(I)
            GO TO 100
         END IF
         IF (MODE.GT.1) THEN
C           see if we can use a redescending M-estimator:
            HUBER=.FALSE.
            DO 672 I = 1,NPARMS
               IF(ABS(DP(I)).GT.ABS(P(I)/128.)) HUBER=.TRUE.
  672       CONTINUE
         END IF
         IF (.NOT.HUBER) THEN
C           may be ready to switch.
            CALL SORT1(S1,NUSED)
            IF (IDF.GT.NREJ) THEN
              SCLNEW=SQRT(((S1((NUSED+1)/2) + S1(NUSED/2+1))*NEQS)/IDF)
            ELSE
              SCLNEW=SCALE*1.2
            END IF
            IF (ABS(SCLNEW-SCALE).GT.SCALE/16. .AND. RESTART.LT.3.)THEN
C               not converged.  Keep Huber.
                HUBER=.TRUE.
                DAMP=DAMP+1.E-9
            ELSE
C               converged, so can switch.
                DAMP=1.E-4
                IF (IPR.EQ.1) CALL TV('Using BIWEIGHT')
            END IF
C           re-scale, in any case.
            IF (RESTART.LT.3.) SCALE=SCLNEW
                        IF (IPR.GT.0) CALL TVN('Restarted at 672')
            MODE=0
            OLDVAR=-1.
         END IF
      ELSE
C        using biweight.
         IF (SCLDEN.GT.1. .AND. NPTS-NREJ.GT.NEQS) THEN
            SCLNEW=SQRT(NUSED*SCLNUM/(SCLDEN*(SCLDEN-1.)))
         ELSE
            SCLNEW=SCALE*2.
         END IF
         IF (ABS(SCLNEW-SCALE).GT.SCALE*RTTEST)THEN
C           need more iterations.
            IF (MODE.EQ.4) THEN
               IF(SCLNEW.GT.0. .AND. ABS(SCLNEW-SCALE).GT.SCALE/16.)THEN
                  DAMP=DAMP+1.E-9
                  MODE=0
C                       IF (IPR.GT.0) CALL TVN('Restarted at 673')
                  OLDVAR=-1.
               END IF
C		   debug biweight:
C                  write(card,*)SCALE,SCLNEW,NUSED
C                  call tv(card)
               IF (SCLNEW.GT.SCALE/2.) THEN
C                 fudge to allow for effects of scale change:
                  ABSMIN=ABSMIN*1.00000001D0
                  IF (SCLNEW.GT.SCALE) THEN
C                    Scale is increasing; adopt full change.
                     SCALE=SCLNEW
                     ABSMIN=ABSMIN*1.00000001D0
                  ELSE
C                    Scale is decreasing; adopt half change.
                     SCALE=(SCLNEW + SCALE)*0.5
                  END IF
               ELSE IF(SCLNEW.GT.0.) THEN
c                 Allow cautious decrease.
                  SCALE=SCALE*0.9
               END IF
            END IF
            GO TO 100
         END IF
      END IF
C
C     HUBER is now false if we can use biweight.
                IF (IPR.GT.1) THEN
                        IF(HUBER) THEN
                           CALL TV('HUBER')
                        ELSE
                           CALL TV('BIWEIGHT')
                        END IF
                END IF
      IF (HUBER) GO TO 100
C
C  TEST PARAMETERS for CONVERGENCE.  *******  (AT END OF ITERATION.)
      DO 775 I = 1,NPARMS
C                  debugging:
C                  write(card,*)I,DP(I),P(I)
C        IF(ABS(DP(I)).GT.ABS(TEST*P(I))) call tv(card)
         IF(ABS(DP(I)).GT.ABS(TEST*P(I))) GO TO 100
C  skip to 100 for another iteration if TEST is not met.
  775 CONTINUE
C  here if TEST is met by all non-zero P's.
      IF(DAMP.GT.RTTEST .OR. MODE.LT.4) GO TO 100
C  do not quit if overdamped.
C                            ******  FINAL RESULTS NOW  *********
  790 MODE=5
C  use lowest values for final results.
      DO 791 I=1,NPARMS
      P(I)=PLOW(I)
  791 PC(I)=PLOW(I)
C  do inverse of undamped matrix to find SP's.
      GO TO 101
C  get normal eqns.for final values.  return here from 135.
  800 DO 801 I=1,NPARMS
      K=I*(I+1)/2
C       Diagonal element: ignore missing parms...
      IF (A(K).EQ.0.D0) A(K)=1.D0
  801 CONTINUE
      CALL PARTIT(A,NPARMS,SP,NBLOX,MBLOK,*880)
C  USE SP AS SCRATCH SPACE TO INVERT A.
C
C  ************* CLEAN UP FINAL LEFTOVERS WHEN PROBLEM IS SOLVED.****
C  get final residuals and covariances, then
C  compute variances of parameters.
  802 I=0
      DO 810 NPT=1,NPTS
      IF(W(NPT).EQ.0.) I=I+1
C  remove deleted points from degrees of freedom.
  810 CONTINUE
      IDF = NPTS - NEQS - I
      IF(NEQS.EQ.0 .OR. IPR.LT.2) GO TO 820
C  print inverse matrix if IPR < 2.
      WRITE(CARD,812)
  812 FORMAT(3X,'I',30X,'Inverse of A(I,J)')
      CALL TV(' ')
      CALL TV(CARD)
      CALL TV(' ')
      K=0
      DO 815 I=1,NPARMS
      IF (.NOT.USED(I)) GO TO 815
      WRITE(PAGE,132)I,(A(J),J=K+1,K+I)
      CALL PRPAGE(PAGE,(I+4)/5)
      CALL TV (' ')
  815 K=K+I
      CALL TV (' ')
      WRITE(CARD,818)IT
  818 FORMAT(I6,' Iterations')
      CALL TV(CARD)
C
  820 IF(IDF.GT.0) WVAR=VAR/FLOAT(IDF)
C    clear SP.
      DO 830 I=1,NPARMS
  830 SP(I)=0.
C  compute STD.ERR's of P'S.
      DO 850 J=1,NPARMS
      I=J
      IF(A(J*(J+1)/2).LE.0.) A(J*(J+1)/2)=1.E38
C   fix up bad results from singular matrices.
      IF (.NOT.USED(I)) THEN
          SP(I) =0.D0
      ELSE
          SP(I) =SQRT(A(J*(J+1)/2)*WVAR)
      END IF
  850 CONTINUE
      RETURN
C  NORMAL TERMINATION. *****
C
C   singular matrix on cleanup:
  880 WRITE(CARD,496) 0.,IT
      CALL TV(CARD)
      IF (MODE.EQ.6) THEN
        DAMP=DAMP*4.
      ELSE
        DAMP=1./2.**48
        MODE=6
      END IF
      WRITE(CARD,901) DAMP
 901  FORMAT('Final matrix inversion damped by ',
     2       1PE9.2,' for stability')
C
C       ...but we have to regenerate the A matrix now.
      CALL TVN(CARD)
      GO TO 101
C
  905 WRITE(CARD,906) NEQS,NIX
  906 FORMAT('DLSQ FAILED.  No. of parameters VARIED =',I3,6X,
     1'No. held fixed =',I3)
      CALL TV(CARD)
      call SPACE2
      call tv('     checking derivatives...')
      CALL SPACE
C                     ******* CHECK DERIVATIVES NUMERICALLY. ******
C  DO ONCE ONLY.
      IF(PARTSW) GO TO 999
C     PARTSW=.TRUE.
      DO 910 I=1,NPARMS
C       denominator is 16777216.d0 for double-precision P's.
        DP(I)=P(I)/16777216.D0
        IF(DP(I).EQ.0.) DP(I)=1./16777216.D0
  910 CONTINUE
      TERM=0.
      VAR=0.
      DO 930 NPT=1,NPTS
      DO 915 J=1,NARGS
  915 Z(J)=X(J,NPT)
      DO 916 J=1,NKARGS
  916 KZ(J)=KX(J,NPT)
      NONEQ=0
C   COMPUTE PART.DERIV.'S IN YP.
      CALL YP
      YC(NPT)=YT
      NONEQ=1
      DO 920 J=1,NPARMS
      I=J
      PC(I)=P(I)
C     skip unused P's:
      IF (.NOT.USED(J)) GO TO 920
      P(I)=P(I)+DP(I)
      CALL YP
      DY=YT-YC(NPT)
      DY=DY/DP(I)
      IF(DY.EQ.0.) GO TO 920
C    SKIP ZERO TERMS.
      TERM=TERM+DY*DY
      VAR=VAR+(PART(I)-DY)**2
  918 IF(IPR.GT.1) THEN
         WRITE(CARD,919)NPT,I,PART(I),DY
  919    FORMAT(' Point no.',I5,5X,'PART(',I3,') =',1P,E15.7,5X,
     1   'Numer.Est.=',E15.7)
         IF (ABS((PART(I)-DY)/DY).GT.0.1)THEN
             CARD(79:)='*'
             IF (IPR.GT.1) CALL TVN(CARD)
         ELSE
             CARD(79:)=' '
             IF (IPR.GT.2) CALL TVN(CARD)
         END IF
      END IF
  920 P(I)=PC(I)
  930 CONTINUE
      IF (TERM.GT.0.) THEN
         VAR=SQRT(VAR/TERM)
         WRITE(CARD,935)VAR
  935    FORMAT('R.M.S. relative error of Partial Derivatives =',
     1    1PE10.2)
      ELSE
         CARD='Zero divisor; numerical test failed.'
      END IF
      CALL TV(CARD)
      IF(VAR.GT.0.1)THEN
        WRITE(CARD,'(''BAD DERIVATIVES'')')
        CALL TV(CARD)
      ELSE IF(VAR.LT.1.E-5) THEN
        WRITE(PAGE,990)
  990 FORMAT(/'Derivatives seem to be satisfactory.  Trouble may be'//
     110X,'A) incompatible model and data, due either to bad data'/
     215X,'or inappropriate model.'//10X,'B) found wrong minimum.'/)
        CALL TV(PAGE(2))
        CALL TV(PAGE(4))
        CALL TVN(PAGE(5))
        CALL TV(PAGE(7))
      ELSE
      END IF
      CALL TV('   Solution abandoned.')
  999 RETURN
      END
      SUBROUTINE PRPAGE(PAGE,N)
C
C       prints page loaded by fmt. 132 in dlsq.
C
        IMPLICIT INTEGER (I-N)
C
        INCLUDE 'MID_REL_INCL:kingfit.inc'
        INCLUDE 'MID_REL_INCL:dlsq.inc'
        CHARACTER*79 PAGE(MXP)
C
      DO 1320 I=1,N
 1320 CALL TVN(PAGE(I))
      RETURN
      END
      SUBROUTINE FIXP(L,VALUE,NIX)
C
C       Fixes P(L) at VALUE, whether already fixed or not.
C
      IMPLICIT NONE
C
C       Arguments:
C
      INTEGER L,NIX
      REAL VALUE
C
      INTEGER II
C
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
C
C
      PG(L)=VALUE
C
      DO 671 II=1,NIX
C           make sure it is not already fixed:
            IF(IX(II).EQ.L) RETURN
  671       CONTINUE
C     it is NOT already fixed.
      NIX=NIX+1
      IX(NIX)=L
C
      RETURN
C
      END
      SUBROUTINE UNFIXP(L,NIX)
C
C       Varies P(L), whether already fixed or not.
C
      IMPLICIT NONE
C
C       Arguments:
C
      INTEGER L,NIX
C
      INTEGER II
C
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
C
C
      DO 671 II=1,NIX
C           see if it is already fixed:
            IF(IX(II).EQ.L) GO TO 675
  671       CONTINUE
C     It is NOT already fixed.  We are OK.
      RETURN
C
C     It IS already fixed.  Remove it from NIX list:
  675 CONTINUE
      DO 680 II=II,NIX-1
C        (Just move rest of list down 1 slot.)
         IX(II)=IX(II+1)
  680 CONTINUE
      NIX=NIX-1
C
      RETURN
C
      END
      SUBROUTINE QFIX(L,VALUE,NIX)
C
C       Asks if user wants to hold P(L) fixed at VALUE.
C
      IMPLICIT NONE
C
C       Arguments:
C
      INTEGER L,NIX
      REAL VALUE
C
C
      INCLUDE 'MID_REL_INCL:kingfit.inc'
      INCLUDE 'MID_REL_INCL:dlsq.inc'
C
      CHARACTER*79 PAGE(21)
      COMMON /SCREEN/PAGE
C
      CHARACTER A1
C
C
   10 IF (ABS(VALUE).GT.0.01) THEN
        WRITE(PAGE(1),'(A,F6.3)')'Want to hold this fixed at ',VALUE
      ELSE
        WRITE(PAGE(1),'(A,1PE9.2)')'Want to hold this fixed at ',VALUE
      END IF
      CALL TV(PAGE(1))
      CALL ASKN('?',A1)
C
      IF (A1.EQ.'Y' .OR. A1.EQ.'O') THEN
         CALL FIXP(L,VALUE,NIX)
      ELSE IF (A1.EQ.'H') THEN
         CALL TV('The value suggested is a reasonable one.  Unless you')
         CALL TVN('know of a better value, reply YES.')
         GO TO 10
      ELSE
         IF (P(L).GT.0.01) THEN
           WRITE(PAGE(1),'(A,F6.3,A)')'Is',P(L),' acceptable?'
         ELSE
           WRITE(PAGE(1),'(A,1PE9.2,A)')'Is',P(L),' acceptable?'
         END IF
   20    CALL ASK(PAGE(1),A1)
         IF (A1.EQ.'Y')THEN
C            OK.
         ELSE IF (A1.EQ.'H') THEN
             CALL TV('That''s the current value; PEPSYS thinks it is a')
             WRITE (PAGE(1),'(A,G9.3)')'bad one.  A reasonable value is'
     1             ,VALUE
             CALL TVN(PAGE(1))
             GO TO 20
         ELSE
C            ask for advice.
             CALL QF('Please enter a reasonable value:',PG(L))
             CALL FIXP(L,PG(L),NIX)
         END IF
      END IF
C
      RETURN
      END
      SUBROUTINE IPRSET(IPR)
C
C       Sets IPR.
C
      INTEGER IPR
C
      REAL ARG
C
C
      CALL SPACE2
      CALL TV('Please specify how much ITERATION OUTPUT you want:')
      CALL SPACE
      CALL TV('     1:  no information about iterations')
      CALL TV('     2:  only iteration number and variance')
      CALL TV('     3:  values of parameters at each iteration')
      CALL TV('     4:  additional details')
      CALL TV('     5:  everything possible')
      CALL SPACE
      CALL QF('Please enter the NUMBER of your choice:',ARG)
      IPR=ARG-0.75
C
C
      RETURN
      END