Attribute VB_Name = "Module5" '****************************************************************************** '* 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 1.1 By J-P Moreau, Paris. * '* (PART 3/3) * '* (www.jpmoreau.fr) * '* Release 1.1: corrected bug in sub qrdecomp * '****************************************************************************** DefDbl A-H, O-Z DefInt I-N Public output As Integer Public trmcod As Integer Public jactyp As Integer Public lam0 As Double Public maxexp As Integer Public maxit As Integer Public njetot As Integer Public help As String Public epsmch As Double Sub qform(n, a(), hhpi(), jac() As Double) '----------------------------------------------------- ' FEB. 14, 1991 ' ' FORM Q^ FROM THE HOUSEHOLDER MATRICES STORED IN ' MATRICES A AND HHPI AND STORE IT IN JAC. '----------------------------------------------------- Dim V1(ISIZE), V2(ISIZE) For i = 1 To n For j = 1 To n jac((i - 1) * ISIZE + j) = zero Next j Next i For j = 1 To n jac((j - 1) * ISIZE + j) = one Next j For k = 1 To n - 1 If hhpi(k) <> zero Then For j = 1 To n For i = k To n V1(i) = a((i - 1) * ISIZE + k) V2(i) = jac((i - 1) * ISIZE + j) Next i tau = DotProduct(n, V1, V2) tau = tau / hhpi(k) For i = k To n jac((i - 1) * ISIZE + j) = jac((i - 1) * ISIZE + j) - tau * a((i - 1) * ISIZE + k) Next i Next j End If Next k End Sub Sub qrdcom(qrsing As Integer, n, epsmch, a(), hhpi(), rdiag()) '--------------------------------------------------------------------- ' FEB. 23, 1992 ' ' THIS SUBROUTINE COMPUTES THE QR DECOMPOSITION OF THE MATRIX A. ' THE DECOMPOSITION IS COMPLETED EVEN IF A SINGULARITY IS DETECTED ' (WHEREUPON QRSING IS SET TO TRUE). '--------------------------------------------------------------------- Dim V1(ISIZE), V2(ISIZE) qrsing = 0 For k = 1 To n - 1 eta = zero For i = k To n eta = XMAX(eta, Abs(a((i - 1) * ISIZE + k))) Next i If eta < epsmch Then qrsing = 1 hhpi(k) = zero rdiag(k) = zero Else For j = k To n a((j - 1) * ISIZE + k) = a((j - 1) * ISIZE + k) / eta Next j sigma = zero For j = k To n sigma = sigma + a((j - 1) * ISIZE + k) ^ 2 Next j sigma = Sign(Sqr(sigma), a((k - 1) * ISIZE + k)) a((k - 1) * ISIZE + k) = a((k - 1) * ISIZE + k) + sigma hhpi(k) = sigma * a((k - 1) * ISIZE + k) rdiag(k) = -eta * sigma For j = k + 1 To n For i = k To n V1(i - k + 1) = a((i - 1) * ISIZE + k) 'modified 04/10/2007 V2(i - k + 1) = a((i - 1) * ISIZE + j) ' idem Next i tau = DotProduct(n - k + 1, V1, V2) tau = tau / hhpi(k) For i = k To n a((i - 1) * ISIZE + j) = a((i - 1) * ISIZE + j) - tau * a((i - 1) * ISIZE + k) Next i Next j End If Next k rdiag(n) = a((n - 1) * ISIZE + n) If Abs(rdiag(n)) < epsmch Then qrsing = 1 End Sub Sub utumul(nradec, ncadec, nraact, ncaact, nrbdec, ncbdec, amat(), bmat()) '----------------------------------------------------------------------------- ' FEB. 8, 1991 ' ' MATRIX MULTIPLICATION: A^A=B WHERE A IS UPPER TRIANGULAR ' ' VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. ' ' NRADEC IS NUMBER OF ROWS IN A DECLARED ' NCADEC IS NUMBER OF COLUMNS IN A DECLARED ' NRAACT IS THE LIMIT FOR THE 1ST INDEX IN A ' NCAACT IS THE LIMIT FOR THE 2ND INDEX IN A ' NRBDEC IS NUMBER OF ROWS IN B DECLARED ' NCBDEC IS NUMBER OF COLUMNS IN B DECLARED ' ' MODIFIED VERSION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ' QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA '------------------------------------------------------------------------------ ' FIND ENTRY IN MATRIX B For i = 1 To ncaact ' FIND NUMBER OF GROUPS OF SIZE 32, 16... nend = IMIN(i, nraact) ncc32 = nend / 32 ncc32r = nend - 32 * ncc32 ncc16 = ncc32r / 16 ncc16r = ncc32r - 16 * ncc16 ncc8 = ncc16r / 8 ncc8r = ncc16r - 8 * ncc8 ncc4 = ncc8r / 4 ncc4r = ncc8r - 4 * ncc4 For j = i To ncaact Sum = zero k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 Sum = Sum + amat((k - 32) * ISIZE + i) * amat((k - 32) * ISIZE + j) + amat((k - 31) * ISIZE + i) * amat((k - 31) * ISIZE + j) _ + amat((k - 30) * ISIZE + i) * amat((k - 30) * ISIZE + j) + amat((k - 29) * ISIZE + i) * amat((k - 29) * ISIZE + j) + _ amat((k - 28) * ISIZE + i) * amat((k - 28) * ISIZE + j) + amat((k - 27) * ISIZE + i) * amat((k - 27) * ISIZE + j) + _ amat((k - 26) * ISIZE + i) * amat((k - 26) * ISIZE + j) + amat((k - 25) * ISIZE + i) * amat((k - 25) * ISIZE + j) Sum = Sum + amat((k - 24) * ISIZE + i) * amat((k - 24) * ISIZE + j) + amat((k - 23) * ISIZE + i) * amat((k - 23) * ISIZE + j) _ + amat((k - 22) * ISIZE + i) * amat((k - 22) * ISIZE + j) + amat((k - 21) * ISIZE + i) * amat((k - 21) * ISIZE + j) + _ amat((k - 20) * ISIZE + i) * amat((k - 20) * ISIZE + j) + amat((k - 19) * ISIZE + i) * amat((k - 19) * ISIZE + j) + _ amat((k - 18) * ISIZE + i) * amat((k - 18) * ISIZE + j) + amat((k - 17) * ISIZE + i) * amat((k - 17) * ISIZE + j) Sum = Sum + amat((k - 16) * ISIZE + i) * amat((k - 16) * ISIZE + j) + amat((k - 15) * ISIZE + i) * amat((k - 15) * ISIZE + j) _ + amat((k - 14) * ISIZE + i) * amat((k - 14) * ISIZE + j) + amat((k - 13) * ISIZE + i) * amat((k - 13) * ISIZE + j) + _ amat((k - 12) * ISIZE + i) * amat((k - 12) * ISIZE + j) + amat((k - 11) * ISIZE + i) * amat((k - 11) * ISIZE + j) + _ amat((k - 10) * ISIZE + i) * amat((k - 10) * ISIZE + j) + amat((k - 9) * ISIZE + i) * amat((k - 9) * ISIZE + j) Sum = Sum + amat((k - 8) * ISIZE + i) * amat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * amat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * amat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * amat((k - 5) * ISIZE + j) + amat((k - 4) * ISIZE + i) _ * amat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * amat((k - 3) * ISIZE + j) + amat((k - 2) * ISIZE + i) * _ amat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * amat((k - 1) * ISIZE + j) Next kk End If If ncc16 > 0 Then For kk = 1 To ncc16 k = k + 16 Sum = Sum + amat((k - 16) * ISIZE + i) * amat((k - 16) * ISIZE + j) + amat((k - 15) * ISIZE + i) * amat((k - 15) * ISIZE + j) _ + amat((k - 14) * ISIZE + i) * amat((k - 14) * ISIZE + j) + amat((k - 13) * ISIZE + i) * amat((k - 13) * ISIZE + j) + _ amat((k - 12) * ISIZE + i) * amat((k - 12) * ISIZE + j) + amat((k - 11) * ISIZE + i) * amat((k - 11) * ISIZE + j) + _ amat((k - 10) * ISIZE + i) * amat((k - 10) * ISIZE + j) + amat((k - 9) * ISIZE + i) * amat((k - 9) * ISIZE + j) Sum = Sum + amat((k - 8) * ISIZE + i) * amat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * amat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * amat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * amat((k - 5) * ISIZE + j) + amat((k - 4) * ISIZE + i) _ * amat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * amat((k - 3) * ISIZE + j) + amat((k - 2) * ISIZE + i) * _ amat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * amat((k - 1) * ISIZE + j) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 Sum = Sum + amat((k - 8) * ISIZE + i) * amat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * amat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * amat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * amat((k - 5) * ISIZE + j) + amat((k - 4) * ISIZE + i) _ * amat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * amat((k - 3) * ISIZE + j) + amat((k - 2) * ISIZE + i) * _ amat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * amat((k - 1) * ISIZE + j) Next kk End If If ncc4 > 0 Then For kk = 1 To ncc4 k = k + 4 Sum = Sum + amat((k - 4) * ISIZE + i) * amat((k - 4) * ISIZE + j) + _ amat((k - 3) * ISIZE + i) * amat((k - 3) * ISIZE + j) + amat((k - 2) * ISIZE + i) * _ amat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * amat((k - 1) * ISIZE + j) Next kk End If If ncc4r > 0 Then For kk = 1 To ncc4r k = k + 1 Sum = Sum + amat((k - 1) * ISIZE + i) * amat((k - 1) * ISIZE + j) Next kk End If bmat((i - 1) * ISIZE + j) = Sum If i <> j Then bmat((j - 1) * ISIZE + i) = bmat((i - 1) * ISIZE + j) Next j Next i End Sub Sub rtrmul(n, a(), h(), rdiag()) '-------------------------------------------------------- ' SEPT. 4, 1991 ' ' FIND R^R FOR QR-DECOMPOSED JACOBIAN. ' ' R IS STORED IN STRICT UPPER TRIANGLE OF A AND RDIAG. '-------------------------------------------------------- Dim wv1(ISIZE) ' TEMPORARILY REPLACE DIAGONAL OF R IN A (A IS RESTORED LATER). For i = 1 To n wv1(i) = a((i - 1) * ISIZE + i) a((i - 1) * ISIZE + i) = rdiag(i) Next i utumul n, n, n, n, n, n, a, h For i = 1 To n a((i - 1) * ISIZE + i) = wv1(i) Next i End Sub Sub onenrm(abort As Integer, pertrb As Integer, n, nunit, output As Integer, _ epsmch, h1norm, h(), scalex()) '---------------------------------------------------------------------------- ' FEB. 23, 1992 ' ' FIND 1-NORM OF H MATRIX IF PERTURBATION IS DESIRED AND PERTURB DIAGONAL. '---------------------------------------------------------------------------- Dim s As String sqrtep = Sqr(epsmch) If output > 4 Then Line0 Msg " DIAGONAL OF MATRIX H (= JAC^JAC) BEFORE BEING PERTURBED:" Line0 For i = 1 To n s = " H(" + Str$(i) + "," + Str$(i) + ") = " + Str$(h((i - 1) * ISIZE + i)): Msg s Next i End If h1norm = zero For j = 1 To n h1norm = h1norm + Abs(h(j)) / scalex(j) Next j h1norm = h1norm / scalex(1) For i = 2 To n temp = zero For j = 1 To i temp = temp + Abs(h((j - 1) * ISIZE + i)) / scalex(j) Next j For j = i + 1 To n temp = temp + Abs(h((i - 1) * ISIZE + j)) / scalex(j) Next j h1norm = XMAX(h1norm, temp / scalex(i)) Next i If output > 4 Then Line0 s = " 1-NORM OF MATRIX H: " + Str$(h1norm): Msg s End If If h1norm < epsmch Then If output > 0 Then Line1 Line0 Msg " PROGRAM FAILS AS 1-NORM OF JACOBIAN IS TOO SMALL." Line0 Line1 End If abort = 1 Exit Sub Else ' PERTURB DIAGONAL OF MATRIX H - USE THIS TO FIND "SN" pertrb = 1 For i = 1 To n h((i - 1) * ISIZE + i) = h((i - 1) * ISIZE + i) + Sqr(1# * n) * sqrtep * h1norm * scalex(i) * scalex(i) Next i If output > 4 Then Line0 Line0 Msg " PERTURBED H MATRIX:" matprt n, n, h End If abort = 0 End If End Sub Sub nstpun(abort As Integer, linesr, overch As Integer, overfl As Integer, _ qrsing As Integer, sclfch As Integer, sclxch As Integer, itnum, _ maxexp, n, nunit, output As Integer, epsmch, a(), delf(), fvecc(), _ h(), hhpi(), jac() As Double, rdiag(), scalef(), scalex(), sn()) '-------------------------------------------------------------------------------- ' FEB. 23, 1992 ' ' THIS SUBROUTINE FINDS THE NEWTON STEP. ' ' IF THE JACOBIAN IS DETECTED AS SINGULAR OR IF THE ESTIMATED CONDITION ' NUMBER IS TOO HIGH (GREATER THAN EPSMCH^(-2/3) ) THEN H=J^J IS FORMED ' AND THE DIAGONAL IS PERTURBED BY ADDING SQRT(N*EPSMCH)*H1NORM*SCALEX(I)^2 ' TO THE CORRESPONDING ELEMENT. A CHOLESKY DECOMPOSITION IS PERFORMED ON ' THIS MODIFIED MATRIX PRODUCING A PSEUDO-NEWTON STEP. ' NOTE: THIS PROCEDURE MAY BE BE BYPASSED FOR ILL-CONDITIONED JACOBIANS ' BY SETTING THE LOGICAL VARIABLE BYPASS TO TRUE IN THE DRIVER. ' ' ABORT IF THE 1-NORM OF MATRIX H BECOMES TOO SMALL ' ALTHOUGH BUT NOT AT A SOLUTION ' BYPASS ALLOWS BYPASSING OF THE SPECIAL TREATMENT FOR ' BADLY CONDITIONED JACOBIANS ' QRSING INDICATES SINGULAR JACOBIAN DETECTED. '-------------------------------------------------------------------------------- Dim maxadd As Double Dim maxffl As Double Dim wv1(ISIZE) Dim pertrb As Integer Dim s As String newstm = 0 abort = 0 overfl = 0 pertrb = 0 sqrtep = Sqr(epsmch) If output > 3 Then Line0 Line0 Msg " SOLUTION OF LINEAR SYSTEM FOR NEWTON STEP, SN" End If ' STORE (POSSIBLY SCALED) JACOBIAN IN MATRIX A If sclfch = 0 Then matcop n, n, n, n, n, n, jac, a Else For i = 1 To n If scalef(i) <> one Then scalfi = scalef(i) For j = 1 To n a((i - 1) * ISIZE + j) = jac((i - 1) * ISIZE + j) * scalfi Next j Else For j = 1 To n a((i - 1) * ISIZE + j) = jac((i - 1) * ISIZE + j) Next j End If Next i End If ' SCALED JACOBIAN IS PRINTED ONLY IF AT LEAST ONE SCALING ' FACTOR IS NOT 1.0 If output > 4 And sclfch <> 0 Then Line0 Line0 Msg " SCALED JACOBIAN MATRIX:" matprt n, n, a End If ' QR DECOMPOSITION OF (POSSIBLY SCALED) JACOBIAN qrdcom qrsing, n, epsmch, a, hhpi, rdiag ' SAVE MATRIX A FOR BACK SUBSTITUTION TO CHECK DEUFLHARDS'S SECOND ' STEP ACCEPTANCE CRITERION IN LINE SEARCH OR TRUST REGION METHOD. matcop n, n, n, n, n, n, a, h If output > 4 And n > 1 And matsup = 0 Then Line0 If sclfch = 0 Then Msg " QR DECOMPOSITION OF JACOBIAN MATRIX:" Else Msg " QR DECOMPOSITION OF SCALED JACOBIAN MATRIX:" End If matprt n, n, a Line0 Msg " DIAGONAL OF R PI FACTORS FROM QR DECOMPOSITION" Line0 For i = 1 To n - 1 s = " RDIAG(" + Str$(i) + ") = " + Str$(rdiag(i)) + " HHPI(" + Str$(i) + ") = " + Str$(hhpi(i)): Msg s Next i s = " RDIAG(" + Str$(n) + ") = " + Str$(rdiag(n)): Msg s Line0 If itnum = 1 Then Msg " NOTE: R IS IN STRICT UPPER TRIANGLE OF MATRIX A PLUS RDIAG." Line0 Msg " THE COLUMNS OF THE LOWER TRIANGLE OF MATRIX A PLUS" Msg " THE ELEMENTS OF VECTOR HHPI FORM THE HOUSEHOLDER" Msg " MATRICES WHICH, WHEN MULTIPLIED TOGETHER, FORM Q." End If End If ' ESTIMATE CONDITION NUMBER IF (SCALED) JACOBIAN IS NOT SINGULAR If bypass = 0 And qrsing = 0 And n > 1 Then If sclxch <> 0 Then ' SET UP FOR CONDITION NUMBER ESTIMATOR - SCALE R WRT X'S For j = 1 To n If scalex(j) <> one Then scalxj = scalex(j) rdiag(j) = rdiag(j) / scalxj For i = 1 To j - 1 a((i - 1) * ISIZE + j) = a((i - 1) * ISIZE + j) / scalxj Next i End If Next j End If condno overch, overfl, maxexp, n, nunit, output, connum, a, rdiag ' IF OVERFLOW DETECTED IN CONDITION NUMBER ESTIMATOR, ASSIGN QRSING AS TRUE If overfl <> 0 Then qrsing = 1 ' NOTE: OVERFL SWITCHED TO FALSE BEFORE FORMATION OF H LATER Else ' ASSIGN DUMMY TO CONNUM FOR SINGULAR JACOBIAN UNLESS N=1 If n = 1 Then connum = one Else connum = zero End If ' MATRIX H=JAC^JAC IS FORMED IN THREE CASES: ' 1) THE (SCALED) JACOBIAN IS SINGULAR ' 2) THE CONDITION NUMBER IS TOO HIGH AND THE OPTION TO BYPASS THE ' PERTURBATION OF THE JACOBIAN IS NOT BEING USED. ' 3) REQUESTED BY USER THROUGH NEWSTM. End If If qrsing <> 0 Or (bypass = 0 And connum > one / sqrtep ^ 1.333) Or newstm = 77 Then ' FORM H=(DF*JAC]^(DF*JAC) WHERE DF=DIAG(SCALEF). USE ' PREVIOUSLY COMPUTED QR DECOMPOSITION OF (SCALED) JACOBIAN ' WHERE R IS STORED IN THE UPPER TRIANGLE OF A AND RDIAG. If overch <> 0 Then overfl = 0 ataov overfl, maxexp, n, nunit, output, jac, h, scalef Else rtrmul n, a, h, rdiag End If If output > 3 Then Line0 If qrsing <> 0 And overfl = 0 Then Msg " SINGULAR JACOBIAN DETECTED: JACOBIAN PERTURBED." Else ' NOTE: IF OVERFL IS TRUE THEN QRSING MUST BE TRUE If overfl <> 0 Then Msg " POTENTIAL OVERFLOW DETECTED IN CONDITION NUMBER ESTIMATOR." Msg " MATRIX AS SINGULAR AND JACOBIAN PERTURBED." Else s = " CONDITION NUMBER TOO HIGH: " + Str$(connum) + ", JACOBIAN PERTURBED.": Msg s End If End If End If overfl = 0 If newstm <> 77 Then ' FIND 1-NORM OF H MATRIX AND PERTURB DIAGONAL onenrm abort, pertrb, n, nunit, output, epsmch, h1norm, h, scalex If abort <> 0 Then Exit Sub End If ' CHOLESKY DECOMPOSITION OF H MATRIX - MAXFFL=0 IMPLIES ' THAT H IS KNOWN TO BE POSITIVE DEFINITE. maxffl = zero cholde n, maxadd, maxffl, sqrtep, h, a If output > 4 Then Line0 Line0 Msg " CHOLESKY DECOMPOSITION OF H MATRIX:" matprt n, n, a End If ' FIND NEWTON STEP FROM CHOLESKY DECOMPOSITION. IF THE DIAGONAL HAS ' BEEN PERTURBED THEN THIS IS NOT THE ACTUAL NEWTON STEP BUT ' ONLY AN APPORXIMATION THEREOF. For i = 1 To n wv1(i) = -delf(i) Next i chsolv overch, overfl, maxexp, n, nunit, output, a, wv1, sn overfl = 0 If output > 3 Then If sclxch = 0 Then Line0 If pertrb = 0 Then Msg " NEWTON STEP FROM CHOLESKY DECOMPOSITION" Else Msg " APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN" End If Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)): Msg s Next i Else Line0 If pertrb = 0 Then Msg " NEWTON STEP FROM CHOLESKY DECOMPOSITION IN SCALED UNITS:" Else Msg " APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN IN SCALED UNITS:" End If Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)) + " SN(" + Str$(i) + ") = " + Str$(scalex(i) * sn(i)): Msg s Next i End If End If ' SET QRSING TO TRUE SO THAT THE CORRECT MATRIX FACTORIZATION IS USED ' IN THE BACK-CALCULATION OF SBAR FOR DEUFLHARD RELAXATION FACTOR ' INITIALIZATION IN LINE SEARCH (ONLY MATTERS WHEN JACOBIAN IS ILL- ' CONDITIONED BUT NOT SINGULAR). qrsing = 1 Else If output > 3 And n > 1 Then If bypass = 0 And connum <= one / (sqrtep ^ 1.333) Then Line0 s = " CONDITION NUMBER ACCEPTABLE: " + Str$(Int(connum * 10000) / 10000) + ", JACOBIAN NOT PERTURBED" Msg s End If If bypass <> 0 And connum > one / (sqrtep ^ 1.333) Then Line0 s = " CONDITION NUMBER HIGH: " + Str$(Int(connum * 10000) / 10000) + ", JACOBIAN NOT PERTURBED AS" Msg s Msg " BYPASS IS TRUE." End If End If ' NOTE: HERE SN STORES THE R.H.S. - IT IS OVERWRITTEN For i = 1 To n sn(i) = -fvecc(i) * scalef(i) Next i If bypass = 0 And sclxch <> 0 Then ' R WAS SCALED BEFORE THE CONDITION NUMBER ESTIMATOR - ' THIS CONVERTS IT BACK TO THE UNSCALED FORM. For j = 1 To n If scalex(j) <> one Then scalxj = scalex(j) rdiag(j) = rdiag(j) * scalxj For i = 1 To j - 1 a((i - 1) * ISIZE + j) = a((i - 1) * ISIZE + j) * scalxj Next i End If Next j End If ' ACCEPTABLE CONDITION NUMBER - USE BACK SUBSTITUTION TO ' FIND NEWTON STEP FROM QR DECOMPOSITION. qrsolv overch, overfl, maxexp, n, nunit, output, a, hhpi, rdiag, sn overfl = 0 If output > 3 Then If sclxch = 0 Then Line0 Msg " NEWTON STEP FROM QR DECOMPOSITION:" Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)): Msg s Next i Else Line0 Msg " NEWTON STEP FROM QR DECOMPOSITION IN SCALED UNITS:" Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)) + " SN(" + Str$(i) + ") = " + Str$(scalex(i) * sn(i)) Msg s Next i End If End If ' TRANSFORM MATRICES FOR SUBSEQUENT CALCULATIONS IN TRUST ' REGION METHODS (A IS STORED ABOVE IN H). If linesr = 0 Then For i = 1 To n a((i - 1) * ISIZE + i) = rdiag(i) For j = 1 To i - 1 a((i - 1) * ISIZE + j) = a((j - 1) * ISIZE + i) Next j Next i End If End If End Sub Sub nstpfa(abort As Integer, linesr, overch As Integer, overfl As Integer, qrsing As Integer, _ restrt As Integer, sclfch As Integer, sclxch As Integer, itnum, maxexp, n, newstm, _ nunit, output As Integer, epsmch, a(), delf(), fvecc(), h(), hhpi(), jac() As Double, _ rdiag(), scalef(), scalex(), sn()) '-------------------------------------------------------------------------------- ' ' FEB. 23, 1992 ' ' THIS SUBROUTINE FINDS THE NEWTON STEP. ' ' IF THE JACOBIAN IS DETECTED AS SINGULAR OR IF THE ESTIMATED CONDITION ' NUMBER IS TOO HIGH ^[GREATER THAN EPSMCH**^[-2/3]] THEN H:=J^J IS FORMED ' AND THE DIAGONAL IS PERTURBED BY ADDING SQRT^[N*EPSMCH]*H1NORM*SCALEX^[I]**2 ' TO THE CORRESPONDING ELEMENT. A CHOLESKY DECOMPOSITION IS PERFORMED ON ' THIS MODIFIED MATRIX PRODUCING A PSEUDO-NEWTON STEP. ' ' IF THE CONDITION NUMBER IS SMALL THEN THE NEWTON STEP, SN, ' IS FOUND DIRECTLY BY BACK SUBSTITUTION. ' ' ABORT IF THE 1-NORM OF MATRIX H BECOMES TOO SMALL ' ALTHOUGH BUT NOT AT A SOLUTION ' BYPASS ALLOWS BYPASSING OF THE SPECIAL TREATMENT FOR ' BADLY CONDITIONED JACOBIANS ' QRSING INDICATES SINGULAR JACOBIAN DETECTED '--------------------------------------------------------------------------------- Dim pertrb As Integer Dim maxadd As Double Dim maxffl As Double Dim wv1(ISIZE) Dim s As String abort = 0 overfl = 0 sqrtep = Sqr(epsmch) If restrt <> 0 Then If output > 3 Then Line0 Line0 Msg " SOLUTION OF LINEAR SYSTEM FOR NEWTON STEP, SN" End If ' STORE (POSSIBLY SCALED) JACOBIAN IN MATRIX A If sclfch = 0 Then matcop n, n, n, n, n, n, jac, a Else For i = 1 To n If scalef(i) <> one Then scalfi = scalef(i) For j = 1 To n a((i - 1) * ISIZE + j) = jac((i - 1) * ISIZE + j) * scalfi Next j Else For j = 1 To n a((i - 1) * ISIZE + j) = jac((i - 1) * ISIZE + j) Next j End If Next i End If ' SCALED JACOBIAN IS PRINTED ONLY IF AT LEAST ONE SCALING FACTOR IS NOT 1.0 If output > 4 And sclfch <> 0 Then Line0 Line0 Msg " SCALED JACOBIAN MATRIX:" matprt n, n, a End If ' QR DECOMPOSITION OF (POSSIBLY SCALED) JACOBIAN qrdcom qrsing, n, epsmch, a, hhpi, rdiag If output > 4 And n > 1 And matsup = 0 Then Line0 If sclfch = 0 Then Msg " QR DECOMPOSITION OF JACOBIAN MATRIX:" Else Msg " QR DECOMPOSITION OF SCALED JACOBIAN MATRIX:" End If matprt n, n, a Line0 Msg " DIAGONAL OF R PI FACTORS FROM QR DECOMPOSITION" Line0 For i = 1 To n - 1 s = " RDIAG(" + Str$(i) + ") = " + Str$(rdiag(i)) + " HHPI(" + Str$(i) + ") = " + Str$(hhpi(i)): Msg s Next i s = " RDIAG(" + Str$(n) + ") = " + Str$(rdiag(n)): Msg s Line0 If itnum = 1 Then Msg " NOTE: R IS IN STRICT UPPER TRIANGLE OF MATRIX A PLUS RDIAG" Line0 Msg " THE COLUMNS OF THE LOWER TRIANGLE OF MATRIX A PLUS" Msg " THE ELEMENTS OF VECTOR HHPI FORM THE HOUSEHOLDER" Msg " MATRICES WHICH, WHEN MULTIPLIED TOGETHER, FORM Q." End If End If ' FORM THE ACTUAL Q^ MATRIX FROM THE HOUSEHOLDER TRANSFORMATIONS STORED ' IN THE LOWER TRIANGLE OF A AND THE FACTORS IN HHPI: STORE IT IN JAC qform n, a, hhpi, jac ' COMPLETE THE UPPER TRIANGULAR R MATRIX BY REPLACING THE ' DIAGONAL OF A. THE QR DECOMPOSITION IS NOW AVAILABLE For i = 1 To n a((i - 1) * ISIZE + i) = rdiag(i) Next i Else ' USING UPDATED FACTORED FORM OF JACOBIAN - CHECK FOR SINGULARITY qrsing = 0 For i = 1 To n If a((i - 1) * ISIZE + i) = zero Then qrsing = 1 Next i End If ' ESTIMATE CONDITION NUMBER IF JACOBIAN IS NOT SINGULAR If bypass = 0 And qrsing = 0 And n > 1 Then If sclxch <> 0 Then ' SET UP FOR CONDITION NUMBER ESTIMATION - SCALE R WRT X'S For j = 1 To n If scalex(j) <> one Then scalxj = scalex(j) rdiag(j) = rdiag(j) / scalxj For i = 1 To j - 1 a((i - 1) * ISIZE + j) = a((i - 1) * ISIZE + j) / scalxj Next i End If Next j End If condno overch, overfl, maxexp, n, nunit, output, connum, a, rdiag ' UNSCALE R IF IT WAS SCALED BEFORE THE CALL TO CONDNO If sclxch <> 0 Then For j = 1 To n If scalex(j) <> one Then scalxj = scalex(j) rdiag(j) = rdiag(j) * scalxj For i = 1 To j - 1 a((i - 1) * ISIZE + j) = a((i - 1) * ISIZE + j) * scalxj Next i End If Next j End If ' IF OVERFLOW DETECTED IN CONDITION NUMBER ESTIMATOR ASSIGN ' QRSING AS TRUE SO THAT THE JACOBIAN WILL BE PERTURBED. If overfl <> 0 Then qrsing = 1 ' NOTE: OVERFL SWITCHED TO FALSE BEFORE FORMATION OF H Else If n = 1 Then connum = one Else connum = zero End If End If ' MATRIX H=JAC^JAC IS FORMED IN TWO CASES: ' 1) THE ^[SCALED] JACOBIAN IS SINGULAR ' 2) THE CONDITION NUMBER IS TOO HIGH AND THE OPTION TO BYPASS ' THE PERTURBATION OF THE JACOBIAN IS NOT BEING USED ' 3) REQUESTED BY THE USER THROUGH NEWSTM. If qrsing <> 0 Or (bypass = 0 And connum > one / (sqrtep ^ 1.333)) Or newstm = 77 Then ' FORM H=(DF*JAC)^(DF*JAC) WHERE DF=DIAG(SCALEF). USE ' PREVIOUSLY COMPUTED QR DECOMPOSITION OF (SCALED) JACOBIAN ' WHERE R IS STORED IN THE UPPER TRIANGLE OF A AND RDIAG. overfl = 0 If overch <> 0 Then ataov overfl, maxexp, n, nunit, output, jac, h, scalef Else rtrmul n, a, h, rdiag End If If output > 3 Then Line0 If qrsing <> 0 And overfl = 0 Then Msg " SINGULAR JACOBIAN DETECTED: JACOBIAN PERTURBED." Else If overfl <> 0 Then Msg " POTENTIAL OVERFLOW DETECTED IN CONDITION NUMBER ESTIMATOR," Msg " MATRIX AS SINGULAR AND JACOBIAN PERTURBED." Else s = " CONDITION NUMBER TOO HIGH: " + Str$(connum) + " JACOBIAN PERTURBED.": Msg s End If End If End If overfl = 0 If newstm <> 77 Then ' FIND 1-NORM OF H MATRIX AND PERTURB DIAGONAL onenrm abort, pertrb, n, nunit, output, epsmch, h1norm, h, scalex End If ' CHOLESKY DECOMPOSITION OF H MATRIX - MAXFFL=0 INDICATES ' THAT H IS KNOWN TO BE POSITIVE DEFINITE. maxffl = zero cholde n, maxadd, maxffl, sqrtep, h, a If output > 4 Then Line0 Line0 Msg " CHOLESKY DECOMPOSITION OF H MATRIX:" matprt n, n, a End If ' FIND NEWTON STEP FROM CHOLESKY DECOMPOSITION. IF THE ' DIAGONAL HAS BEEN PERTURBED THEN THIS IS NOT THE ACTUAL ' NEWTON STEP BUT ONLY AN APPROXIMATION THEREOF. For i = 1 To n wv1(i) = -delf(i) Next i chsolv overch, overfl, maxexp, n, nunit, output, a, wv1, sn overfl = 0 If output > 3 Then If sclxch = 0 Then If pertrb = 0 Then Msg " NEWTON STEP FROM CHOLESKY DECOMPOSITION" Else Msg " APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN" End If Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)): Msg s Next i Else Line0 If pertrb = 0 Then Msg " NEWTON STEP FROM CHOLESKY DECOMPOSITION IN SCALED UNITS" Else Msg " APPROXIMATE NEWTON STEP FROM PERTURBED JACOBIAN IN SCALED UNITS" End If Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)) + " SN(" + Str$(i) + ") = " + Str$(scalex(i) * sn(i)): Msg s Next i End If End If ' SET QRSING TO TRUE SO THAT THE CORRECT MATRIX FACTORIZATION ' IS USED IN THE BACK-CALCULATION OF SBAR FOR DEUFLHARD ' RELAXATION FACTOR INITIALIZATION. qrsing = 1 Else If output > 3 And n > 1 Then If bypass = 0 And connum <= one / (sqrtep ^ 1.333) Then Line0 s = " CONDITION NUMBER ACCEPTABLE: " + Str$(connum) + ", JACOBIAN NOT PERTURBED": Msg s End If If bypass <> 0 And connum > one / (sqrtep ^ 1.333) Then Line0 s = " CONDITION NUMBER HIGH: " + Str$(connum) + ", JACOBIAN NOT PERTURBED AS": Msg s Msg " BYPASS IS TRUE." End If End If For i = 1 To n Sum = zero For j = 1 To n Sum = Sum - jac((i - 1) * ISIZE + j) * scalef(j) * fvecc(i) Next j sn(i) = Sum Next i rsolv overch, overfl, maxexp, n, nunit, output, a, rdiag, sn overfl = 0 If output > 3 Then If sclxch = 0 Then Line0 Msg " NEWTON STEP FROM QR DECOMPOSITION" Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)): Msg s Next i Else Line0 Msg " NEWTON STEP FROM QR DECOMPOSITION IN SCALED UNITS" Line0 For i = 1 To n s = " SN(" + Str$(i) + ") = " + Str$(sn(i)) + " SN(" + Str$(i) + ") = " + Str$(scalex(i) * sn(i)): Msg s Next i End If End If ' TRANSFORM MATRICES FOR SUBSEQUENT CALCULATIONS IN TRUST REGION METHOD If linesr = 0 Then For i = 2 To n For j = 1 To i - 1 a((i - 1) * ISIZE + j) = a((j - 1) * ISIZE + i) Next j Next i End If End If End Sub Sub title(cauchy As Integer, deuflh As Integer, geoms As Integer, linesr, newton, _ overch As Integer, acptcr As Integer, contyp As Integer, itsclf, itsclx, _ jactyp, jupdm, maxit, maxns, maxqns, mgll, minqns, n, narmij, ninitn, _ njacch, nunit, output As Integer, qnupdm As Integer, stopcr As Integer, _ trupdm As Integer, alpha, confac, delfac, delta, epsmch, etafac, fcnold, ftol, _ lam0 As Double, maxstp As Double, mstpf As Double, nsttol As Double, omega, _ ratiof, sigma, stptol, boundl(), boundu(), fvecc(), scalef(), scalex(), xc()) '-------------------------------------------------------- ' APR. 13, 1991 ' ' THIS PROCEDURE WRITES TITLE AND RECORDS PARAMETERS. '-------------------------------------------------------- 'Label: 100 Dim s As String If output < 2 Then Exit Sub Print #1, Print #1, Print #1, Line1 Line1 Line0 Msg " N N E S" Line0 Msg " NONMONOTONIC NONLINEAR EQUATION SOLVER VERSION 1.05" Line0 Msg " COPYRIGHT 1991, BY R.S. BAIN" Line0 Line1 Line1 Print #1, Print #1, If output < 3 Then GoTo 100 Line1 Line1 If newton <> 0 Then If jupdm = 0 Then Msg " METHOD: NEWTON (NO LINE SEARCH)" ElseIf jupdm = 1 Then Msg " METHOD: QUASI-NEWTON (NO LINE SEARCH) USING BROYDEN UPDATE" ElseIf jupdm = 2 Then Msg " METHOD: QUASI-NEWTON (NO LINE SEARCH) USING LEE AND LEE UPDATE" End If Line0 If overch Then Msg " OVERLOW CHECKING IN USE." Else Msg " OVERFLOW CHECKING NOT IN USE." End If If jactyp = 0 Then If njacch > 0 Then s = " ANALYTICAL JACOBIAN USED, CHECKED NUMERICALLY, NJACCH: " + Str$(njacch): Msg s Else Msg " ANALYTICAL JACOBIAN USED; NOT CHECKED." End If ElseIf jactyp = 1 Then Msg " JACOBIAN ESTIMATED USING FORWARD DIFFERENCES." ElseIf jactyp = 2 Then Msg " JACOBIAN ESTIMATED USING BACKWARD DIFFERENCES." Else Msg " JACOBIAN ESTIMATED USING CENTRAL DIFFERENCES." End If Line0 Line1 Else If linesr <> 0 Then Line0 If (deuflh) Then Msg " DEUFLHARD RELAXATION FACTOR INITIALIZATION IN EFFECT." Else Msg " DEUFLHARD RELAXATION FACTOR INITIALIZATION NOT IN EFFECT." End If Else If etafac = one Then Msg " METHOD: TRUST REGION USING SINGLE DOGLEG STEPS." Else Msg " METHOD: TRUST REGION USING DOUBLE DOGLEG STEPS." End If Line0 If cauchy <> 0 Then Msg " INITIAL STEP CONSTRAINED BY SCALED CAUCHY STEP." Else Msg " INITIAL STEP CONSTRAINED BY SCALED NEWTON STEP." End If End If If geoms <> 0 Then Msg " METHOD: GEOMETRIC SEARCH." Else Msg " METHOD: SEARCH BASED ON SUCCESSIVE MINIMIZATIONS." End If If overch <> 0 Then Msg " OVERLOW CHECKING IN USE." Else Msg " OVERFLOW CHECKING NOT IN USE." End If If jupdm = 0 Then _ Msg " NO QUASI-NEWTON UPDATE USED." If jupdm = 1 Then If qnupdm = 0 Then Msg " BROYDEN QUASI-NEWTON UPDATE OF UNFACTORED JACOBIAN." Else Msg " BROYDEN QUASI-NEWTON UPDATE OF FACTORED JACOBIAN." End If ElseIf jupdm = 2 Then If qnupdm = 0 Then Msg " LEE AND LEE QUASI-NEWTON UPDATE OF UNFACTORED JACOBIAN." Else Msg " LEE AND LEE QUASI-NEWTON UPDATE OF FACTORED JACOBIAN." End If End If If jactyp = 0 Then If njacch > 0 Then s = " ANALYTICAL JACOBIAN USED, CHECKED NUMERICALLY, NJACCH: " + Str$(njacch): Msg s Else Msg " ANALYTICAL JACOBIAN USED; NOT CHECKED." End If ElseIf jactyp = 1 Then Msg " JACOBIAN ESTIMATED USING FORWARD DIFFERENCES." ElseIf jactyp = 2 Then Msg " JACOBIAN ESTIMATED USING BACKWARD DIFFERENCES." Else Msg " JACOBIAN ESTIMATED USING CENTRAL DIFFERENCES." End If If linesr = 0 Then Line0 If trupdm = 0 And jupdm > 0 Then Msg " TRUST REGION UPDATED USING POWELL STRATEGY." Else Msg " TRUST REGION UPDATED USING DENNIS AND SCHNABEL STRATEGY." End If End If Line0 Line1 Line0 If itsclf <> 0 Then s = " ADAPTIVE FUNCTION SCALING STARTED AT ITERATION: .........." + Str$(itsclf): Msg s Line0 End If If itsclx <> 0 Then s = " ADAPTIVE VARIABLE SCALING STARTED AT ITERATION: .........." + Str$(itsclx): Msg s Line0 End If If linesr <> 0 Then If jupdm = 0 Then s = " MAXIMUM NUMBER OF STEPS IN LINE SEARCH, MAXNS: ..........." + Str$(maxns): Msg s Else s = " MAXIMUM NUMBER OF NEWTON LINE SEARCH STEPS, MAXNS: ......." + Str$(maxns): Msg s s = " MAXIMUM NUMBER OF QUASI-NEWTON LINE SEARCH STEPS, MAXQNS: " + Str$(maxqns): Msg s End If Else If jupdm = 0 Then s = " MAXIMUM NUMBER OF TRUST REGION UPDATES, MAXNS: ..........." + Str$(maxns): Msg s Else s = " MAXIMUM NO. OF NEWTON TRUST REGION UPDATES, MAXNS: ......." + Str$(maxns): Msg s s = " MAXIMUM NO. OF QUASI-NEWTON TRUST REGION UPDATES, MAXQNS: " + Str$(maxqns): Msg s End If End If If narmij < maxit Then Line0 s = " NUMBER OF OBJECTIVE FUNCTION VALUES COMPARED, MGLL: ......" + Str$(mgll): Msg s End If If jupdm > 0 Then If narmij = maxit Then Line0 s = " MINIMUM NUMBER OF STEPS BETWEEN JACOBIAN UPDATES, MINQNS: " + Str$(minqns): Msg s s = " NUMBER OF NON-QUASI-NEWTON STEPS AT START, NINITN: ......." + Str$(ninitn): Msg s End If s = " NUMBER OF ARMIJO STEPS AT START, NARMIJ: ................." + Str$(narmij): Msg s End If Line0 If stopcr = 3 Then s = " FUNCTION AND STEP SIZE STOPPING CRITERIA, STOPCR: ........" + Str$(stopcr): Msg s ElseIf stopcr = 12 Then s = " FUNCTION OR STEP SIZE STOPPING CRITERIA, STOPCR: ........." + Str$(stopcr): Msg s ElseIf stopcr = 1 Then s = " STEP SIZE STOPPING CRITERION, STOPCR: ...................." + Str$(stopcr): Msg s Else s = " FUNCTION STOPPING CRITERION, STOPCR: ....................." + Str$(stopcr): Msg s End If If newton = 0 Then Line0 If acptcr = 12 Then s = " FUNCTION AND STEP SIZE ACCEPTANCE CRITERIA, ACPTCR: ......" + Str$(acptcr): Msg s ElseIf acptcr = 2 Then s = " STEP SIZE ACCEPTANCE CRITERION, ACPTCR: .................." + Str$(acptcr): Msg s Else s = " FUNCTION ACCEPTANCE CRITERION, ACPTCR: ..................." + Str$(acptcr): Msg s End If If contyp <> 0 Then Line0 s = " CONSTRAINTS IN USE, CONTYP: .............................." + Str$(contyp): Msg s End If End If Line0 Line1 Line0 s = " ESTIMATED MACHINE EPSILON, EPSMCH: ......" + Str$(epsmch): Msg s Line0 s = " FACTOR TO ESTABLISH MAXIMUM STEP SIZE, MSTPF : ......." + Str$(mstpf): Msg s s = " CALCULATED MAXIMUM STEP SIZE, MAXSTP: ................" + Str$(Int(maxstp * 1000) / 1000): Msg s If linesr = 0 Then If delta < zero Then Line0 Msg " INITIAL TRUST REGION NOT PROVIDED." Else Line0 s = " INITIAL TRUST REGION SIZE, DELTA: ...................." + Str$(delta): Msg s End If If etafac < one Then s = " FACTOR TO SET DIRECTION OF TRUST REGION STEP, ETAFAC: ... " + Str$(etafac): Msg s End If s = " TRUST REGION UPDATING FACTOR, DELFAC: ................" + Str$(delfac): Msg s End If If newton = 0 Then Line0 s = " FACTOR IN OBJECTIVE FUNCTION COMPARISON, ALPHA: ......" + Str$(alpha): Msg s If linesr <> 0 And newton = 0 Then Line0 s = " REDUCTION FACTOR FOR RELAXATION FACTOR, SIGMA: ......." + Str$(sigma): Msg s End If If jupdm <> 0 Then Line0 s = " REDUCTION REQUIRED IN OBJ. FUNCTION FOR QN STEP, RATIOF: " + Str$(ratiof): Msg s End If If jupdm = 2 Then Line0 s = " FACTOR IN LEE AND LEE UPDATE, OMEGA: ....................." + Str$(omega): Msg s End If End If Line0 If stopcr <> 2 Then s = " STOPPING TOLERANCE FOR STEP SIZE, STPTOL: ....." + Str$(stptol): Msg s s = " STOPPING TOLERANCE FOR NEWTON STEP, NSTTOL: ..." + Str$(nsttol): Msg s End If If stopcr <> 1 Then s = " STOPPING TOLERANCE FOR OBJECTIVE FUNCTION:....." + Str$(ftol): Msg s End If If linesr <> 0 And newton = 0 And lam0 < one Then Line0 s = " INITIAL LAMBDA IN LINE SEARCH, LAM0: ................." + Str$(lam0): Msg s End If If contyp > 0 Then Line0 s = " FACTOR TO ENSURE STEP WITHIN CONSTRAINTS, CONFAC: ........" + Str$(confac): Msg s End If Line0 Line1 Line0 Msg " SCALING FACTORS" Line0 Msg " COMPONENT VALUES FUNCTION VALUES" Line0 For i = 1 To n s = " SCALEX(" + Str$(i) + ") = " + Str$(scalex(i)) + " SCALEF(" + Str$(i) + ") = " + Str$(scalef(i)): Msg s Next i If contyp > 0 Then Line0 Line1 Line0 Msg " LOWER AND UPPER BOUNDS" Line0 Msg " LOWER BOUNDS UPPER BOUNDS" Line0 For i = 1 To n s = " BOUNDL(" + Str$(i) + ") = " + Str$(boundl(i)) + " BOUNDU(" + Str$(i) + ") = " + Str$(boundu(i)): Msg s Next i End If Line0 100 If output = 2 Then Line1 Line1 Line0 Msg " INITIAL ESTIMATES INITIAL FUNCTION VALUES" Line0 For i = 1 To n s = " X(" + Str$(i) + ") = " + Str$(xc(i)) + " F(" + Str$(i) + ") = " + Str$(fvecc(i)): Msg s Next i Line0 s = " INITIAL OBJECTIVE FUNCTION VALUE = " + Str$(fcnold): Msg s Line0 Line1 Line1 End Sub Sub qmin(nunit, output As Integer, delfts, delta, deltaf, stplen) '----------------------------------------------------------------------- ' FEB. 9, 1991 ' ' SET THE NEW TRUST REGION SIZE, DELTA, BASED ON A QUADRATIC ' MINIMIZATION WHERE DELTA IS THE INDEPENDENT VARIABLE. ' ' DELTAF IS THE DIFFERENCE IN THE SUM-OF-SQUARES OBJECTIVE ' FUNCTION VALUE AND DELFTS IS THE DIRECTIONAL DERIVATIVE IN ' THE DIRECTION OF THE CURRENT STEP, S, WHICH HAS STEP LENGTH STPLEN. '----------------------------------------------------------------------- Const point1 = 0.1 Dim s As String If deltaf - delfts <> zero Then ' CALCULATE DELTA WHERE MINIMUM WOULD OCCUR - DELTMP. ' THIS IS PROVISIONAL AS IT MUST BE WITHIN CERTAIN LIMITS TO BE ACCEPTED deltmp = -delfts * stplen / (two * (deltaf - delfts)) If output > 4 Then Line0 s = " TEMPORARY DELTA FROM QUADRATIC MINIMIZATION: " + Str$(deltmp): Msg s s = " VERSUS CURRENT DELTA: " + Str$(delta): Msg s End If ' REDUCE DELTA DEPENDING ON THE MAGNITUDE OF DELTMP. ' IT MUST BE WITHIN (0.1*DELTA,0.5*DELTA) TO BE ACCEPTED - ' OTHERWISE THE NEAREST ENDPOINT OF THE INTERVAL IS USED If deltmp < point1 * delta Then delta = point1 * delta If output > 4 Then Line0 Msg " NEW DELTA SET TO 0.1 CURRENT DELTA" End If ElseIf deltmp > delta / two Then delta = delta / two If output > 4 Then Line0 Msg " NEW DELTA SET TO 0.5 CURRENT DELTA" End If Else delta = deltmp If output > 4 Then Line0 Msg " NEW DELTA SET TO DELTMP" End If End If Else If output > 4 Then Line0 Msg " TO AVOID OVERFLOW NEW DELTA SET TO 0.5 CURRENT DELTA." End If delta = delta / two End If End Sub Sub trstup(geoms As Integer, newtkn, overch As Integer, overfl As Integer, qrsing As Integer, _ sclfch As Integer, sclxch As Integer, acpcod As Integer, acpstr As Integer, acptcr As Integer, _ contyp As Integer, isejac, jupdm, maxexp, mgll, mnew, n, narmij, nfunc, notrst, _ qnupdm As Integer, retcod As Integer, trupdm As Integer, alpha, confac, _ delfac, delstr, delta, epsmch, fcnmax, fcnnew, fcnold, fcnpre, maxstp As Double, _ newlen As Double, newmax As Double, powtau, rellen, stptol, a(), astore(), _ boundl(), boundu(), delf(), fplpre(), ftrack(), fvec(), fvecc(), hhpi(), jac() As Double, _ rdiag(), rhs(), s(), sbar(), scalef(), scalex(), strack(), xc(), xplpre(), xplus()) '------------------------------------------------------------- ' FEB. 28, 1992 ' ' THIS SUBROUTINE CHECKS FOR ACCEPTANCE OF A TRUST REGION ' STEP GENERATED BY THE DOUBLE DOGLEG METHOD. '------------------------------------------------------------- Const pt5 = 0.5, threeq = 0.75, onept1 = 1.1 Dim convio As Integer Dim V1(ISIZE), wv3(ISIZE) Dim temp1(ISIZE * ISIZE), temp2(ISIZE * ISIZE) Dim s1 As String ' NOTE: ACCEPTANCE CODE, ACPCOD, IS 0 ON ENTRANCE TO TRSTUP convio = 0 overfl = 0 retcod = 0 If output > 3 Then Line0 Line0 If (Not sclfch) And (Not sclxch) Then Msg " TRUST REGION UPDATING" Else Msg " TRUST REGION UPDATING (ALL X''S AND F''S IN UNSCALED UNITS)" End If Line0 End If ' CHECK TO MAKE SURE "S" IS A DESCENT DIRECTION - FIND DIRECTIONAL ' DERIVATIVE AT CURRENT XC USING S GENERATED BY DOGLEG SUBROUTINE. innerp overch, overfl, maxexp, n, n, n, nunit, output, delfts, delf, s If output > 3 Then Line0 s1 = " INNER PRODUCT OF DELF AND S, DELFTS: ........" + Str$(delfts): Msg s1 End If If delfts > zero Then If output > 3 Then Line0 Msg " DIRECTIONAL DERIVATIVE POSITIVE, SEARCH DIRECTION REVERSED." End If For i = 1 To n s(i) = -s(i) Next i End If ' FIND MAXIMUM OBJECTIVE FUNCTION VALUE AND MAXIMIUM STEP ' LENGTH FOR NONMONOTONIC SEARCH. THIS HAS TO BE DONE ONLY ' ONCE DURING EACH ITERATION (WHERE NOTRST=1). If notrst = 1 Then 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 End If ' TEST TRIAL POINT - FIND XPLUS AND TEST FOR CONSTRAINT ' VIOLATIONS IF CONTYP DOES NOT EQUAL 0. For i = 1 To n wv3(i) = -one ' WV3 IS A MARKER FOR "VIOLATORS" - IT CHANGES TO 1 OR 2 xplus(i) = xc(i) + s(i) If contyp > 0 Then If xplus(i) < boundl(i) Then convio = 1 wv3(i) = one ElseIf xplus(i) > boundu(i) Then convio = 1 wv3(i) = two End If End If Next i ' IF CONSTRAINT IS VIOLATED... If convio <> 0 Then If output > 3 Then Line0 Msg " CONSTRAINT VIOLATED:" Line0 Msg " TRIAL ESTIMATES (VIOLATIONS MARKED):" Line0 For i = 1 To n If wv3(i) > zero Then s1 = " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)): Msg s1 Else s1 = " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)): Msg s1 End If Next i End If ' FIND STEP WITHIN CONSTRAINED REGION. ' FIND THE RATIO OF THE DISTANCE FROM THE (I)TH COMPONENT TO ITS ' CONSTRAINT TO THE LENGTH OF THE PROPOSED STEP, XPLUS(I)-XC(I). ' MULTIPLY THIS BY CONFAC (DEFAULT 0.95) TO ENSURE THE NEW STEP STAYS ' WITHIN THE ACCEPTABLE REGION UNLESS XC IS CLOSE TO THE BOUNDARY ' (RATIO <= 1/2). IN SUCH CASES A FACTOR OF 0.5*CONFAC IS USED. ' NOTE: ONLY THE VIOLATING COMPONENTS ARE REDUCED. ratiom = one ' RATIOM STORES THE MINIMUM VALUE OF RATIO For i = 1 To n If wv3(i) = one Then ratio = (boundl(i) - xc(i)) / s(i) ElseIf wv3(i) = two Then ratio = (boundu(i) - xc(i)) / s(i) End If If wv3(i) > zero Then ' NOTE: RATIO IS STORED IN WV3 FOR OUTPUT ONLY wv3(i) = ratio ratiom = XMIN(ratiom, ratio) If ratio > pt5 Then xplus(i) = xc(i) + confac * ratio * s(i) Else ' WITHIN BUFFER ZONE xplus(i) = xc(i) + confac * ratio * s(i) / two End If s(i) = xplus(i) - xc(i) End If Next i If output > 3 Then Line0 Line0 Msg " NEW S AND XPLUS VECTORS (WITH RATIOS FOR VIOLATIONS)" Line0 Msg " NOTE: RATIOS ARE RATIO OF LENGTH TO BOUNDARY FROM CURRENT" Msg " X VECTOR TO MAGNITUDE OF CORRESPONDING PROPOSED STEP." Line0 For i = 1 To n If wv3(i) < zero Then s1 = " S(" + Str$(i) + ") = " + Str$(s(i)) + " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)): Msg s1 Else s1 = " S(" + Str$(i) + ") = " + Str$(s(i)) + " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)) + " " + Str$(wv3(i)) Msg s1 End If Next i Line0 s1 = " MINIMUM OF RATIOS, RATIOM: " + Str$(ratiom): Msg s1 End If ' THE NEW POINT, XPLUS, IS NOT NECESSARILY IN A DESCENT DIRECTION. ' CHECK DIRECTIONAL DERIVATIVE FOR MODIFIED STEP, DLFTSM. innerp overch, overfl, maxexp, n, n, n, nunit, output, dlftsm, delf, s If output > 3 Then Line0 s1 = " INNER PRODUCT OF DELF AND MODIFIED S, DL FTSM: " + Str$(dlftsm): Msg s1 End If ' IF DLFTSM IS POSITIVE REDUCE TRUST REGION. IF NOT, TEST NEW POINT If dlftsm > zero Then delta = confac * ratiom * delta retcod = 8 Exit Sub End If End If 'if convio ' CONSTRAINTS NOT (OR NO LONGER) VIOLATED - TEST NEW POINT fcn1 overfl, n, fvec, xplus nfunc = nfunc + 1 ' IF OVERFLOW AT NEW POINT REDUCE TRUST REGION AND RETURN If overfl <> 0 Then ' IF THE OVERFLOW COMES AS A RESULT OF INCREASING DELTA WITHIN THE ' CURRENT ITERATION (IMPLYING DELSTR IS POSITIVE) AND DIVIDING DELTA ' BY 10 WOULD PRODUCE A DELTA WHICH IS SMALLER THAN THAT AT THE ' STORED POINT, THEN USE STORED POINT AS THE UPDATED RESULT. If delstr > delta / ten Then For i = 1 To n temp1((i - 1) * ISIZE + 1) = xplpre(i): temp2((i - 1) * ISIZE + 1) = xplus(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n xplus(i) = temp2((i - 1) * ISIZE + 1) Next i For i = 1 To n temp1((i - 1) * ISIZE + 1) = fplpre(i): temp2((i - 1) * ISIZE + 1) = fvec(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n fvec(i) = temp2((i - 1) * ISIZE + 1) Next i acpcod = acpstr delta = delstr fcnnew = fcnpre retcod = 1 Else delta = delta / ten retcod = 9 End If Exit Sub End If ' NO OVERFLOW IN RESIDUAL VECTOR If output > 3 Then Line0 Msg " TRIAL ESTIMATES FUNCTION VALUES" Line0 For i = 1 To n s1 = " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)) + " FVEC(" + Str$(i) + ") = " + Str$(fvec(i)): Msg s1 Next i End If ' IF NO OVERFLOW WITHIN RESIDUAL VECTOR FIND OBJECTIVE FUNCTION fcnevl overfl, maxexp, n, nunit, output, epsmch, fcnnew, fvec, scalef ' IF OVERFLOW IN OBJECTIVE FUNCTION EVALUATION REDUCE TRUST REGION AND RETURN If overfl <> 0 Then ' IF THE OVERFLOW COMES AS A RESULT OF INCREASING DELTA WITHIN THE ' CURRENT ITERATION (SO THAT DELSTR IS POSITIVE) AND DIVIDING DELTA ' BY 10 WOULD PRODUCE A DELTA WHICH IS SMALLER THAN THAT AT THE STORED ' POINT THEN USE STORED POINT AS THE UPDATED RESULT. If delstr > delta / ten Then For i = 1 To n temp1((i - 1) * ISIZE + 1) = xplpre(i): temp2((i - 1) * ISIZE + 1) = xplus(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n xplus(i) = temp2((i - 1) * ISIZE + 1) Next i For i = 1 To n temp1((i - 1) * ISIZE + 1) = fplpre(i): temp2((i - 1) * ISIZE + 1) = fvec(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n fvec(i) = temp2((i - 1) * ISIZE + 1) Next i acpcod = acpstr delta = delstr fcnnew = fcnpre retcod = 2 Else delta = delta / ten retcod = 10 End If Exit Sub Else ' NO OVERFLOW AT TRIAL POINT - COMPARE OBJECTIVE FUNCTION TO FCNMAX If output > 3 Then Line0 If (Not sclfch) Then s1 = " OBJECTIVE FUNCTION AT XPLUS, FCNNEW: ........." + Str$(fcnnew): Msg s1 Else s1 = " SCALED OBJECTIVE FUNCTION AT XPLUS, FCNNEW: .." + Str$(fcnnew): Msg s1 End If s1 = " COMPARE TO FCNMAX+ALPHA*DELFTS: .............." + Str$(fcnmax + alpha * delfts): Msg s1 End If End If ' IF ACPTCR=12 CHECK SECOND DEUFLHARD STEP ACCEPTANCE TEST ' BY FINDING 2-NORM OF SBAR. THERE ARE FOUR POSSIBILITIES ' DEPENDING ON WHETHER THE JACOBIAN IS OR IS NOT SINGULAR ' AND WHETHER QNUPDM IS 0 OR 1 If acptcr = 12 Then If (qrsing) Then ' FORM -J^F AS RIGHT HAND SIDE - METHOD DEPENDS ON ' WHETHER QNUPDM EQUALS 0 OR 1 (UNFACTORED OR FACTORED) If qnupdm = 0 Then ' UNSCALED JACOBIAN IS IN MATRIX JAC For k = 1 To n wv3(k) = -fvec(k) * scalef(k) * scalef(k) Next k mmulv n, jac, wv3, rhs Else ' R IN UPPER TRIANGLE OF A PLUS RDIAG AND Q^ IN JAC ' FROM QR DECOMPOSITION OF SCALED JACOBIAN. For i = 1 To n Sum = zero For j = 1 To n Sum = Sum - jac((i - 1) * ISIZE + j) * fvec(j) * scalef(j) Next j wv3(i) = Sum Next i rhs(1) = rdiag(1) * wv3(1) For j = 2 To n For k = 1 To j - 1 V1(k) = a((k - 1) * ISIZE + j) Next k Sum = DotProduct(j - 1, V1, wv3) rhs(j) = Sum + rdiag(j) * wv3(j) Next j End If chsolv overch, overfl, maxexp, n, nunit, output, a, rhs, sbar Else ' RIGHT HAND SIDE IS -FVEC If qnupdm = 0 Or jupdm = 0 Then ' QR DECOMPOSITION OF SCALED JACOBIAN STORED IN ASTORE For k = 1 To n sbar(k) = -fvec(k) * scalef(k) Next k qrsolv overch, overfl, maxexp, n, nunit, output, astore, hhpi, rdiag, sbar Else ' SET UP RIGHT HAND SIDE - MULTIPLY -FVEC BY Q^ ' (STORED IN JAC). RHS IS A WORK VECTOR ONLY HERE For k = 1 To n wv3(k) = -fvec(k) * scalef(k) Next k avmul n, n, n, n, jac, wv3, sbar rsolv overch, overfl, maxexp, n, nunit, output, a, rdiag, sbar End If End If ' NORM OF (SCALED) SBAR IS NEEDED FOR SECOND ACCEPTANCE TEST For k = 1 To n wv3(k) = scalex(k) * sbar(k) Next k twonrm overfl, maxexp, n, epsmch, sbrnrm, wv3 If output > 4 Then Line0 If sclxch = 0 Then Msg " DEUFLHARD SBAR VECTOR" Line0 For i = 1 To n s1 = "sbar(" + Str$(i) + ") = " + Str$(sbar(i)): Msg s1 Next i Else Msg " DEUFLHARD SBAR VECTOR IN SCALED X UNITS" Line0 For i = 1 To n s1 = " SBAR(" + Str$(i) + ") = " + Str$(sbar(i)) + " SBAR('+Str$(i)" + ") = " _ + Str$(scalex(i) * sbar(i)): Msg s1 Next i End If End If If output > 3 Then Line0 If sclxch = 0 Then s1 = " VALUE OF SBRNRM AT XPLUS: ................." + Str$(sbrnrm): Msg s1 Else s1 = " VALUE OF SCALED SBRNRM AT XPLUS: .........." + Str$(sbrnrm): Msg s1 End If s1 = " NEWMAX: ..................................." + Str$(ewmax): Msg s1 End If If sbrnrm < newmax Then acpcod = 2 ' FUNCTION VALUE ACCEPTANCE IS ALSO CHECKED REGARDLESS ' OF WHETHER SECOND STEP ACCEPTANCE CRITERION WAS MET. End If 'if acptor=12 ' ESTABLISH DELTAF FOR USE IN COMPARISON TO PREDICTED ' CHANGE IN OBJECTIVE FUNCTION, DELFPR, LATER deltaf = fcnnew - fcnold If fcnnew >= fcnmax + alpha * delfts Then ' FAILURE OF FIRST STEP ACCEPTANCE TEST. TEST LENGTH OF ' STEP TO ENSURE PROGRESS IS STILL BEING MADE rellen = zero For i = 1 To n rellen = XMAX(rellen, Abs(s(i)) / XMAX((Abs(xplus(i))), one / scalex(i))) Next i If rellen < stptol Then ' NO PROGRESS BEING MADE - RETCOD = 7 STOPS PROGRAM For i = 1 To n temp1((i - 1) * ISIZE + 1) = xc(i): temp2((i - 1) * ISIZE + 1) = xplus(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n xplus(i) = temp2((i - 1) * ISIZE + 1) Next i retcod = 7 Exit Sub Else ' FAILURE OF STEP BY OBJECTIVE FUNCTION CRITERION. ' ESTABLISH A NEW DELTA FROM EITHER SIMPLE DIVISION ' BY DELFAC OR BY FINDING THE MINIMUM OF A QUADRATIC MODEL If (geoms) Then delta = delta / delfac Else ' FIRST FIND LENGTH OF TRUST REGION STEP For i = 1 To n wv3(i) = s(i) * scalex(i) Next i twonrm overfl, maxexp, n, epsmch, stplen, wv3 qmin nunit, output, delfts, delta, deltaf, stplen End If If delta < delstr Then ' IF DELTA HAS BEEN INCREASED AT THIS ITERATION AND THE ' DELTA FROM QMIN IS LESS THAN THE DELTA AT THE PREVIOUSLY ' ACCEPTED (STORED) POINT THEN RETURN TO THAT POINT AND ' ACCEPT IT AS THE UPDATED ITERATE. For i = 1 To n temp1((i - 1) * ISIZE + 1) = xplpre(i): temp2((i - 1) * ISIZE + 1) = xplus(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n xplus(i) = temp2((i - 1) * ISIZE + 1) Next i For i = 1 To n temp1((i - 1) * ISIZE + 1) = fplpre(i): temp2((i - 1) * ISIZE + 1) = fvec(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n fvec(i) = temp2((i - 1) * ISIZE + 1) Next i acpcod = acpstr delta = delstr fcnnew = fcnpre retcod = 3 Exit Sub End If ' IF THE SECOND ACCEPTANCE TEST HAS BEEN PASSED RETURN ' WITH NEW TRUST REGION AND CONTINUE ON TO NEXT ITERATION; ' OTHERWISE TRY A NEW STEP WITH REDUCED DELTA. If acpcod = 2 Then retcod = 4 Else ' FAILURE OF FIRST STEP ACCEPTANCE TEST retcod = 11 End If Exit Sub End If Else ' OBJECTIVE FUNCTION MEETS FIRST ACCEPTANCE CRITERION. ' IN NONMONOTONIC SEARCHES IT MAY BE GREATER THAN THE ' PREVIOUS OBJECTIVE FUNCTION VALUE - CONSIDER THIS CASE FIRST. If deltaf >= alpha * delfts Then ' AN ACCEPTABLE STEP HAS BEEN FOUND FOR THE NONMONOTONIC SEARCH ' BUT THE OBJECTIVE FUNCTION VALUE IS NOT A "DECREASE" FROM THE ' PREVIOUS ITERATION (ACTUALLY IT MIGHT BE BETWEEN ZERO AND ' ALPHA*DELFTS). ACCEPT STEP BUT REDUCE DELTA. delta = delta / delfac retcod = 5 If acpcod = 2 Then acpcod = 12 Else acpcod = 1 End If Exit Sub End If ' COMPARE DELTAF TO DELTAF PREDICTED, DELFPR, TO DETERMINE NEXT ' TRUST REGION SIZE. NOTE: DELTAF MUST BE LESS THAN ALPHA*DELFTS ' (IN ESSENCE NEGATIVE) TO HAVE REACHED THIS POINT IN TRSTUP. ' R IS IN UPPER TRIANGLE OF MATRIX A SO THE FOLLOWING CODE FINDS: ' DELFPR = DELF^S + 1/2 S^J^JS = DELF^S + 1/2 S^R^RS uvmul n, n, n, n, a, s, wv3 delfpr = delfts + DotProduct(n, wv3, wv3) / two If output > 4 Then Line0 s1 = " PREDICTED CHANGE IN OBJECTIVE FUNCTION, DELFPR: " + Str$(delfpr): Msg s1 s1 = " ACTUAL CHANGE IN OBJECTIVE FUNCTION, DELTAF: " + Str$(deltaf): Msg s1 End If If (retcod <= 6 And Abs(delfpr - deltaf) <= Abs(deltaf) / ten) Or (deltaf <= delfts And newtkn = 0) And convio = 0 _ And delstr = zero Then If XMIN(newlen, maxstp) / delta > onept1 Then ' PROMISING STEP - INCREASE TRUST REGION. ' STORE CURRENT POINT. For i = 1 To n temp1((i - 1) * ISIZE + 1) = xplus(i): temp2((i - 1) * ISIZE + 1) = xplpre(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n xplpre(i) = temp2((i - 1) * ISIZE + 1) Next i For i = 1 To n temp1((i - 1) * ISIZE + 1) = fvec(i): temp2((i - 1) * ISIZE + 1) = fplpre(i) Next i matcop n, n, 1, 1, n, 1, temp1, temp2 For i = 1 To n fplpre(i) = temp2((i - 1) * ISIZE + 1) Next i delstr = delta fcnpre = fcnnew ' IF NONMONOTONIC STEPS ARE BEING USED EXPAND TRUST ' REGION TO NEWLEN, OTHERWISE EXPAND BY DELFAC. If isejac > narmij Then delta = XMIN(newlen, maxstp) Else delta = XMIN(delfac * delta, maxstp) End If retcod = 12 If acpcod = 2 Then acpstr = 12 Else acpstr = 1 End If acpcod = 0 Else retcod = 0 If acpcod = 2 Then acpcod = 12 Else acpcod = 1 End If End If Exit Sub Else ' CHANGE TRUST REGION SIZE DEPENDING ON DELTAF AND DELFPR retcod = 6 If acpcod = 2 Then acpcod = 12 Else acpcod = 1 End If If deltaf >= delfpr / ten Then delta = delta / delfac If output > 3 Then Line0 Msg " CHANGE IN F, DELTAF, IS > 0.1 DELFPR - REDUCE DELTA." End If ElseIf (trupdm = 0) And (jupdm > 0) Then ' POWELL'S UPDATING SCHEME - FIND JAC S FIRST If qnupdm = 0 Then ' UNSCALED JACOBIAN IN JAC For i = 1 To n rhs(i) = s(i) * scalef(i) Next i avmul n, n, n, n, jac, rhs, wv3 Else ' MULTIPLY BY R FIRST uvmul n, n, n, n, a, s, rhs ' THEN Q (IN JAC^) mmulv n, jac, rhs, wv3 End If dmult = delfpr / ten - deltaf sp = zero ss = zero For k = 1 To n wv3(k) = wv3(k) + fvecc(k) sp = sp + Abs(fvec(k) * (fvec(k) - wv3(k))) ss = ss + (fvec(k) - wv3(k)) * (fvec(k) - wv3(k)) Next k If sp + Sqr(sp * sp + dmult * ss) < epsmch Then powlam = ten Else powlam = one + dmult / (sp + Sqr(sp * sp + dmult * ss)) End If powlam = Sqr(powlam) powmu = XMIN(XMIN(delfac, powlam), powtau) powtau = powlam / powmu If output > 3 Then Line0 Msg " FACTORS IN POWELL UPDATING SCHEME:" Line0 s1 = " LAMBDA: " + Str$(powlam) + " MU: " + Str$(powmu) + " TAU: " + Str$(powtau): Msg s1 Msg " DELTA IS MINIMUM OF MU*DELTA AND MAXSTP." End If delta = XMIN(powmu * delta, maxstp) Else If deltaf < threeq * delfpr Then delta = XMIN(delfac * delta, maxstp) If output > 3 Then Line0 Msg " CHANGE IN F, DELTAF, IS LESS THAN 0.75 X PREDICTED." s1 = " DELTA INCREASED TO: " + Str$(delta): Msg s1 End If Else If output > 3 Then Line0 Msg " DELTAF BETWEEN 0.1 AND 0.75 DELFPR - LEAVE DELTA UNCHANGED" End If End If End If End If End If 'if fcnnew... End Sub Sub rcdprt(nunit, retcod As Integer, delta, rellen, stptol) '------------------------------------------------------------------------- ' FEB. 14, 1991 ' ' DESCRIBE MEANING OF RETURN CODES, RETCOD, FROM TRUST REGION UPDATING. '------------------------------------------------------------------------- Dim s As String Line0 s = " RETCOD, FROM TRUST REGION UPDATING: " + Str$(retcod): Msg s Line0 If retcod = 1 Then Msg " PROMISING STEP FOUND; DELTA HAS BEEN INCREASED TO NEWLEN BUT" Msg " BECAUSE OF OVERFLOWS IN THE FUNCTION VECTOR(S) IN SUBSEQUENT" Msg " STEP(S) THE PROJECTED DELTA IS LESS THAN THAT AT THE ALREADY" Msg " SUCCESSFUL STEP - RETURN TO SUCCESSFUL STEP AND ACCEPT AS" Msg " NEW POINT." ElseIf retcod = 2 Then Msg " PROMISING STEP FOUND; DELTA HAS BEEN INCREASED TO NEWLEN BUT" Msg " BECAUSE OF OVERFLOWS IN THE OBJECTIVE FUNCTION IN SUBSEQUENT" Msg " STEP(S) THE PROJECTED DELTA IS LESS THAN THAT AT THE ALREADY" Msg " SUCCESSFUL STEP - RETURN TO SUCCESSFUL STEP AND ACCEPT AS" Msg " NEW POINT.')" ElseIf retcod = 3 Then Msg " PROMISING STEP FOUND; DELTA HAS BEEN INCREASED TO NEWLEN BUT" Msg " BECAUSE OF SUBSEQUENT FAILURES IN THE STEP ACCEPTANCE TEST(S)" Msg " THE PROJECTED DELTA IS LESS THAN THAT AT THE ALREADY" Msg " SUCCESSFUL STEP - RETURN TO SUCCESSFUL STEP AND ACCEPT." ElseIf retcod = 4 Then Msg " STEP ACCEPTED BY STEP SIZE CRITERION ONLY - DELTA REDUCED" ElseIf retcod = 5 Then Msg " STEP ACCEPTED - NEW FUNCTION VALUE > PREVIOUS =>" Msg " REDUCE TRUST REGION." ElseIf retcod = 6 Then Msg " STEP ACCEPTED - DELTA CHANGED AS DETAILED ABOVE." ElseIf retcod = 7 Then Msg " NO PROGRESS MADE: RELATIVE STEP SIZE IS TOO SMALL" s = " REL. STEP SIZE, RELLEN = " + Str$(rellen) + ", STPTOL = " + Str$(stptol): Msg s ElseIf retcod = 8 Then Msg " POINT MODIFIED BY CONSTRAINTS NOT A DESCENT DIRECTION" Msg " DELTA REDUCED TO CONFAC*RATIOM*DELTA." ElseIf retcod = 9 Then Msg " OVERFLOW DETECTED IN FUNCTION VECTOR - DELTA REDUCED." ElseIf retcod = 10 Then Msg " OVERFLOW IN OBJECTIVE FUNCTION - DELTA REDUCED." ElseIf retcod = 11 Then Msg " STEP NOT ACCEPTED - REDUCE TRUST REGION SIZE BY MINIMIZATION" Msg " OF QUADRATIC MODEL IN STEP DIRECTION." Else Msg " PROMISING STEP - INCREASE DELTA TO NEWLEN AND TRY A NEW STEP." End If Line0 s = " DELTA ON RETURN FROM TRUST REGION UPDATING: " + Str$(delta): Msg s End Sub Sub nnes(absnew As Integer, cauchy As Integer, deuflh As Integer, geoms As Integer, linesr, newton, _ overch As Integer, acptcr As Integer, itsclf, itsclx, jupdm, maxns, maxqns, minqns, n, narmij, _ niejev, njacch, qnupdm As Integer, stopcr As Integer, supprs As Integer, trupdm As Integer, _ alpha, confac, delta, delfac, etafac, fcnnew, fdtolj, ftol, mstpf As Double, _ nsttol As Double, omega, ratiof, sigma, stptol, a(), boundl(), boundu(), delf(), fsave(), ftrack(), _ fvec(), fvecc(), h(), hhpi(), jac() As Double, plee(), rdiag(), s(), sbar(), scalef(), scalex(), _ sn(), ssdhat(), strack(), vhat(), xc(), xplus(), xsave()) '---------------------------------------------------------------------------- ' FEB. 28, 1992 ' *** Solve a non-linear system - Main procedure *** '---------------------------------------------------------------------------- 'Labels: 310, 320, 330 Const pt01 = 0.01 Dim maxstp As Double Dim newlen As Double Dim newmax As Double Dim wv1(ISIZE), wv2(ISIZE), wv4(ISIZE) Dim acpcod As Integer Dim acpstr As Integer Dim contyp As Integer Dim countr As Integer Dim outtmp As Integer Dim retcod As Integer Dim abort As Integer Dim checkj As Integer Dim frstdg As Integer Dim overfl As Integer Dim qnfail As Integer Dim qrsing As Integer Dim restrt As Integer Dim savest As Integer Dim sclfch As Integer Dim sclxch As Integer Dim ss As String newstm = 0: mgll = 10 ' PRINT HELP IF REQUESTED If help <> "NONE" Then ' olhelp nunit, help Exit Sub End If qnfail = 0 restrt = 1 savest = 0 countr = 0 isejac = 0 mnew = 0 nac1 = 0 nac2 = 0 nac12 = 0 nfunc = 1 njetot = 0 outtmp = output trmcod = 0 retcod = 0 ' ESTABLISH INITIAL FUNCTION VALUE AND CHECK FOR STARTING ESTIMATE WHICH ' IS A SOLUTION. ALSO, CHECK FOR INCOMPATIBILITIES IN INPUT PARAMETERS. initch instop, linesr, newton, overfl, sclfch, sclxch, acptcr, contyp, _ jactyp, jupdm, maxexp, n, nunit, output, qnupdm, stopcr, _ trupdm, epsmch, fcnold, ftol, boundl, boundu, fvecc, scalef, _ scalex, xc ' IF A FATAL ERROR IS DETECTED IN INITCH RETURN TO MAIN PROGRAM. ' NOTE: SOME INCOMPATIBILITIES ARE CORRECTED WITHIN INITCH AND EXECUTION ' CONTINUES. WARNINGS ARE GENERATED WITHIN INITCH. If instop <> 0 Then Exit Sub ' ESTABLISH MAXIMUM STEP LENGTH ALLOWED (USUALLY THIS IS MUCH ' LARGER THAN ACTUAL STEP SIZES - IT IS ONLY TO PREVENT ' EXCESSIVELY LARGE STEPS). THE FACTOR MSTPF CONTROLS THE ' MAGNITUDE OF MAXSTP AND IS SET BY THE USER (DEFAULT=1000). maxst overfl, maxexp, n, nunit, output, epsmch, maxstp, mstpf, scalex, xc ' WRITE TITLE AND RECORD PARAMETERS title cauchy, deuflh, geoms, linesr, newton, overch, acptcr, contyp, _ itsclf, itsclx, jactyp, jupdm, maxit, maxns, maxqns, mgll, _ minqns, n, narmij, niejev, njacch, nunit, output, qnupdm, stopcr, _ trupdm, alpha, confac, delfac, delta, epsmch, etafac, fcnold, _ ftol, lam0, maxstp, mstpf, nsttol, omega, ratiof, sigma, stptol, _ boundl, boundu, fvecc, scalef, scalex, xc ' INITIALIZE FTRACK AND STRACK VECTORS (FTRACK STORES (TRACKS) THE FUNCTION ' VALUES FOR THE NONMONOTONIC COMPARISON AND STRACK, SIMILARLY, STORES THE ' LENGTH OF THE NEWTON STEPS - WHICH ARE USED IN CONJUNCTION WITH ' DEUFLHARD 'S SECOND ACCEPTANCE CRITERION). For j = 0 To mgll - 1 ftrack(j) = zero strack(j) = zero Next j ' MAIN ITERATIVE LOOP - MAXIT IS SPECIFIED BY USER. ' ITNUM COUNTS OVERALL ITERATIONS. ' ISEJAC COUNTS ITERATIONS SINCE LAST EXPLICIT JACOBIAN ' EVALUATION IF A QUASI-NEWTON METHOD IS BEING USED. For itnum = 1 To maxit ' SUPPRESS OUTPUT IF DESIRED - USED IF DETAILED OUTPUT IS DESIRED ' FOR LATER ITERATIONS ONLY (DEFAULT VALUE FOR SUPPRS IS 0 - ' I.E. NO SUPPRESSION). If itnum < supprs Then output = 3 Else output = outtmp End If If output > 2 Then Line1 Line0 ss = " ITERATION NUMBER: " + Str$(itnum): Msg ss End If ' UPDATE ITERATION COUNTER, ISEJAC, FOR QUASI-NEWTON METHODS ' ONLY (I.E. JUPDM > 0). If jupdm > 0 And restrt <> 0 And newton = 0 Then isejac = 1 Else isejac = isejac + 1 End If If output > 4 And jupdm > 0 And newton = 0 Then Line0 If restrt <> 0 Then If itnum > niejev Then Msg " RESTRT IS TRUE, ISEJAC SET TO 1." Else ss = " # OF ITERATIONS SINCE EXPLICIT JACOBIAN, ISEJAC, INCREASED TO " + Str$(isejac): Msg ss End If End If ' WHEN AN EXPLICIT JACOBIAN IS BEING USED IN QUASI-NEWTON METHODS THEN ' MAXNS STEPS ARE ALLOWED IN THE LINE SEARCH. FOR STEPS BASED ON A ' QUASI-NEWTON APPROXIMATION MAXQNS ARE ALLOWED. SIMILARLY MAXNS AND ' MAXQNS STEPS ARE ALLOWED IN TRUST REGION METHODS RESPECTIVELY. ' THIS IS AN ATTEMPT TO AVOID AN EXCESSIVE NUMBER OF FUNCTION ' EVALUATIONS IN A DIRECTION WHICH WOULD NOT LEAD TO A SIGNIFICANT ' REDUCTION. If newton = 0 Then If linesr <> 0 Then If jupdm > 0 Then If isejac = 1 Then ' JACOBIAN UPDATED EXPLICITLY maxlin = maxns Else ' QUASI-NEWTON UPDATE maxlin = maxqns End If Else maxlin = maxns End If Else If jupdm > 0 Then If isejac = 1 Then ' JACOBIAN UPDATED EXPLICITLY maxtrs = maxns Else ' QUASI-NEWTON UPDATE maxtrs = maxqns End If Else maxtrs = maxns End If End If If jupdm > 0 And output > 4 Then Line0 If (linesr) Then ss = " MAXLIN SET TO: " + Str$(maxlin): Msg ss Else ss = " MAXTRS SET TO: " + Str$(maxtrs): Msg ss End If End If End If ' ESTABLISH WHETHER JACOBIAN IS TO BE CHECKED NUMERICALLY - ' NJACCH ESTABLISHES THE NUMBER OF ITERATIONS FOR WHICH JACOBIAN ' CHECKING IS DESIRED. IF CHECKJ IS TRUE, A FORWARD DIFFERENCE NUMERICAL ' APPROXIMATION OF THE JACOBIAN IS COMPARED TO THE ANALYTICAL VERSION. ' STORE THE NUMBER OF FUNCTION EVALUATIONS SO THAT THESE "EXTRA" ARE NOT ' INCLUDED IN OVERALL STATISTICS. If jactyp = 0 Then If itnum > njacch Then checkj = 0 Else checkj = 1 nfestr = nfetot End If Else checkj = 0 End If ' EVALUATE JACOBIAN AT FIRST STEP OR IF NO QUASI-NEWTON ' UPDATE IS BEING USED (RESTRT IS FALSE ONLY IN QUASI- ' NEWTON METHODS WHEN THE QUASI-NEWTON UPDATE IS BEING USED). If restrt <> 0 Or jupdm = 0 Then ' IF MORE THAN ONE DAMPED NEWTON STEP HAS BEEN REQUESTED AT THE START, ' IDENTIFY THIS AS THE REASON FOR EXPLICIT JACOBIAN EVALUATION. If jupdm > 0 And itnum <= niejev And itnum > 1 And output > 4 Then Line0 Msg " AS ITNUM <= NIEJEV JACOBIAN EVALUATED EXPLICITLY." End If ' OTHERWISE JUPDM IS 0 OR RESTRT IS TRUE AND ITNUM IS > NIEJEV If jupdm > 0 And itnum > 1 And output > 4 And itnum > niejev Then Line0 Msg " RESTRT IS TRUE - JACOBIAN EVALUATED EXPLICITLY." End If ' NOTE: MATRIX H IS USED HERE TO HOLD THE FINITE DIFFERENCE ' ESTIMATION USED IN CHECKING THE ANALYTICAL JACOBIAN IF CHECKJ ' IS TRUE. VECTORS WV1 AND, FOR CENTRAL DIFFERENCES, WV2 ' TEMPORARILY HOLD THE FINITE DIFFERENCE FUNCTION EVALUATIONS. jacobi checkj, jacerr, overfl, jactyp, n, nunit, output, epsmch, _ fdtolj, boundl, boundu, fvecc, wv1, wv2, jac, h, scalex, xc njetot = njetot + 1 ' RETURN IF ANALYTICAL AND NUMERICAL JACOBIANS DON'T AGREE ' (APPLICABLE ONLY IF CHECKJ IS TRUE AND A DISCREPANCY IS FOUND ' WITHIN SUBROUTINE JACOBI). A WARNING IS GIVEN FROM WITHIN JACOBI If jacerr <> 0 Then Exit Sub ' RESET TOTAL NUMBER OF FUNCTION EVALUATIONS TO NEGLECT ' THOSE USED IN CHECKING ANALYTICAL JACOBIAN. If checkj <> 0 Then nfetot = nfestr If jupdm > 0 Then ' FCNMIN IS THE MINIMUM OF THE OBJECTIVE FUNCTION FOUND SINCE ' THE LAST EXPLICIT JACOBIAN EVALUATION. IT IS USED TO ESTABLISH ' WHICH STEP THE PROGRAM RETURNS TO WHEN A QUASI-NEWTON STEP FAILS. fcnmin = fcnold ' POWTAU IS THE TAU FROM POWELL'S TRUST REGION UPDATING SCHEME ' (USED WHEN JUPDM=1). IT IS RESET TO 1.0 AT EVERY EXPLICIT ' JACOBIAN EVALUATION IN QUASI-NEWTON METHODS. powtau = one ' AT EVERY EXPLICIT JACOBIAN EVALUATION EXCEPT POSSIBLY THE FIRST, ' IN QUASI-NEWTON METHODS, A NEW TRUST REGION IS CALCULATED ' INTERNALLY USING EITHER THE NEWTON STEP OR THE CAUCHY STEP AS ' SPECIFIED BY THE USER IN LOGICAL VARIABLE "CAUCHY" IN ' SUBROUTINE DELCAU. THIS IS FORCED BY SETTING DELTA TO A ' NEGATIVE NUMBER. If itnum > 1 Then delta = -one ' RESET COUNTER FOR NUMBER OF FAILURES OF RATIO TEST (IS ' FCNNEW/FCNOLD < RATIOF?) SINCE LAST EXPLICIT JACOBIAN UPDATE nfail = 0 If output > 4 And itnum > niejev Then Line0 Msg " NUMBER OF FAILURES OF RATIO TEST, NFAIL, SET BACK TO 0." End If ' ESTABLISH A NEW MAXIMUM STEP LENGTH ALLOWED If itnum > 1 Then maxst overfl, maxexp, n, nunit, output, epsmch, maxstp, mstpf, scalex, xc ' SET "P" MATRIX TO IDENTITY FOR LEE AND LEE UPDATE (JUPDM=2). If jupdm = 2 And itnum >= niejev Then For j = 1 To n For i = 1 To n plee((i - 1) * ISIZE + j) = zero Next i plee((j - 1) * ISIZE + j) = one Next j End If End If End If If (restrt <> 0 Or qnupdm = 0) And output > 4 And matsup = 0 Then ' WRITE JACOBIAN MATRIX Line0 Msg " JACOBIAN MATRIX:" matprt n, n, jac End If ' ESTABLISH SCALING MATRICES IF DESIRED (ITSCLF=0 => NO ADAPTIVE ' SCALING, WHILE ITSCLF > 0 MEANS ADAPTIVE SCALING STARTS AFTER THE ' (ITSCLF]TH ITERATION). SIMILARLY FOR COMPONENT SCALING... ' NOTE: SCALING FACTORS ARE UPDATED ONLY WHEN THE JACOBIAN IS UPDATED ' EXPLICITLY IN QUASI-NEWTON METHODS. ' FUNCTION SCALING. If restrt <> 0 And itsclf > 0 And itsclf <= itnum Then ascalf n, epsmch, fvecc, jac, scalef If output > 4 Then Line0 Msg " FUNCTION SCALING VECTOR:" Line0 For i = 1 To n ss = " SCALEF(" + Str$(i) + ") = " + Str$(scalef(i)): Msg ss Next i End If ' RECALCULATE OBJECTIVE FUNCTION WITH NEW SCALING FACTORS. ' THIS AVOIDS PREMATURE FAILURES IF THE CHANGE IS SCALING FACTORS ' WOULD MAKE THE PREVIOUS OBJECTIVE FUNCTION VALUE SMALLER. fcnevl overfl, maxexp, n, nunit, output, epsmch, fcnold, fvecc, scalef End If ' COMPONENT SCALING If jupdm = 0 Or qnupdm = 0 Then ' NEWTON STEP - UNFACTORED FORM BEING UPDATED nstpun abort, linesr, overch, overfl, qrsing, sclfch, sclxch, _ itnum, maxexp, n, nunit, output, epsmch, a, delf, fvecc, h, _ hhpi, jac, rdiag, scalef, scalex, sn Else ' NEWTON STEP - FACTORED FORM BEING UPDATED nstpfa abort, linesr, overch, overfl, qrsing, restrt, sclfch, _ sclxch, itnum, maxexp, n, newstm, nunit, output, epsmch, a, _ delf, fvecc, h, hhpi, jac, rdiag, scalef, scalex, sn End If ' RUN IS ABORTED IF THE JACOBIAN BECOMES ESSENTIALLY ALL ZEROS. If abort <> 0 Then Exit Sub ' CHECK FOR CONVERGENCE ON LENGTH OF NEWTON STEP IF ' STOPCR = 1, 12 OR 3. If stopcr <> 2 Then stpmax = zero For i = 1 To n stpmax = XMAX(stpmax, Abs(sn(i)) / XMAX(xc(i), scalex(i))) wv1(i) = scalex(i) * xc(i) Next i twonrm overfl, maxexp, n, epsmch, xnorm, wv1 If stpmax <= nsttol * (one + xnorm) Then trmcod = 1 ' IF STOPCR=3 THEN OBJECTIVE FUNCTION VALUE MUST BE ' DETERMINED AS WELL - OTHERWISE A SOLUTION HAS BEEN FOUND. If stopcr <> 3 Then GoTo 330 ' NOTE: STATEMENT 202 PRECEDES CONVERGENCE CHECKING SUBROUTINE End If End If ' FIND LENGTH OF (SCALED) NEWTON STEP, NEWLEN For i = 1 To n wv1(i) = scalex(i) * sn(i) Next i twonrm overfl, maxexp, n, epsmch, newlen, wv1 ' FOR ITERATIONS AFTER THE ARMIJO STEPS HAVE BEEN COMPLETED ' AT THE BEGINNING (IN OTHER WORDS THE MONOTONIC STEPS) ' STORE THE FUNCTION AND, POSSIBLY, THE NEWTON STEP LENGTHS ' IN THE FTRACK AND STRACK VECTORS, RESPECTIVELY. If isejac >= narmij Then If isejac = 1 Then strack(0) = newlen ' NEWMAX IS USED TO KEEP A BOUND ON THE ENTRIES IN THE STRACK VECTOR. newmax = newlen Else strack(countr) = XMIN(newmax, newlen) End If ' THE OBJECTIVE FUNCTION VALUE IS STORED EVEN IF IT IS ' GREATER THAN ANY PRECEEDING FUNCTION VALUE. ftrack(countr) = fcnold ' WRITE FTRACK AND STRACK VECTORS IF DESIRED. SINCE ONLY THE LAST ' MGLL VALUES ARE NEEDED THE COUNTER CIRCULATES THROUGH THE VECTOR ' CAUSING ONLY THE MGLL MOST RECENT VALUES TO BE KEPT. ' NOTE: THESE VECTORS ARE NOT APPLICABLE IF NEWTON'S METHOD IS BEING ' USED. If newton = 0 And output > 4 Then Line0 Line0 ' IF ONLY THE FUNCTION VALUE ACCEPTANCE TEST IS BEING USED, ' THUS ACPTCR=1, THEN ONLY THE FTRACK VECTOR IS APPLICABLE. If acptcr = 1 Then ss = " CURRENT FTRACK VECTOR; LATEST CHANGE: ELEMENT " + Str$(countr): Msg ss Line0 For j = 0 To mgll - 1 ss = " FTRACK(" + Str$(j) + ") = " + Str$(ftrack(j)): Msg ss Next j Else ' BOTH THE FUNCTION VALUE AND THE STEP SIZE ' ACCEPTANCE TESTS ARE BEING USED, ACPTCR=12. ss = " CURRENT FTRACK AND STRACK VECTORS; LATEST CHANGE: ELEMENT " + Str$(countr): Msg ss Line0 For j = 0 To mgll - 1 ss = " FTRACK(" + Str$(j) + ") = " + Str$(ftrack(j)) + " STRACK(" + Str$(j) + ") = " + Str$(strack(j)) Msg ss Next j End If End If ' UPDATE COUNTING INTEGER, COUNTR. RECYCLE IF COUNTR ' HAS REACHED MGLL-1. If countr = mgll - 1 Then countr = 0 Else countr = countr + 1 End If End If ' RESET STEP ACCEPTANCE CODE AND DELSTR acpcod = 0 If linesr = 0 Then delstr = zero ' RESET QNFAIL TO FALSE TO AVOID PREMATURE STOPPING qnfail = 0 If linesr <> 0 Then ' THE MAIN LINE SEARCH IS CALLED IN SUBROUTINE LINE lineS abort, absnew, deuflh, geoms, newton, overch, overfl, qnfail, _ qrsing, restrt, sclfch, sclxch, acpcod, acptcr, contyp, _ isejac, itnum, jupdm, maxexp, maxlin, mgll, mnew, n, narmij, _ nfunc, nunit, output, qnupdm, stopcr, trmcod, alpha, confac, _ epsmch, fcnmax, fcnnew, fcnold, lam0, maxstp, newlen, sbrnrm, _ sigma, a, boundl, boundu, delf, ftrack, fvec, h, hhpi, jac, _ rdiag, wv1, s, sbar, scalef, scalex, sn, strack, xc, xplus If abort <> 0 Then 'Line not ok If output > 0 Then Line0 Line1 End If Exit Sub End If ' NOTE: 201 PRECEDES PRINTING OF ITERATION RESULTS If (newton) Then GoTo 320 Else ' TRUST REGION METHOD ' ESTABLISH INITIAL TRUST REGION SIZE, DELTA, AND/OR ' FIND LENGTH OF SCALED DESCENT STEP, CAULEN. delcau cauchy, overch, overfl, isejac, maxexp, n, nunit, output, _ beta, caulen, delta, epsmch, maxstp, newlen, sqrtz, a, delf, scalex frstdg = 1 If output > 3 Then Line0 Line0 Msg " SUMMARY OF TRUST REGION METHOD USING (DOUBLE) DOGLEG STEP." End If ' MAIN INTERNAL LOOP FOR TRUST REGION METHOD. ' THE TRUST REGION SIZE IS STORED FOR COMPARISON ' LATER TO SET THE PARAMETER POWTAU USED IN POWELL'S ' TRUST REGION UPDATING SCHEME (QUASI-NEWTON, TRUPDM=1). delta0 = delta For notrst = 1 To maxtrs dogleg frstdg, newtkn, overch, overfl, maxexp, n, notrst, nunit, _ output, beta, caulen, delta, etafac, newlen, sqrtz, delf, _ s, scalex, sn, ssdhat, vhat ' NOTE: WV1 AND WV4 HOLD THE COMPONENT AND RESIDUAL VECTOR ' RESPECTIVELY FOR A TRIAL POINT WHICH HAS BEEN FOUND ' TO BE ACCEPTABLE WHILE THE TRUST REGION IS EXPANDED ' AND A NEW TRIAL POINT TESTED. ' WV2 AND WV3 ARE WORK VECTORS. ' H IS CALLED ASTORE INSIDE TRSTUP. trstup geoms, newtkn, overch, overfl, qrsing, sclfch, sclxch, _ acpcod, acpstr, acptcr, contyp, isejac, jupdm, maxexp, _ mgll, mnew, n, narmij, nfunc, notrst, _ qnupdm, retcod, trupdm, alpha, confac, delfac, delstr, _ delta, epsmch, fcnmax, fcnnew, fcnold, fcnpre, maxstp, _ newlen, newmax, powtau, rellen, stptol, a, h, boundl, _ boundu, delf, wv1, ftrack, fvec, fvecc, hhpi, jac, rdiag, _ wv2, s, sbar, scalef, scalex, strack, xc, wv4, xplus If (output > 4 Or retcod = 7) And output > 2 Then rcdprt nunit, retcod, delta, rellen, stptol End If ' IF NO PROGRESS WAS BEING MADE (RETCOD=7) IN A QUASI-NEWTON ' STEP RETRY WITH AN EXPLICIT JACOBIAN EVALUATION. If retcod = 7 And restrt = 0 Then qnfail = 1 ' RETURN CODE LESS THAN 8 EXITS FROM TRUST REGION LOOP If retcod < 8 Then GoTo 310 Next notrst ' IF NO SUCCESSFUL STEP FOUND IN A QUASI-NEWTON STEP, ' RETRY WITN AN EXPLICIT JACOBIAN EVALUATION. If restrt = 0 Then qnfail = 1 End If 310 If linesr = 0 Then If delta < delta0 Then powtau = one delta = XMAX(delta, 0.0000000001) End If ' IF RETCOD=7 AND STOPCR=2, RESET STOPPING CRITERION TO ' AVOID HANGING IN TRUST REGION METHOD. (RETCOD=7 MEANS ' THE RELATIVE STEP LENGTH WAS LESS THAN STPTOL). If linesr = 0 And retcod = 7 And stopcr = 2 And qnfail = 0 Then stopcr = 12 ' RETAIN NUMBER OF STEPS ACCEPTED BY EACH CRITERION FOR PERFORMANCE EVALUATION. If newton = 0 Then If acpcod = 1 Then nac1 = nac1 + 1 ElseIf acpcod = 2 Then nac2 = nac2 + 1 ElseIf acpcod = 12 Then nac12 = nac12 + 1 End If End If ' *** PRINT RESULTS OF ITERATION *** 320 If output > 2 Then nersl newton, restrt, sclfch, sclxch, acpcod, jupdm, n, nunit, output, fcnnew, fvec, xplus End If ' CHECK FOR CONVERGENCE. STATEMENT 202 IS USED IF THE ' STEP SIZE OF THE NEWTON STEP IS FOUND TO BE WITHIN ' THE SPECIFIED TOLERANCE AND STOPCR IS 1 OR 12 330 If qnfail = 0 Then ' IF QNFAIL IS TRUE THE QUASI-NEWTON SEARCH FAILED TO ' FIND A SATISFACTORY STEP - SINCE THE JACOBIAN IS TO ' BE RE-EVALUATED AVOID PREMATURE STOPPAGES IN NESTOP nestop absnew, linesr, newton, sclfch, sclxch, acptcr, itnum, n, _ nac1, nac2, nac12, nfunc, njetot, nunit, output, stopcr, _ trmcod, fcnnew, ftol, nsttol, stpmax, stptol, fvec, scalef, _ scalex, xc, xplus End If ' IF THE TERMINATION CODE, TRMCOD, IS GREATER THAN 0 THEN ' CONVERGENCE HAS BEEN REACHED. If trmcod > 0 Then Exit Sub 'normal exit ' QUASI-NEWTON UPDATING - JUPDM > 0 If jupdm > 0 Then ' QNFAIL MEANS A FAILURE IN THE QUASI-NEWTON SEARCH. ' RE-EVALUATE JACOBIAN AND TRY A DAMPED NEWTON STEP. ' MAXLIN IS CHANGED FROM MAXQNS TO MAXNS OR MAXTRS IS ' CHANGED, SIMILARLY, AT THE START OF THE LOOP. If qnfail <> 0 Then If output > 4 Then Line0 Msg " FAILURE IN QUASI-NEWTON SEARCH: QNFAIL IS TRUE." End If ElseIf newton = 0 Then ' QNFAIL IS FALSE - THE STEP HAS BEEN ACCEPTED If output > 4 Then Line0 ss = " FCNNEW= " + Str$(fcnnew) + " FCNOLD= " + Str$(fcnold) + " RATIO= " + Str$(fcnnew / fcnold): Msg ss End If If fcnnew / fcnold > ratiof Then ' STEP ACCEPTED BUT NOT A SIGNIFICANT IMPROVEMENT nfail = nfail + 1 If output > 4 Then Line0 ss = " RATIO > RATIOF SO NFAIL INCREASED TO: " + Str$(nfail): Msg ss Line0 End If Else ' STEP ACCEPTED WITH A SIGNIFICANT IMPROVEMENT If fcnnew / fcnold > pt01 Then nosubt = 0 Else nosubt = 1 End If ' ITEMPT IS USED LOCALLY FOR OUTPUT CONTROL itemp = nfail nfail = IMAX(nfail - nosubt, 0) If output > 4 Then Line0 If itemp = nfail Then ss = " NFAIL STAYS AT: " + Str$(nfail): Msg ss Else ss = " NFAIL CHANGED TO: " + Str$(nfail): Msg ss End If End If End If ' SAVE THE RESULTS FOR RESTART IF A FAILURE IN THE ' QUASI-NEWTON METHOD OCCURS - ESSENTIALLY THIS ' FINDS THE BEST POINT SO FAR If isejac = 1 Or nfail = 1 Or (nfail <= minqns And fcnnew / fcnmin < one) Then savest = 1 fcnmin = fcnnew End If If savest <> 0 Then savest = 0 fstore = fcnnew If output > 4 Then Line0 Msg " STEP IS SAVED" Line0 Msg " SAVED COMPONENT AND FUNCTION VALUES" Line0 End If For i = 1 To n xsave(i) = xplus(i) fsave(i) = fvec(i) If output > 4 Then ss = " XSAVE(" + Str$(i) + ") = " + Str$(xsave(i)) + " FSAVE(" + Str$(i) + ") = " + Str$(fsave(i)): Msg ss End If Next i If output > 4 Then Line0 itstr = itnum End If End If ' NOTE: IF QNFAIL IS TRUE THEN NFAIL CANNOT HAVE INCREASED IMPLYING ' THAT NFAIL CANNOT NOW BE GREATER THAN MINQNS If qnfail <> 0 Or nfail > minqns Then ' RESTART FROM BEST POINT FOUND SO FAR restrt = 1 If output > 4 Then Line0 Line0 Msg " RESTRT IS TRUE." End If For j = 0 To mgll - 1 ftrack(j) = zero Next j If acptcr = 12 Then For j = 0 To mgll - 1 strack(j) = zero Next j End If countr = 0 isejac = 0 mnew = 0 trmcod = 0 fcnold = fstore If output > 4 Then Line0 ss = " RETURN TO ITERATION: " + Str$(itstr) + " WHERE:": Msg ss Line0 End If For i = 1 To n xc(i) = xsave(i) fvecc(i) = fsave(i) If output > 4 Then ss = " XC(" + Str$(i) + ") = " + Str$(xc(i)) + " FVECC(" + Str$(i) + ") = " + Str$(fvecc(i)): Msg ss End If Next i If output > 4 Then Line0 Else If itnum >= niejev Then restrt = 0 End If End If ' UPDATE JACOBIAN IF DESIRED: ' QNUPDM = 0 ==> ACTUAL JACOBIAN BEING UPDATED ' QNUPDM = 1 ==> FACTORED JACOBIAN BEING UPDATED If restrt = 0 Then If qnupdm = 0 Then If jupdm = 1 Then ' USE BROYDEN UPDATE broyun overfl, maxexp, n, nunit, output, epsmch, fvec, fvecc, jac, scalex, xc, xplus ElseIf jupdm = 2 Then ' USE LEE AND LEE UPDATE llun overch, overfl, isejac, maxexp, n, nunit, output, epsmch, _ omega, fvec, fvecc, jac, plee, s, scalex, xc, xplus End If ' THE FACTORED FORM OF THE JACOBIAN IS UPDATED If jupdm = 1 Then broyfa overch, overfl, sclfch, sclxch, maxexp, n, nunit, _ output, epsmch, a, delf, fvec, fvecc, jac, rdiag, _ s, scalef, scalex, wv1, wv2, xc, xplus ElseIf jupdm = 2 Then llfa overch, overfl, sclfch, sclxch, isejac, maxexp, n, nunit, _ output, epsmch, omega, a, delf, fvec, fvecc, jac, plee, _ rdiag, s, scalef, scalex, xc, xplus End If End If End If ' UPDATE CURRENT VALUES - RESET TRMCOD TO ZERO. ' UPDATE M "VECTOR" (ACTUALLY ONLY THE LATEST VALUE IS NEEDED). mold = mnew If isejac < narmij Then mnew = 0 Else mnew = IMIN(mold + 1, mgll - 1) End If If jupdm = 0 Or (jupdm > 0 And restrt = 0) Then update mnew, mold, n, trmcod, fcnnew, fcnold, fvec, fvecc, xc, xplus End If Next itnum '*** itnum main loop *** If output > 0 Then Line1 Line0 ss = " NO SOLUTION FOUND AFTER " + Str$(itnum - 1) + " ITERATION(S).": Msg ss Line0 Msg " FINAL ESTIMATES FINAL FUNCTION VALUES" Line0 For i = 1 To n ss = " X(" + Str$(i) + ") = " + Str$(Int(xplus(i) * 100000) / 100000) + " F(" + Str$(i) + ") = " + _ Str$(fvec(i)) Msg ss Next i Line0 ss = " FINAL OBJECTIVE FUNCTION VALUE = " + Str$(fcnnew): Msg ss Line0 Line1 End If End Sub