Attribute VB_Name = "Module4"
'******************************************************************************
'*              RESOLUTION OF A SET OF NON-LINEAR EQUATIONS                   *
'* -------------------------------------------------------------------------- *
'* COPYRIGHT 1991, BY R.S. BAIN                                               *
'* Originally wrapped by bain@luther.che.wisc.edu on Tue Mar 30 16:27:48 1993 *
'*                                                                            *
'* As of this writing, we do not know how to contact Rod Bain, the            *
'* author of NNES.                                                            *
'*                                                                            *
'* -- David M. Gay (dmg@bell-labs.com)                                        *
'*    Bell Labs, Murray Hill                                                  *
'*    8 February 1999                                                         *
'*                                                                            *
'*                            Visual Basic 4.0 Release By J-P Moreau, Paris.  *
'*                                            (PART 2/3)                      *
'******************************************************************************
DefDbl A-H, O-Z
DefInt I-N


Sub dogleg(frstdg As Integer, newtkn, overch As Integer, overfl As Integer, maxexp, n, notrst, nunit, _
            output As Integer, beta, caulen, delta, etafac, newlen As Double, sqrtz, delf(), s(), _
            scalex(), sn(), ssdhat(), vhat())
'--------------------------------------------------------------------------
'    FEB. 23, 1992
'
'    THIS SUBROUTINE FINDS A TRUST REGION STEP USING THE
'    (DOUBLE) DOGLEG METHOD.
'--------------------------------------------------------------------------
Dim lambda As Double
Dim ss As String

  overfl = 0
  eta = one

  If output > 3 Then
    Line0
    ss = "    TRUST REGION STEP: " + Str$(notrst) + "  TRUST REGION LENGTH, DELTA: " + Str$(delta): Msg ss
    Line0
    ss = "    LENGTH OF NEWTON STEP, NEWLEN: " + Str$(newlen): Msg ss
  End If

'   CHECK FOR NEWTON STEP WITHIN TRUST REGION - IF SO USE NEWTON STEP

  If newlen <= delta Then
    For i = 1 To n
      s(i) = sn(i)
    Next i
    newtkn = 1
    temp = delta
    delta = newlen
    If output > 3 Then
      Line0
      Msg "       NEWTON STEP WITHIN ACCEPTABLE RANGE ( <= THAN DELTA)"
      If temp = delta Then
        ss = "       DELTA STAYS AT LENGTH OF NEWTON STEP: " + Str$(delta): Msg ss
      Else
        Msg "       DELTA SET TO LENGTH OF NEWTON STEP."
      End If
      Line0
      Msg "       FULL NEWTON STEP ATTEMPTED."
    End If
    Exit Sub
  Else

'   NEWTON STEP NOT WITHIN TRUST REGION - APPLY ^[DOUBLE)
'   DOGLEG PROCEDURE.  ^[IF ETAFAC EQUALS 1.0 THEN THE SINGLE
'   DOGLEG PROCEDURE IS BEING APPLIED).

    newtkn = 0
    If frstdg <> 0 Then

'   SPECIAL SECTION FOR FIRST DOGLEG STEP - CALCULATES
'   CAUCHY POINT (MINIMIZER OF MODEL FUNCTION IN STEEPEST
'   DESCENT DIRECTION OF THE OBJECTIVE FUNCTION).

      frstdg = 0
      If output > 4 Then
        Line0
        If etafac = one Then
          Msg "       FIRST SINGLE DOGLEG STEP."
        Else
          Msg "       FIRST DOUBLE DOGLEG STEP."
        End If
        Line0
        Msg "          SCALED CAUCHY STEP."
        Line0
      End If

'   NOTE: BETA AND SQRTZ WERE CALCULATED IN SUBROUTINE DELCAU

      zeta = sqrtz * sqrtz

'   FIND STEP TO CAUCHY POINT

      FACTOR = -(zeta / beta)
      For i = 1 To n
        ssdhat(i) = FACTOR * (delf(i) / scalex(i))
        If output > 4 Then
          ss = "           SSDHAT(" + Str$(i) + ") = " + Str$(ssdhat(i)): Msg ss
        End If
      Next i
      innerp overch, overfl, maxexp, n, n, n, nunit, output, delfts, delf, sn
      overfl = 0

'   PROTECT AGAINST (RARE) CASE WHEN CALCULATED DIRECTIONAL
'   DERIVATIVE EQUALS ZERO.

      If delfts <> zero Then
'   STANDARD EXECUTION
        gamma = (zeta / Abs(delfts)) * (zeta / beta)
        eta = etafac + (one - etafac) * gamma
      Else
        If output > 1 And wrnsup = 0 Then
          Line0
          Msg "    WARNING: DELFTS=0,  ETA SET TO 1.0 TO AVOID DIVISION BY ZERO."
        End If
        eta = one
      End If
      If output > 4 Then
        Line0
        ss = "          ETA = " + Str$(eta): Msg ss
        Line0
        Msg "        VHAT VECTOR     VHAT(I) = ETA*SN(I)*SCALEX(I) - SSDHAT(I)"
        Line0
      End If
      For i = 1 To n
        vhat(i) = eta * scalex(i) * sn(i) - ssdhat(i)
        If output > 4 Then
          ss = "            VHAT(" + Str$(i) + ") = " + Str$(vhat(i)): Msg ss
        End If
      Next i
    End If 'if frstdg
  End If 'if newlen

' ETA*NEWLEN <= DELTA MEANS TAKE STEP IN NEWTON DIRECTION TO TRUST REGION BOUNDARY.

  If eta * newlen <= delta Then
    If output > 4 Then
      Line0
      Msg "          ETA*NEWLEN <= DELTA     S(I) = (DELTA/NEWLEN)*SN(I)"
    End If
    For i = 1 To n
      s(i) = (delta / newlen) * sn(i)
    Next i
  Else

'   DISTANCE TO CAUCHY POINT >= DELTA MEANS TAKE STEP IN
'   STEEPEST DESCENT DIRECTION TO TRUST REGION BOUNDARY.
    If caulen >= delta Then
      If output > 4 Then
        Line0
        Msg "          CAULEN >= DELTA,  S(I)=(DELTA/CAULEN)*(SSDHAT(I)/SCALEX(I))"
      End If
      For i = 1 To n
        s(i) = (delta / caulen) * (ssdhat(i) / scalex(i))
      Next i
    Else
'     TAKE (DOUBLE) DOGLEG STEP

      innerp overch, overfl, maxexp, n, n, n, nunit, output, temp, ssdhat, vhat
      innerp overch, overfl, maxexp, n, n, n, nunit, output, tempv, vhat, vhat
      overfl = 0
      lambda = (-temp + Sqr(temp * temp - tempv * (caulen * caulen - delta * delta))) / tempv
      If output > 4 Then
        Line0
        Msg "          S(I) = (SSDHAT(I)+LAMBDA*VHAT(I))/SCALEX(I)"
        Line0
        ss = "          WHERE LAMBDA = " + Str$(lambda): Msg ss
      End If
      For i = 1 To n
        s(i) = (ssdhat(i) + lambda * vhat(i)) / scalex(i)
      Next i
    End If
  End If
  If output > 3 Then
    Line0
    Line0
    Msg "          REVISED STEP FROM SUBROUTINE DOGLEG:"
    Line0
    For i = 1 To n
      ss = "           S(" + Str$(i) + ") = " + Str$(s(i)): Msg ss
    Next i
  End If
End Sub

Sub gradf(overch As Integer, overfl As Integer, restrt As Integer, sclfch As Integer, _
           sclxch As Integer, jupdm, maxexp, n, nunit, output As Integer, qnupdm As Integer, _
           delf(), fvecc(), jac() As Double, scalef(), scalex())
'----------------------------------------------------------------------------
'    FEB. 23, 1992
'
'    THIS SUBROUTINE COMPUTES THE GRADIENT OF THE FUNCTION
'
'           F=1/2{SCALEF*FVECC)^(SCALEF*FVECC)
'
'    WHICH IS USED AS THE OBJECTIVE FUNCTION FOR MINIMIZATION.
'
'    NOTE: WHEN THE FACTORED FORM OF THE JACOBIAN IS UPDATED IN QUASI-NEWTON
'          METHODS THE GRADIENT IS UPDATED AS WELL IN THE SAME SUBROUTINE -
'          IT IS PRINTED HERE THOUGH.  IN THESE CASES QNUPDM > 0.
'----------------------------------------------------------------------------
Dim wv1(ISIZE)
Dim s As String

  If restrt <> 0 Or jupdm = 0 Or qnupdm = 0 Then

'   GRADIENT NOT ALREADY UPDATED:    FIND DELF = J^F

    For i = 1 To n
      If sclfch <> 0 Then
        wv1(i) = fvecc(i) * scalef(i) * scalef(i)
      Else
        wv1(i) = fvecc(i)
      End If
    Next i

    If (overch) Then
'   check each entry individually
      atvov overfl, maxexp, n, nunit, output, jac, wv1, delf
    Else
      mmulvt n, jac, wv1, delf  'modified 04/10/2007
    End If

'   PRINT GRADIENT VECTOR, DELF

    If output > 3 Then
      If sclxch = 0 Then
        Line0
        If (Not sclfch) Then
          Msg "    GRADIENT OF OBJECTIVE FUNCTION:"
        Else
          Msg "    GRADIENT OF SCALED OBJECTIVE FUNCTION:"
        End If
        Line0
        For i = 1 To n
          s = "    DELF(" + Str$(i) + ") = " + Str$(delf(i)): Msg s
        Next i
      Else
        Line0
        Msg "    GRADIENT OF OBJECTIVE FUNCTION IN SCALED X UNITS:"
        Line0
        For i = 1 To n
          s = "    DELF(" + Str$(i) + ") = " + Str$(delf(i)) + "         DELF(" + Str$(i) + ") = " + _
              Str$(scalef(i) * scalef(i) * delf(i) / scalex(i)): Msg s
        Next i
      End If
    End If
  End If

End Sub

Sub initch(instop, linesr, newton, overfl As Integer, sclfch As Integer, sclxch As Integer, _
            acptcr As Integer, contyp As Integer, jactyp, jupdm, maxexp, n, nunit, output As Integer, _
            qnupdm As Integer, stopcr As Integer, trupdm As Integer, epsmch, fcnold, ftol, _
            boundl(), boundu(), fvecc(), scalef(), scalex(), xc())
'-----------------------------------------------------------------------------
'    AUG. 27, 1991
'
'    THIS SUBROUTINE FIRST CHECKS TO SEE IF N IS WITHIN THE ACCEPTABLE RANGE.
'
'    THE SECOND CHECK IS TO SEE IF THE INITIAL ESTIMATE IS
'    ALREADY A SOLUTION BY THE FUNCTION VALUE CRITERION, FTOL.
'
'    THE THIRD CHECK IS MADE TO SEE IF THE NEWTON OPTION IS BEING
'    USED WITH THE LINE SEARCH.  IF NOT A WARNING IS GIVEN AND
'    THE LINE SEARCH OPTION IS INVOKED.
'
'    THE FOURTH CHECK IS TO ENSURE APPLICABILITY OF SELECTED
'    VALUES FOR INTEGER CONSTANTS.
'
'    THE FIFTH CHECK IS TO WARN THE USER IF INITIAL ESTIMATES ARE NOT WITHIN
'    THE RANGES SET BY THE BOUNDL AND BOUNDU VECTORS.
'    CONTYP IS CHANGED FROM 0 TO 1 IF ANY BOUND HAS BEEN SET BY THE USER
'
'    THE SIXTH CHECK ENSURES BOUNDL^[I) < BOUNDU^[I) FOR ALL I.
'----------------------------------------------------------------------------
'Label: 130
Dim frster As Integer
Dim s As String

  instop = 0
  temp1 = -(ten ^ maxexp)
  temp2 = ten ^ maxexp

' CHECK FOR N IN RANGE

  If n <= 0 Then
    instop = 1
    Line1
    Line0
    Msg "  N IS OUT OF RANGE - RESET TO POSITIVE INTEGER."
    Line0
    Line1
  End If

' CHECK FOR SCALING FACTORS POSITIVE

  frster = 1
  sclfch = 0
  sclxch = 0
  For i = 1 To n
    If scalef(i) <= zero Then
      If frster <> 0 Then
        instop = 1
        frster = 0
        Line1
      End If
      Line0
      s = "       SCALEF(" + Str$(i) + ") = " + Str$(scalef(i)) + "    SHOULD BE POSITIVE.": Msg s
    End If
    If scalef(i) <> one Then sclfch = 1
    If scalex(i) <= zero Then
      If frster <> 0 Then
        instop = 1
        frster = 0
        Line1
      End If
      Line0
      s = "       SCALEX(" + Str$(i) + ") = " + Str$(scalex(i)) + "    SHOULD BE POSITIVE.": Msg s
    End If
    If scalex(i) <> one Then sclxch = 1
  Next i
  If frster = 0 Then
    Line0
    Line1
  End If

' EVALUATE INITIAL RESIDUAL VECTOR AND OBJECTIVE FUNCTION AND
' CHECK TO SEE IF THE INITIAL GUESS IS ALREADY A SOLUTION.

  fcn1 overfl, n, fvecc, xc

' NOTE: NUMBER OF LINE SEARCH FUNCTION EVALUATIONS, NFUNC,
' INITIALIZED AT 1 WHICH REPRESENTS THIS EVALUATION.

  If overfl <> 0 Then
    Line1
    Line0
    Msg "       OVERFLOW IN INITIAL FUNCTION VECTOR EVALUATION."
    Line0
    Line1
    instop = 1
    Exit Sub
  End If

  fcnevl overfl, maxexp, n, nunit, output, epsmch, fcnold, fvecc, scalef

  If overfl <> 0 Then
    Line1
    Line0
    Msg "       OVERFLOW IN INITIAL OBJECTIVE FUNCTION EVALUATION."
    Line0
    Line1
    instop = 1
    Exit Sub
  End If

' CHECK FOR SOLUTION USING SECOND STOPPING CRITERION

  For i = 1 To n
    If Abs(fvecc(i)) > ftol Then GoTo 130
  Next i

  instop = 1
  Line1
  Line0
  Msg "  WARNING: THIS IS ALREADY A SOLUTION BY THE CRITERIA OF THE SOLVER"
  Line0

'    IF THE PROBLEM IS BADLY SCALED THE OBJECTIVE FUNCTION MAY MEET THE
'    TOLERANCE ALTHOUGH THE INITIAL ESTIMATE IS NOT THE SOLUTION.

  Msg "  THIS MAY POSSIBLY BE ALLEVIATED BY RESCALING THE PROBLEM IF THE"
  Msg "  INITIAL ESTIMATE IS KNOWN NOT TO BE A SOLUTION."
  Line0
  Line1

' CHECK FOR NEWTON'S METHOD REQUESTED BUT LINE SEARCH NOT BEING USED

130 If newton <> 0 And linesr = 0 Then
    linesr = 1
    Line1
    Line0
    Msg "  WARNING: INCOMPATIBLE OPTIONS: NEWTON=TRUE AND LINESR=FALSE."
    Msg "  LINESR SET TO TRUE; EXECUTION OF NEWTON METHOD CONTINUING..."
    Line0
    Line1
  End If

' CHECK INTEGER CONSTANTS

  If acptcr <> 1 And acptcr <> 12 Then
    Line1
    Line0
    s = "  ACPTCR NOT AN ACCEPTABLE VALUE: " + Str$(acptcr): Msg s
    instop = 1
    Line0
    Line1
  End If
  If jactyp < 0 Or jactyp > 3 Then
    Line1
    Line0
    s = "  JACTYP: " + Str$(jactyp) + " - NOT IN PROPER RANGE.": Msg s
    instop = 1
    Line0
    Line1
  End If
  If stopcr <> 1 And stopcr <> 12 And (stopcr <> 2 And stopcr <> 3) Then
    Line1
    Line0
    s = "  STOPCR NOT AN ACCEPTABLE VALUE: " + Str$(stopcr): Msg s
    instop = 1
    Line0
    Line1
  End If
  If qnupdm < 0 Or qnupdm > 1 Then
    Line1
    Line0
    s = "  QNUPDM: " + Str$(qnupdm) + " - NOT IN PROPER RANGE.": Msg s
    instop = 1
    Line0
    Line1
  End If
  If trupdm < 0 Or trupdm > 1 Then
    Line1
    Line0
    s = "  TRUPDM: " + Str$(trupdm) + " - NOT IN PROPER RANGE.": Msg s
    instop = 1
    Line0
    Line1
  End If
  If jupdm < 0 Or jupdm > 2 Then
    Line1
    Line0
    s = "  JUPDM: " + Str$(jupdm) + " - NOT IN PROPER RANGE.": Msg s
    instop = 1
    Line0
    Line1
  End If

'  CHECK FOR INITIAL ESTIMATES NOT WITHIN SPECIFIED BOUNDS AND
'  SET CONTYP TO 1 => AT LEAST ONE BOUND IS IN EFFECT.

  contyp = 0
  For i = 1 To n
    If boundl(i) <> temp1 Or boundu(i) <> temp2 Then
      contyp = 1
      Exit Sub
    End If
  Next i

  frster = 1
  If contyp <> 0 Then
    For i = 1 To n

' CHECK FOR INITIAL ESTIMATES OUT OF RANGE AND LOWER
' BOUND GREATER THAN OR EQUAL TO THE UPPER BOUND.

      If xc(i) < boundl(i) Or xc(i) > boundu(i) Then
        If (frster) Then
          instop = 1
          frster = 0
          Line1
          Line0
          Msg "       COMPONENTS MUST BE WITHIN BOUNDS:"
          Line0
          Msg "      NO         XC               BOUNDL          BOUNDU"
          Line0
        End If
        s = "     " + Str$(i) + "   " + Str$(xc(i)) + "      " + Str$(boundl(i)) + "    " + Str$(boundu(i)): Msg s
      End If
    Next i

    If frster = 0 Then
      Line0
      Line1
    End If

    frster = 1
    For i = 1 To n
      If boundl(i) >= boundu(i) Then
        If frster <> 0 Then
          frster = 0
          Line1
          Line0
          Msg "       LOWER BOUND MUST BE LESS THAN UPPER BOUND - VIOLATIONS LISTED."
          Line0
        End If
        s = "       BOUNDL(" + Str$(i) + ") = " + Str$(boundl(i)) + "    BOUNDU(" + Str$(i) + ") = " + Str$(boundu(i)): Msg s
      End If
    Next i

    If frster = 0 Then
      Line0
      Line1
    End If
  
  End If

End Sub


Sub jaccd(n, nunit, output As Integer, epsmch, fvecj1(), fvecj2(), _
           jacfdm() As Double, scalex(), xc())
'----------------------------------------------------------------------
'    FEB. 11, 1991
'
'    THIS SUBROUTINE EVALUATES THE JACOBIAN USING CENTRAL DIFFERENCES.
'
'    FVECJ1 AND FVECJ2 ARE TEMPORARY VECTORS TO HOLD THE
'    RESIDUAL VECTORS FOR THE CENTRAL DIFFERENCE CALCULATION.
'----------------------------------------------------------------------
Dim overfl As Integer

  overfl = 0
  curtep = epsmch ^ 0.33

  For j = 1 To n
    deltaj = curtep * Sign((XMAX(Abs(xc(j)), one / scalex(j))), xc(j))
    tempj = xc(j)
    xc(j) = xc(j) + deltaj

'   NOTE: THIS STEP IS FOR FLOATING POINT ACCURACY ONLY

    deltaj = xc(j) - tempj

    fcn1 overfl, n, fvecj1, xc
    xc(j) = tempj - deltaj
    fcn1 overfl, n, fvecj2, xc
    If (overfl <> 0 And output > 2) And wrnsup = 0 Then
      Line0
      Msg "    WARNING: OVERFLOW IN FUNCTION VECTOR IN <JACCD>."
    End If
    For k = 1 To n
      jacfdm((k - 1) * ISIZE + j) = (fvecj1(k) - fvecj2(k)) / (two * deltaj)
    Next k
    xc(j) = tempj
  Next j
  
End Sub

Sub jacfd(jactyp, n, nunit, output As Integer, epsmch, boundl(), boundu(), _
            fvecc(), fvecj1(), jacfdm() As Double, scalex(), xc())
'-----------------------------------------------------------------------------
'   FEB. 15, 1991
'
'   THIS SUBROUTINE EVALUATES THE JACOBIAN USING ONE-SIDED FINITE DIFFERENCES.
'
'   JACTYP "1" SIGNIFIES FORWARD DIFFERENCES
'   JACTYP "2" SIGNIFIES BACKWARD DIFFERENCES
'
'   FVECJ1 IS A TEMPORARY VECTOR WHICH STORES THE RESIDUAL
'   VECTOR FOR THE FINITE DIFFERENCE CALCULATION.
'-----------------------------------------------------------------------------
Dim overfl As Integer
Dim s As String

  sqrtep = Sqr(epsmch)

'   FINITE-DIFFERENCE CALCULATION BY COLUMNS

  For j = 1 To n

'   DELTAJ IS THE STEP SIZE - IT IS ALWAYS POSITIVE

    deltaj = sqrtep * XMAX(Abs(xc(j)), one / scalex(j))
    tempj = xc(j)
'   temporary storage of xc(j)
    If jactyp = 1 Then
      If xc(j) + deltaj <= boundu(j) Then
'     STEP WITHIN BOUNDS - COMPLETE FORWARD DIFFERENCE
        xc(j) = xc(j) + deltaj
        deltaj = xc(j) - tempj
        fordif overfl, j, n, deltaj, fvecc, fvecj1, jacfdm, xc
      Else
'     STEP WOULD VIOLATE BOUNDU - TRY BACKWARD DIFFERENCE
        If xc(j) - deltaj >= boundl(j) Then
          xc(j) = xc(j) - deltaj
          bakdif overfl, j, n, deltaj, tempj, fvecc, fvecj1, jacfdm, xc
        Else
'       STEP WOULD ALSO VIOLATE BOUNDL - IF THE DIFFERENCE IN THE
'       BOUNDS, (BOUNDU-BOUNDL), IS GREATER THAN DELTAJ CALCULATE
'       THE FUNCTION VECTOR AT EACH BOUND AND USE THIS DIFFERENCE -
'       THIS REQUIRES ONE EXTRA FUNCTION EVALUATION.
'       THE CURRENT FVECC IS STORED IN WV3, THEN REPLACED.
          If boundu(j) - boundl(j) >= deltaj Then
            bnddif overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc
            If output > 2 And wrnsup = 0 And overfl = 0 Then
              Line0
              Msg "    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES."
              s = "           LOWER AND UPPER BOUNDS USED FOR JACOBIAN COLUMN: " + Str$(j): Msg s
              Msg "           THIS REQUIRED ONE EXTRA FUNCTION EVALUATION."
            End If
          Else
'         BOUNDS ARE EXTREMELY CLOSE (BUT NOT EQUAL OR
'         THE PROGRAM WOULD HAVE STOPPED IN INITCH).
            If output > 2 And wrnsup = 0 And overfl = 0 Then
              Line0
              Msg "    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES."
              s = "           BOUNDS ARE EXTREMELY CLOSE FOR COMPONENT: " + Str$(j): Msg s
              Msg "           FINITE DIFFERENCE JACOBIAN IS UNRELIABLE."
            End If
            bnddif overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc
          End If
        End If
      End If
    Else   'jactyp<>1
      If xc(j) - deltaj >= boundl(j) Then
        xc(j) = xc(j) - deltaj
        bakdif overfl, j, n, deltaj, tempj, fvecc, fvecj1, jacfdm, xc
      Else
        If xc(j) + deltaj <= boundu(j) Then
          xc(j) = xc(j) + deltaj
          fordif overfl, j, n, deltaj, fvecc, fvecj1, jacfdm, xc
        Else
          If boundu(j) - boundl(j) >= deltaj Then
            bnddif overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc
            If output > 2 And wrnsup = 0 And overfl = 0 Then
              Line0
              Msg "    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES."
              s = "           LOWER AND UPPER BOUNDS USED FOR JACOBIAN COLUMN: " + Str$(j): Msg s
              Msg "           THIS REQUIRED ONE EXTRA FUNCTION EVALUATION."
            End If
          Else
            bnddif overfl, j, n, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, xc
            For k = 1 To n
              jacfdm((k - 1) * ISIZE + j) = zero
            Next k
            If output > 2 And wrnsup = 0 And overfl = 0 Then
              Line0
              Msg "    WARNING: BOUNDS TOO CLOSE FOR 1-SIDED FINITE-DIFFERENCES."
              s = "           BOUNDS ARE EXTREMELY CLOSE FOR COMPONENT: " + Str$(j): Msg s
              Msg "           FINITE DIFFERENCE JACOBIAN IS UNRELIABLE."
            End If
          End If
        End If
      End If
    End If
    If overfl <> 0 And output > 2 And wrnsup = 0 Then
      Line0
      Msg "    WARNING: OVERFLOW IN FUNCTION VECTOR IN SUBROUTINE JACFD."
      s = "           BOUNDS ARE EXTREMELY CLOSE FOR COMPONENT: " + Str$(j): Msg s
      For k = 1 To n
        jacfdm((k - 1) * ISIZE + j) = zero
      Next k
    End If
    xc(j) = tempj
  Next j
 
End Sub

Sub jacobi(checkj As Integer, jacerr, overfl As Integer, jactyp, n, nunit, output As Integer, _
            epsmch, fdtolj, boundl(), boundu(), fvecc(), fvecj1(), fvecj2(), _
            jac() As Double, jacfdm() As Double, scalex(), xc())
'-------------------------------------------------------------------------
'    APR. 13, 1991
'
'    THIS SUBROUTINE EVALUATES THE JACOBIAN.  IF CHECKJ IS TRUE
'    THEN THE ANALYTICAL JACOBIAN IS CHECKED NUMERICALLY.
'
'    jacob IS A USER-SUPPLIED ANALYTICAL JACOBIAN USED ONLY IF JACTYP=0.
'    THE JACOBIAN NAME MAY BE CHANGED BY USING THE EXTERNAL STATEMENT
'    IN THE MAIN DRIVER.
'
'    JACFD ESTIMATES THE JACOBIAN USING FINITE DIFFERENCES:
'    FORWARD IF JACTYP=1 OR BACKWARD IF JACTYP=2.
'
'    JACCD ESTIMATES THE JACOBIAN USING CENTRAL DIFFERENCES.
'
'    IF THE ANALYTICAL JACOBIAN IS CHECKED THE FINITE DIFFERENCE
'    JACOBIAN IS STORED IN "JACFDM" AND THEN COMPARED.
'
'    FRSTER  INDICATES FIRST ERROR - USED ONLY TO SET BORDERS FOR OUTPUT
'    JACERR  FLAG TO INDICATE TO THE CALLING PROGRAM AN ERROR
'            IN THE ANALYTICAL JACOBIAN
'-------------------------------------------------------------------------
Dim frster As Integer
Dim s As String

  frster = 1
  jacerr = 0
  overfl = 0

  If jactyp = 0 Then
    jacob overfl, n, jac, xc
  ElseIf jactyp = 1 Or jactyp = 2 Then
    jacfd jactyp, n, nunit, output, epsmch, boundl, boundu, fvecc, fvecj1, jac, scalex, xc
  Else
    jaccd n, nunit, output, epsmch, fvecj1, fvecj2, jac, scalex, xc
  End If

  If (jactyp = 0) And (checkj) Then

'   NOTE: JACTYP=0 SENT TO JACFD PRODUCES A FORWARD DIFFERENCE ESTIMATE OF THE JACOBIAN.

    jacfd jactyp, n, nunit, output, epsmch, boundl, boundu, fvecc, fvecj1, jacfdm, scalex, xc

    For j = 1 To n
      For i = 1 To n
        If Abs((jacfdm((i - 1) * ISIZE + j) - jac((i - 1) * ISIZE + j)) / XMAX(XMAX(Abs(jac((i - 1) * ISIZE + j)), _
           Abs(jacfdm((i - 1) * ISIZE + j))), one)) > fdtolj Then
          jacerr = 1
          If output >= 0 Then
            If frster <> 0 Then
              frster = 0
              Line1
            End If
            Line0
            s = "    CHECK JACOBIAN TERM (" + Str$(i) + "," + Str$(j) + "):": Msg s
            Line0
            s = "    ANALYTICAL DERIVATIVE IS " + Str$(jac((i - 1) * ISIZE + j)): Msg s
            s = "     NUMERICAL DERIVATIVE IS " + Str$(jacfdm((i - 1) * ISIZE + j)): Msg s
            Line0
            Line1
          End If
        End If
      Next i
    Next j
  End If

End Sub

Sub lineS(abort As Integer, absnew As Integer, deuflh As Integer, geoms As Integer, newton, _
          overch As Integer, overfl As Integer, qnfail As Integer, qrsing As Integer, _
          restrt As Integer, sclfch As Integer, sclxch As Integer, acpcod As Integer, _
          acptcr As Integer, contyp As Integer, isejac, itnum, jupdm, maxexp, maxlin, mgll, mnew, _
          n, narmij, nfunc, nunit, output As Integer, qnupdm As Integer, stopcr As Integer, _
          trmcod As Integer, alpha, confac, epsmch, fcnmax, fcnnew, fcnold, lam0 As Double, _
          maxstp As Double, newlen As Double, sbrnrm, sigma, a(), boundl(), boundu(), delf(), _
          ftrack(), fvec(), h(), hhpi(), jac() As Double, rdiag(), rhs(), s(), sbar(), scalef(), _
          scalex(), sn(), strack(), xc(), xplus())
'--------------------------------------------------------------
'    SEPT. 9, 1991
'--------------------------------------------------------------
'Label: 90
Dim lambda As Double
Dim mu As Double
Dim newmax As Double
Dim norm As Double
Dim nrmpre As Double
Dim convio As Integer
Dim wv2(ISIZE)
Dim ss As String

  convio = 0
  overfl = 0
  abort = 0
  qnfail = 0

  If newton <> 0 Or absnew <> 0 Then
    
    If (newton) Then

'   FIND NEXT ITERATE FOR PURE NEWTON'S METHOD
      For k = 1 To n
        xplus(k) = xc(k) + sn(k)
      Next k
    Else

'   FIND NEXT ITERATE FOR "ABSOLUTE" NEWTON'S METHOD.
'   IF COMPONENT I WOULD BE OUTSIDE ITS BOUND THEN TAKE ABSOLUTE VALUE OF
'   THE VIOLATION AND GO THIS DISTANCE INTO THE FEASIBLE REGION.
'   ENSURE THAT THIS REFLECTION OFF ONE BOUND DOES NOT VIOLATE THE OTHER

      For i = 1 To n
        wv2(i) = zero
        If sn(i) >= zero Then
          If xc(i) + sn(i) > boundu(i) Then
            convio = 1
            wv2(i) = two
            xplus(i) = XMAX(two * boundu(i) - xc(i) - sn(i), boundl(i))
          Else
            xplus(i) = xc(i) + sn(i)
          End If
        Else
          If xc(i) + sn(i) < boundl(i) Then
            convio = 1
            wv2(i) = -two
            xplus(i) = XMIN(two * boundl(i) - xc(i) - sn(i), boundu(i))
          Else
            xplus(i) = xc(i) + sn(i)
          End If
        End If
      Next i
    End If

    If convio <> 0 And output > 4 Then
      Line0
      Msg "       CONSTRAINT VIOLATORS IN ABSOLUTE NEWTON'S METHOD"
      Line0
      Msg "       COMPONENT  PROPOSED POINT  VIOLATED BOUND  FEASIBLE VALUE"
      Line0
      For i = 1 To n
        If wv2(i) > zero Then
          ss = "     " + Str$(i) + "     " + Str$(xc(i) + sn(i)) + "     " + Str$(boundu(i)) + "    " + Str$(xplus(i))
          Msg ss
        ElseIf wv2(i) < zero Then
          ss = "     " + Str$(i) + "     " + Str$(xc(i) + sn(i)) + "     " + Str$(boundl(i)) + "    " + Str$(xplus(i))
          Msg ss
        End If
      Next i
    End If

    fcn1 overfl, n, fvec, xplus

    nfunc = nfunc + 1
    If overfl <> 0 Then
      overfl = 0
      fcnnew = ten ^ maxexp
      If output > 2 Then
        Line0
        Msg "       POTENTIAL OVERFLOW DETECTED IN FUNCTION EVALUATION."
      End If
      GoTo 90
    End If

    fcnevl overfl, maxexp, n, nunit, output, epsmch, fcnnew, fvec, scalef

'    RETURN FROM PURE NEWTON'S METHOD - OTHERWISE CONDUCT LINE SEARCH

90  Exit Sub
  End If

  If output > 3 Then
    Line0
    Line0
    Msg "    SUMMARY OF LINE SEARCH"
    Line0
  End If

'    SHORTEN NEWTON STEP IF LENGTH IS GREATER THAN MAXSTP

  If newlen > maxstp Then
    If output > 3 Then
      Line0
      Msg "       LENGTH OF NEWTON STEP SHORTENED TO MAXSTP."
    End If
    For k = 1 To n
      sn(k) = sn(k) * maxstp / newlen
    Next k
  End If

'    CHECK DIRECTIONAL DERIVATIVE (MAGNITUDE AND SIGN).

  innerp overch, overfl, maxexp, n, n, n, nunit, output, delfts, delf, sn

  If overfl <> 0 Then
    overfl = 0
    If output > 2 And wrnsup = 0 Then
      Line0
      ss = "    WARNING: DIRECTIONAL DERIVATIVE, DELFTS, SET TO " + Str$(delfts): Msg ss
    End If
  End If

'   REVERSE SEARCH DIRECTION IF DIRECTIONAL DERIVATIVE IS POSITIVE

  If delfts > zero Then
    For k = 1 To n
      sn(k) = -sn(k)
    Next k
    delfts = -delfts
    If output > 2 And wrnsup = 0 Then
      Line0
      Msg "    WARNING: DIRECTIONAL DERIVATIVE IS POSITIVE: DIRECTION REVERSED."
    End If
  End If

'   OUTPUT INFORMATION

  If output > 3 Then
    ss = "       INNER PRODUCT OF DELF AND SN, DELFTS: ......." + Str$(delfts): Msg ss
    If sclxch = 0 Then
      ss = "       LENGTH OF NEWTON STEP, NEWLEN: .............." + Str$(newlen): Msg ss
    Else
      ss = "       LENGTH OF SCALED NEWTON STEP, NEWLEN: ......." + Str$(newlen): Msg ss
    End If
    ss = "       MAXIMUM STEP LENGTH ALLOWED, MAXSTP: ........" + Str$(Int(maxstp * 1000) / 1000): Msg ss
  End If

'    ESTABLISH INITIAL RELAXATION FACTOR

  If deuflh <> 0 Then

'    AT FIRST STEP IN DAMPED NEWTON OR AFTER EXPLICIT JACOBIAN EVALUATION
'    IN QUASI-NEWTON OR IF THE STEP SIZE IS WITHIN STOPPING TOLERANCE BUT
'    STOPCR=3, THE LINE SEARCH IS STARTED AT LAMBDA=1.

    If (isejac = 1 Or trmcod = 1) And stopcr = 3 Then
      lambda = lam0
    Else
      For k = 1 To n
        wv2(k) = (sbar(k) - sn(k)) * scalex(k)
      Next k

      twonrm overfl, maxexp, n, epsmch, norm, wv2

'     PREVENT DIVIDE BY ZERO IF NORM IS ZERO (UNDERFLOWS)

      If norm < epsmch Then
'     START LINE SEARCH AT LAMBDA=LAM0, USE DUMMY MU
        mu = ten
      Else
        mu = nrmpre * lambda / norm
      End If

      If output > 4 Then
        Line0
        ss = "       DEUFLHARD TEST RATIO, MU: " + Str$(mu): Msg ss
      End If

'     SET INITIAL LAMBDA DEPENDING ON MU.   THIS IS A MODIFICATION OF
'     DEUFLHARD 'S METHOD WHERE THE CUTOFF VALUE WOULD BE 0.7 FOR LAM0=1.0

      If mu > lam0 / ten Then
        lambda = lam0
      Else
        lambda = lam0 / ten
      End If
    End If
  Else
    lambda = lam0
  End If

'   STORE LENGTH OF NEWTON STEP.  IF NEWTON STEP LENGTH WAS
'   GREATER THAN MAXSTP IT WAS SHORTENED TO MAXSTP.

  nrmpre = XMIN(maxstp, newlen)

'   ESTABLISH FCNMAX AND NEWMAX FOR NONMONOTONIC LINE SEARCH

  newmax = newlen
  fcnmax = fcnold
  If isejac > narmij Then
    If isejac < narmij + mgll Then
      For j = 1 To mnew
        fcnmax = XMAX(fcnmax, ftrack(j - 1))
        newmax = XMAX(newmax, strack(j - 1))
      Next j
    Else
      For j = 0 To mnew
        fcnmax = XMAX(fcnmax, ftrack(j))
        newmax = XMAX(newmax, strack(j))
      Next j
    End If
  End If

  If output > 3 Then
    Line0
    Line0
    If sclxch = 0 Then
      Msg "       LINE SEARCH"
    Else
      Msg "       LINE SEARCH (X'S GIVEN IN UNSCALED UNITS)"
    End If
  End If

'   CONDUCT LINE SEARCH

  deufls abort, deuflh, geoms, overch, overfl, qnfail, qrsing, restrt, _
         sclfch, sclxch, acpcod, acptcr, contyp, itnum, jupdm, maxexp, _
         maxlin, n, nfunc, nunit, output, qnupdm, stopcr, alpha, confac, _
         delfts, epsmch, fcnmax, fcnnew, fcnold, lambda, newmax, sbrnrm, _
         sigma, a, h, boundl, boundu, delf, fvec, hhpi, jac, rdiag, rhs, _
         s, sbar, scalef, scalex, sn, xc, xplus

End Sub

Sub llfa(overch As Integer, overfl As Integer, sclfch As Integer, sclxch As Integer, isejac, _
          maxexp, n, nunit, output As Integer, epsmch, omega, a(), delf(), fvec(), fvecc(), _
          jac() As Double, plee(), rdiag(), s(), scalef(), scalex(), xc(), xplus())
'---------------------------------------------------------
'    FEB. 23, 1991
'
'    THE LEE AND LEE QUASI-NEWTON METHOD IS APPLIED TO
'    THE FACTORED FORM OF THE JACOBIAN.
'
'    NOTE: T AND W ARE TEMPORARY WORKING VECTORS ONLY.
'---------------------------------------------------------
Dim t(ISIZE), w(ISIZE), wv3(ISIZE)
Dim temp1(ISIZE * ISIZE), temp2(ISIZE * ISIZE)
Dim skipup As Integer
Dim ss As String

  overfl = 0
  sqrtep = Sqr(epsmch)
  skipup = 1

  For i = 1 To n
    a((i - 1) * ISIZE + i) = rdiag(i)
    s(i) = xplus(i) - xc(i)
  Next i

'  R IS IN THE UPPER TRIANGLE OF A.

'       T=RS

  uvmul n, n, n, n, a, s, t

'  FORM PART OF NUMERATOR AND CHECK TO SEE IF A SIGNIFICANT
'  CHANGE WOULD BE MADE TO THE JACOBIAN.

  For i = 1 To n
    For k = 1 To n
      w(k) = jac(k)
    Next k
    innerp overch, overfl, maxexp, n, n, n, nunit, output, Sum, w, t
    w(i) = scalef(i) * (fvec(i) - fvecc(i)) - Sum

'  TEST TO ENSURE VECTOR W IS NONZERO.  IF W^[I]=0 FOR
'  ALL I THEN THE UPDATE IS SKIPPED - SKIPUP IS TRUE.

    If Abs(w(i)) > sqrtep * scalef(i) * (Abs(fvec(i)) + Abs(fvecc(i))) Then
      skipup = 0
'  update to be performed
    Else
      w(i) = zero
    End If
  Next i

  If skipup = 0 Then

'   T=Q^W   Q^ IS STORED IN JAC

    avmul n, n, n, n, jac, w, t

'   FIND DENOMINATOR; FORM W=S^P (P IS SYMMETRIC SO PS IS FOUND)

    avmul n, n, n, n, plee, s, w
    If sclxch <> 0 Then
'   SCALE W TO FIND DENOMINATOR
      For k = 1 To n
        wv3(k) = w(k) * scalex(k) * scalex(k)
      Next k
    Else
      For k = 1 To n
        temp1((k - 1) * ISIZE + 1) = w(k)
        temp2((k - 1) * ISIZE + 1) = wv3(k)
      Next k
      matcop n, n, 1, 1, n, 1, temp1, temp2
      For k = 1 To n
        wv3(k) = temp2((k - 1) * ISIZE + 1)
      Next k
    End If

    innerp overch, overfl, maxexp, n, n, n, nunit, output, denom, wv3, s

'   IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN

    If overfl <> 0 Then
      If output > 3 Then
        Line0
        Msg "    WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF"
        Msg "    LEE AND LEE UPDATE."
      End If
      Exit Sub
    End If

'   IF DENOM IS ZERO THE SOLVER IS PROBABLY NEAR SOLUTION -
'   AVOID OVERFLOW AND CONTINUE WITH SAME JACOBIAN.

    If denom = zero Then Exit Sub

'   THE SCALED VERSION OF S IS TAKEN TO THE UPDATE

    For k = 1 To n
      w(k) = w(k) * scalex(k) * scalex(k) / denom
    Next k

'   UPDATE THE QR DECOMPOSITION USING A SERIES OF GIVENS (JACOBI) ROTATIONS

    qrupda overfl, maxexp, n, epsmch, a, jac, t, w

'   RESET RDIAG AS DIAGONAL OF CURRENT R WHICH IS IN THE UPPER TRIANGLE OF A

    For i = 1 To n
      rdiag(i) = a((i - 1) * ISIZE + i)
    Next i

'   UPDATE P MATRIX

    denom = omega ^ (isejac + 2) + denom
    plee(1) = plee(1) - wv3(1) * wv3(1) / denom
    For j = 2 To n
      For i = 1 To j - 1
        plee((i - 1) * ISIZE + j) = plee((i - 1) * ISIZE + j) - wv3(i) * wv3(j) / denom
        plee((j - 1) * ISIZE + i) = plee((i - 1) * ISIZE + j)
      Next i
      plee((j - 1) * ISIZE + j) = plee((j - 1) * ISIZE + j) - wv3(j) * wv3(j) / denom
    Next j
  End If

'   UPDATE THE GRADIENT VECTOR, DELF.
'   DELF = (QR)^F = R^Q^F = R^JAC F

  If sclfch <> 0 Then
    For k = 1 To n
      w(k) = fvec(k) * scalef(k)
    Next k
  Else
    For k = 1 To n
      temp1((k - 1) * ISIZE + 1) = fvec(k)
      temp2((k - 1) * ISIZE + 1) = w(k)
    Next k
    matcop n, n, 1, 1, n, 1, temp1, temp2
    For k = 1 To n
      w(k) = temp2((k - 1) * ISIZE + 1)
    Next k
  End If
  avmul n, n, n, n, jac, w, t
  mmulv n, a, t, delf

End Sub

Sub llun(overch As Integer, overfl As Integer, isejac, maxexp, n, nunit, _
          output As Integer, epsmch, omega, fvec(), fvecc(), jac() As Double, _
          plee(), s(), scalex(), xc(), xplus())
'--------------------------------------------------------
'
'    FEB. 13, 1991
'
'    UPDATE THE JACOBIAN USING THE LEE AND LEE METHOD.
'--------------------------------------------------------
Dim wv1(ISIZE), temp(ISIZE), temp1(ISIZE)
Dim ss As String

  sqrtep = Sqr(epsmch)
  For i = 1 To n
    s(i) = (xplus(i) - xc(i)) * scalex(i)
  Next i
  For i = 1 To n
    For j = 1 To n
      temp(j) = plee((j - 1) * ISIZE + i)
    Next j
    wv1(i) = DotProduct(n, s, temp)
  Next i

  innerp overch, overfl, maxexp, n, n, n, nunit, output, denom, wv1, s

' IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN

  If overfl <> 0 Then
    If output > 3 Then
      Line0
      Msg "    WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF"
      Msg "    LEE AND LEE UPDATE."
    End If
    Exit Sub
  End If

'  IF DENOM IS ZERO THE SOLVER MUST BE VERY NEAR SOLUTION -
'  AVOID OVERFLOW AND CONTINUE WITH SAME JACOBIAN

  If denom = zero Then Exit Sub
  For i = 1 To n
    For j = 1 To n
      temp(i) = jac((i - 1) * ISIZE + j)
      temp1(i) = xplus(j) - xc(j)
    Next j
    Sum = DotProduct(n, temp, temp1)
    tempi = fvec(i) - fvecc(i) - Sum
    If Abs(tempi) >= sqrtep * (Abs(fvec(i)) + Abs(fvecc(i))) Then
      tempi = tempi / denom
      For j = 1 To n
        jac((i - 1) * ISIZE + j) = jac((i - 1) * ISIZE + j) + tempi * wv1(j) * scalex(j)
      Next j
    End If
  Next i

' UPDATE P MATRIX

  denom = omega ^ (isejac + 2) + denom
  plee(1) = plee(1) - wv1(1) * wv1(1) / denom
  For j = 2 To n
    For i = 1 To j - 1
      plee((i - 1) * ISIZE + j) = plee((i - 1) * ISIZE + j) - wv1(i) * wv1(j) / (denom * scalex(i) * scalex(j))
      plee((j - 1) * ISIZE + i) = plee((i - 1) * ISIZE + j)
    Next i
    plee((j - 1) * ISIZE + j) = plee((j - 1) * ISIZE + j) - wv1(j) * wv1(j) / denom
  Next j

End Sub

Sub matprt(nrowpr, ncolpr, a())
'-----------------------------------------------------------------------------
'    FEB. 6, 1991
'
'    THIS SUBROUTINE PRINTS RECTANGULAR BLOCKS STARTING WITH ELEMENT A(1,1)
'    OF SIZE NROWPR BY NCOLPR FOR MATRIX A (WHICH HAS DECLARED SIZE NROWA BY
'    NCOLA).  THE MATRIX IS PRINTED AS A BLOCK FOR SIZES UP TO 5X5 OR BY
'    COLUMNS IF IT IS LARGER.
'
'    NROWPR IS THE NUMBER OF ROWS TO BE PRINTED
'    NCOLPR IS THE NUMBER OF COLUMNS TO BE PRINTED
'
'    IF MATRIX PRINTING IS TO BE SUPPRESSED THEN LOGICAL
'    VARIABLE MATSUP MUST BE SET TO TRUE BEFORE THE CALL TO NNES.
'-----------------------------------------------------------------------------
Dim s As String
Dim s1 As String
Dim s2 As String

  If matsup <> 0 Then
    Line0
    Msg "       MATRIX PRINTING SUPPRESSED."
    Exit Sub
  End If

' FOR NCOLPR <= 5 WRITE MATRIX AS A WHOLE

  s = "": s1 = "": s2 = ""
  Line0
  If ncolpr <= 5 Then
    s = "           "
    For k = 1 To ncolpr
      s = Str$(k) + "        "
    Next k
    Msg s
    Line0
    For i = 1 To nrowpr
      s = Str$(i) + "     "
      For k = 1 To ncolpr
        s = s + Str$(a((i - 1) * ISIZE + k)) + " "
      Next k
      Msg s
    Next i
    Line0
  Else

'   LIMIT IS THE NUMBER OF GROUPS OF 5 COLUMNS
    limit = ncolpr / 5

'   WRITE COMPLETE BLOCKS FIRST (LEFTOVERS LATER)
    For j = 1 To limit
      s = "           "
      For k = 1 + (j - 1) * 5 To 5 + (j - 1) * 5
        s = s + Str$(k) + "        "
      Next k
      Msg s
      Line0
      For i = 1 To nrowpr
        s = Str$(i) + "     "
        For k = 1 + (j - 1) * 5 To 5 + (j - 1) * 5
          s = s + Str$(a((i - 1) * ISIZE + k)) + " "
        Next k
        Msg s
      Next i
      Line0
    Next j

'   WRITE REMAINING ELEMENTS
    s = "           "
    For k = 5 * limit + 1 To ncolpr
      s = s + Str$(k) + "        "
    Next k
    Msg s
    Line0
    For i = 1 To nrowpr
      s = Str$(i) + "     "
      For k = 5 * limit + 1 To ncolpr
        s = s + Str$(a((i - 1) * ISIZE + k)) + " "
      Next k
      Msg s
    Next i
    Line0
  End If

End Sub

Sub maxst(overfl As Integer, maxexp, n, nunit, output As Integer, _
          epsmch, maxstp As Double, mstpf As Double, scalex(), xc())
'--------------------------------------------------------------------------
'
'    FEB. 11, 1991
'
'    THIS SUBROUTINE ESTABLISHES A MAXIMUM STEP LENGTH BASED ON THE 2-NORMS
'    OF THE INITIAL ESTIMATES AND THE SCALING FACTORS MULTIPLIED BY A
'    FACTOR MSTPF.
'
'              maxstp = mstpf * Max(norm1, norm2)
'
'              WHERE   MSTPF  USER-CHOSEN FACTOR (DEFAULT: 1000)
'                      NORM1  2-NORM OF SCALED STARTING ESTIMATES
'                      NORM2  2-NORM OF COMPONENT SCALING FACTORS
'--------------------------------------------------------------------------
Dim norm1 As Double
Dim norm2 As Double
Dim wv1(ISIZE)
Dim s As String

  overfl = 0
  For k = 1 To n
    wv1(k) = scalex(k) * xc(k)
  Next k

  twonrm overfl, maxexp, n, epsmch, norm1, wv1

  If overfl <> 0 Then
    maxstp = ten ^ maxexp
    If output > 2 And wrnsup = 0 Then
      Line0
      s = "    WARNING: NORM OF SCALED INITIAL ESTIMATE SET TO " + Str$(norm1): Msg s
      Line0
      s = "       MAXIMUM STEP SIZE, MAXSTP, SET TO " + Str$(Int(maxstp * 1000) / 1000): Msg s
    End If
    Exit Sub
  End If

  twonrm overfl, maxexp, n, epsmch, norm2, scalex

  If overfl <> 0 Then
    maxstp = ten ^ maxexp
    If output > 2 And wrnsup = 0 Then
      Line0
      s = "    WARNING: NORM OF SCALING FACTORS SET TO " + Str$(norm2): Msg s
      Line0
      s = "       MAXIMUM STEP SIZE, MAXSTP, SET TO " + Str$(Int(maxstp * 1000) / 1000): Msg s
    End If
    Exit Sub
  End If

  maxstp = mstpf * XMAX(norm1, norm2)

End Sub

Sub nersl(newton, restrt As Integer, sclfch As Integer, sclxch As Integer, _
           acpcod As Integer, jupdm, n, nunit, output As Integer, fcnnew, _
           fvec(), xplus())
'-------------------------------------------------------
'    SEPT. 2, 1991
'
'    THE RESULTS OF EACH ITERATION ARE PRINTED.
'-------------------------------------------------------
Dim s As String
  Line0
  Line0
  If sclxch = 0 Then
    Msg "    SUMMARY OF ITERATION RESULTS"
  Else
    Msg "    SUMMARY OF ITERATION RESULTS (X'S GIVEN IN UNSCALED UNITS)"
  End If
  Line0
  Msg "        UPDATED ESTIMATES               UPDATED FUNCTION VALUES"
  Line0
  For i = 1 To n
    s = "      X(" + Str$(i) + ") = " + Str$(xplus(i)) + "  F(" + Str$(i) + ") = " + Str$(fvec(i)): Msg s
  Next i
  Line0
  If sclfch = 0 Then
    s = "      OBJECTIVE FUNCTION VALUE: " + Str$(fcnnew): Msg s
  Else
    s = "      SCALED OBJECTIVE FUNCTION VALUE: " + Str$(fcnnew): Msg s
  End If
  Line0
  If output > 3 And newton = 0 Then
    s = "      STEP ACCEPTANCE CODE, ACPCOD: " + Str$(acpcod): Msg s
  End If
  If restrt <> 0 And output >= 3 And jupdm <> 0 Then
    If output > 3 Then Line0
    Msg "      NOTE: JACOBIAN EVALUATED EXPLICITLY AT THIS STEP."
    Line0
  End If
End Sub

Sub nestop(absnew As Integer, linesr, newton, sclfch As Integer, sclxch As Integer, _
            acptcr As Integer, itnum, n, nac1, nac2, nac12, nfunc, njetot, nunit, _
            output As Integer, stopcr As Integer, trmcod As Integer, fcnnew, ftol, _
            nsttol As Double, stpmax, stptol, fvec(), scalef(), scalex(), xc(), xplus())
'----------------------------------------------------------------------------
'      FEB. 23, 1992
'
' *** THIS SUBROUTINE CHECKS TO SEE IF THE CONVERGENCE CRITERIA HAVE BEEN MET
'     AND PRINTS MAIN RESULTS TO OUTPUT FILE.  ***
'----------------------------------------------------------------------------
'Label: 100
Dim max1 As Double
Dim max2 As Double
Dim s As String

  If output > 3 Then
    Line0
    Msg "    CONVERGENCE TESTING"
    Line0
  End If

'   IF THE NEWTON STEP WAS WITHIN TOLERANCE THEN TRMCOD IS 1

  If trmcod = 1 Then
    If output > 3 Then
      Line0
      If sclxch = 0 Then
        s = "    MAXIMUM NEWTON STEP LENGTH STPMAX: " + Str$(stpmax): Msg s
      Else
        s = "    MAXIMUM SCALED NEWTON STEP LENGTH STPMAX: " + Str$(stpmax): Msg s
      End If
      s = "      FIRST CONVERGENCE CRITERION MET;  NSTTOL IS: " + Str$(nsttol): Msg s
      Line0

'   SKIP CHECKING OTHER STEP SIZE CRITERION AS TRMCOD IS ALREADY 1
      GoTo 100
    End If
  End If

'   IF THE NEWTON STEP WAS NOT WITHIN TOLERANCE THEN, IF
'   STOPCR IS NOT EQUAL TO 2, THE SECOND STEP SIZE STOPPING
'   CRITERION MUST BE CHECKED.

  If stopcr <> 2 And trmcod <> 1 Then
    max1 = zero
    For i = 1 To n
      ratio1 = (Abs(xplus(i) - xc(i))) / XMAX(Abs(xplus(i)), one / scalex(i))
      max1 = XMAX(max1, ratio1)
      If output > 4 Then
        s = "    RELATIVE STEP SIZE (" + Str$(i) + ") = " + Str$(ratio1): Msg s
      End If
    Next i
    If output > 4 Then Line0
    If output > 3 Then
      If sclxch = 0 Then
        s = "      MAXIMUM STEP SIZE: " + Str$(max1): Msg (s)
        s = "      (STPTOL: " + Str$(stptol) + ")": Msg s
      Else
        s = "      MAXIMUM RELATIVE STEP SIZE: " + Str$(max1): Msg s
        s = "      (STPTOL: " + Str$(stptol) + ")": Msg s
      End If
      Line0
    End If
    If max1 < stptol Then trmcod = 1
  End If

'   NOTE: CONTINUATION AT LABEL 100 MEANS THAT TRMCOD WAS 1 ON ENTRY
'         SO THE STEP SIZE CRITERION ABOVE DID NOT NEED TO BE CHECKED.

'   THE SECOND STOPPING CRITERION IS CHECKED IF NEEDED.

100 If stopcr = 2 Or stopcr = 12 Or (stopcr = 3 And trmcod = 1) Then
    
    max2 = zero

    For i = 1 To n
      max2 = XMAX(max2, scalef(i) * Abs(fvec(i)))
      If output > 4 Then
        If sclfch = 0 Then
          s = "   ABSOLUTE FUNCTION VECTOR (" + Str$(i) + ") = " + Str$(Abs(fvec(i))): Msg s
        Else
          s = "   ABSOLUTE SCALED FUNCTION VECTOR (" + Str$(i) + ") = " + Str$(scalef(i) * Abs(fvec(i))): Msg s
        End If
      End If
    Next i

    If output > 4 Then Line0
    If output > 3 Then
      If sclfch = 0 Then
        s = "      MAXIMUM ABSOLUTE FUNCTION: " + Str$(max2): Msg s
        s = "      (FTOL: " + Str$(ftol) + ")": Msg s
      Else
        s = "      MAX ABSOLUTE SCALED FUNCTION: " + Str$(max2): Msg s
        s = "      (FTOL: " + Str$(ftol) + ")": Msg s
      End If
      Line0
    End If

    If max2 < ftol Then
      If stopcr = 3 And trmcod = 1 Then

'     BOTH NEEDED STOPPING CRITERIA HAVE BEEN MET
        trmcod = 3

      ElseIf stopcr = 12 And trmcod = 1 Then

'     BOTH STOPPING CRITERIA HAVE BEEN MET ALTHOUGH
'     EITHER ONE WOULD BE SATISFACTORY.
        trmcod = 12

      ElseIf stopcr = 2 Or stopcr = 12 Then
        trmcod = 2
      End If
    ElseIf stopcr = 3 Then

'     ONLY THE FIRST STOPPING CRITERION WAS MET - TRMCOD
'     MUST BE RESET FROM 1 BACK TO 0.
      trmcod = 0
    End If
  End If

'    PRINT FINAL RESULTS IF CONVERGENCE REACHED

  If trmcod > 0 Then
    If output > 0 Then
      If output = 1 Then Line1
      Line0
      s = "    CONVERGENCE REACHED; TERMINATION CODE: ..............." + Str$(trmcod): Msg s

'   ITERATION RESULTS NOT PRINTED IN NERSL

      Line0
      Line0
      If sclfch <> 0 Or sclxch <> 0 Then
        Msg "                           UNSCALED RESULTS"
        Line0
      End If
      Msg "          FINAL ESTIMATES                 FINAL FUNCTION VALUES"
      Line0
      For i = 1 To n
        s = "    X(" + Str$(i) + ") = " + Str$(Int(xplus(i) * 100000) / 100000) + "    F(" + Str$(i) + ") = " + _
        Str$(fvec(i))
        Msg s
      Next i
      Line0

      If sclfch <> 0 Then
'       NEED UNSCALED OBJECTIVE FUNCTION
        Sum = DotProduct(n, fvec, fvec)
        s = "    FINAL OBJECTIVE FUNCTION VALUE: " + Str$(Sum / two): Msg s
      Else
        s = "    FINAL OBJECTIVE FUNCTION VALUE: " + Str$(fcnnew): Msg s
      End If

      If sclfch <> 0 Or sclxch <> 0 Then
        Line0
        Line0
        Msg "                           SCALED RESULTS"
        Line0
        Msg "          FINAL ESTIMATES                  FINAL FUNCTION VALUES"
        Line0
        For i = 1 To n
          s = "    X(" + Str$(i) + ") = " + Str$(scalex(i) * xplus(i)) + "  F(" + Str$(i) + ") = " + Str$(scalef(i) * fvec(i)):
          Msg s
        Next i
        Line0
        s = "    FINAL OBJECTIVE FUNCTION VALUE: " + Str$(fcnnew): Msg s
      End If
    End If
  Else
    If output > 3 Then
      Msg "    CONVERGENCE NOT REACHED."
      Line0
    End If
    Exit Sub
  End If

' TERMINATION HAS BEEN REACHED

  If output > 0 Then
    Line0
    Line0
    s = "    TOTAL NUMBER OF ITERATIONS: .........................." + Str$(itnum): Msg s
    If newton = 0 And absnew = 0 Then
      If linesr <> 0 Then
        s = "    TOTAL NUMBER OF LINE SEARCH FUNCTION EVALUATIONS: ...." + Str$(nfunc): Msg s
      Else
        s = "    TOTAL NUMBER OF TRUST REGION FUNCTION EVALUATIONS: ..." + Str$(nfunc): Msg s
      End If
    End If
    s = "    TOTAL NUMBER OF EXPLICIT JACOBIAN EVALUATIONS: ......." + Str$(njetot): Msg s
    s = "    TOTAL NUMBER OF FUNCTION EVALUATIONS: ................" + Str$(nfetot): Msg s
    Line0

    If newton = 0 And absnew = 0 And acptcr <> 1 And output > 2 Then
      s = "    NUMBER OF STEPS ACCEPTED BY FUNCTION VALUE ONLY: ....." + Str$(nac1): Msg s
      s = "    NUMBER OF STEPS ACCEPTED BY STEP SIZE VALUE ONLY: ...." + Str$(nac2): Msg s
      s = "    NUMBER OF STEPS ACCEPTED BY EITHER CRITERION: ........" + Str$(nac12): Msg s
      Line0
    End If
    Line1
  End If

End Sub