Attribute VB_Name = "Module1" '****************************************************************************** '* 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 1/3) * '* (www.jpmoreau.fr) * '****************************************************************************** DefInt I-N DefDbl A-H, O-Z Public Const ISIZE = 25 Public Const one = 1# Public Const ten = 10# Public Const two = 2# Public Const zero = 0# Public bypass As Integer Public matsup As Integer Public wrnsup As Integer Public itnum As Integer Public nfunc As Integer Public smallb As Double Public bigb As Double Public smalls As Double Public bigs As Double Public bigr As Double Sub abmul(nradec, nraact, ncndec, ncbact, ncdec, ncact, amat(), bmat(), cmat(), arow()) '------------------------------------------------------------------------ ' FEB. 8, 1991 ' ' MATRIX MULTIPLICATION AB=C ' ' VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 ' EACH ROW OF MATRIX A IS SAVED AS A COLUMN, AROW, BEFORE USE. ' ' NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX ' ' NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT FOR 2ND INDEX ' NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT ' ' I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED ' NCBDEC IS THE NUMBER OF COLUMNS OF B DECLARED ' NCDEC IS THE COMMON DECLARED DIMENSION ' ' MODIFICATION OF MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ' QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA '------------------------------------------------------------------------ ' FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32 = ncact / 32 ncc32r = ncact - 32 * ncc32 ncc16 = ncc32r / 16 ncc16r = ncc32r - 16 * ncc16 ncc8 = ncc16r / 8 ncc8r = ncc16r - 8 * ncc8 ncc4 = ncc8r / 4 ncc4r = ncc8r - 4 * ncc4 ' REASSIGN ROWS TO VECTOR AROW For i = 1 To nraact k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 arow(k - 31) = amat((i - 1) * ISIZE + k - 31) arow(k - 30) = amat((i - 1) * ISIZE + k - 30) arow(k - 29) = amat((i - 1) * ISIZE + k - 29) arow(k - 28) = amat((i - 1) * ISIZE + k - 28) arow(k - 27) = amat((i - 1) * ISIZE + k - 27) arow(k - 26) = amat((i - 1) * ISIZE + k - 26) arow(k - 25) = amat((i - 1) * ISIZE + k - 25) arow(k - 24) = amat((i - 1) * ISIZE + k - 24) arow(k - 23) = amat((i - 1) * ISIZE + k - 23) arow(k - 22) = amat((i - 1) * ISIZE + k - 22) arow(k - 21) = amat((i - 1) * ISIZE + k - 21) arow(k - 20) = amat((i - 1) * ISIZE + k - 20) arow(k - 19) = amat((i - 1) * ISIZE + k - 19) arow(k - 18) = amat((i - 1) * ISIZE + k - 18) arow(k - 17) = amat((i - 1) * ISIZE + k - 17) arow(k - 16) = amat((i - 1) * ISIZE + k - 16) arow(k - 15) = amat((i - 1) * ISIZE + k - 15) arow(k - 14) = amat((i - 1) * ISIZE + k - 14) arow(k - 13) = amat((i - 1) * ISIZE + k - 13) arow(k - 12) = amat((i - 1) * ISIZE + k - 12) arow(k - 11) = amat((i - 1) * ISIZE + k - 11) arow(k - 10) = amat((i - 1) * ISIZE + k - 10) arow(k - 9) = amat((i - 1) * ISIZE + k - 9) arow(k - 8) = amat((i - 1) * ISIZE + k - 8) arow(k - 7) = amat((i - 1) * ISIZE + k - 7) arow(k - 6) = amat((i - 1) * ISIZE + k - 6) arow(k - 5) = amat((i - 1) * ISIZE + k - 5) arow(k - 4) = amat((i - 1) * ISIZE + k - 4) arow(k - 3) = amat((i - 1) * ISIZE + k - 3) arow(k - 2) = amat((i - 1) * ISIZE + k - 2) arow(k - 1) = amat((i - 1) * ISIZE + k - 1) arow(k) = amat((i - 1) * ISIZE + k) Next kk End If If ncc16 > 0 Then For kk = 1 To ncc16 k = k + 16 arow(k - 15) = amat((i - 1) * ISIZE + k - 15) arow(k - 14) = amat((i - 1) * ISIZE + k - 14) arow(k - 13) = amat((i - 1) * ISIZE + k - 13) arow(k - 12) = amat((i - 1) * ISIZE + k - 12) arow(k - 11) = amat((i - 1) * ISIZE + k - 11) arow(k - 10) = amat((i - 1) * ISIZE + k - 10) arow(k - 9) = amat((i - 1) * ISIZE + k - 9) arow(k - 8) = amat((i - 1) * ISIZE + k - 8) arow(k - 7) = amat((i - 1) * ISIZE + k - 7) arow(k - 6) = amat((i - 1) * ISIZE + k - 6) arow(k - 5) = amat((i - 1) * ISIZE + k - 5) arow(k - 4) = amat((i - 1) * ISIZE + k - 4) arow(k - 3) = amat((i - 1) * ISIZE + k - 3) arow(k - 2) = amat((i - 1) * ISIZE + k - 2) arow(k - 1) = amat((i - 1) * ISIZE + k - 1) arow(k) = amat((i - 1) * ISIZE + k) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 arow(k - 7) = amat((i - 1) * ISIZE + k - 7) arow(k - 6) = amat((i - 1) * ISIZE + k - 6) arow(k - 5) = amat((i - 1) * ISIZE + k - 5) arow(k - 4) = amat((i - 1) * ISIZE + k - 4) arow(k - 3) = amat((i - 1) * ISIZE + k - 3) arow(k - 2) = amat((i - 1) * ISIZE + k - 2) arow(k - 1) = amat((i - 1) * ISIZE + k - 1) arow(k) = amat((i - 1) * ISIZE + k) Next kk End If If ncc4 > 0 Then For kk = 1 To ncc4 k = k + 4 arow(k - 3) = amat((i - 1) * ISIZE + k - 3) arow(k - 2) = amat((i - 1) * ISIZE + k - 2) arow(k - 1) = amat((i - 1) * ISIZE + k - 1) arow(k) = amat((i - 1) * ISIZE + k) Next kk End If If ncc4r > 0 Then For kk = 1 To ncc4r k = k + 1 arow(k) = amat((i - 1) * ISIZE + k) Next kk End If ' FIND ENTRY FOR MATRIX C USING COLUMN VECTOR AROW For j = 1 To ncbact Sum = zero k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 Sum = Sum + arow(k - 31) * bmat((k - 32) * ISIZE + j) + arow(k - 30) * bmat((k - 31) * ISIZE + j) + _ arow(k - 29) * bmat((k - 30) * ISIZE + j) + arow(k - 28) * bmat((k - 29) * ISIZE + j) + _ arow(k - 27) * bmat((k - 28) * ISIZE + j) + arow(k - 26) * bmat((k - 27) * ISIZE + j) + _ arow(k - 25) * bmat((k - 26) * ISIZE + j) + arow(k - 24) * bmat((k - 25) * ISIZE + j) Sum = Sum + arow(k - 23) * bmat((k - 24) * ISIZE + j) + arow(k - 22) * bmat((k - 23) * ISIZE + j) + _ arow(k - 21) * bmat((k - 22) * ISIZE + j) + arow(k - 20) * bmat((k - 21) * ISIZE + j) + _ arow(k - 19) * bmat((k - 20) * ISIZE + j) + arow(k - 18) * bmat((k - 19) * ISIZE + j) + _ arow(k - 17) * bmat((k - 18) * ISIZE + j) + arow(k - 16) * bmat((k - 17) * ISIZE + j) Sum = Sum + arow(k - 15) * bmat((k - 16) * ISIZE + j) + arow(k - 14) * bmat((k - 15) * ISIZE + j) + _ arow(k - 13) * bmat((k - 14) * ISIZE + j) + arow(k - 12) * bmat((k - 13) * ISIZE + j) + _ arow(k - 11) * bmat((k - 12) * ISIZE + j) + arow(k - 10) * bmat((k - 11) * ISIZE + j) + _ arow(k - 9) * bmat((k - 10) * ISIZE + j) + arow(k - 8) * bmat((k - 9) * ISIZE + j) Sum = Sum + arow(k - 7) * bmat((k - 8) * ISIZE + j) + arow(k - 6) * bmat((k - 7) * ISIZE + j) + _ arow(k - 5) * bmat((k - 6) * ISIZE + j) + arow(k - 4) * bmat((k - 5) * ISIZE + j) + _ arow(k - 3) * bmat((k - 4) * ISIZE + j) + arow(k - 2) * bmat((k - 3) * ISIZE + j) + _ arow(k - 1) * bmat((k - 2) * ISIZE + j) + arow(k) * bmat((k - 1) * ISIZE + j) Next kk End If If ncc16 > 0 Then For kk = 1 To ncc16 k = k + 16 Sum = Sum + arow(k - 15) * bmat((k - 16) * ISIZE + j) + arow(k - 14) * bmat((k - 15) * ISIZE + j) + _ arow(k - 13) * bmat((k - 14) * ISIZE + j) + arow(k - 12) * bmat((k - 13) * ISIZE + j) + _ arow(k - 11) * bmat((k - 12) * ISIZE + j) + arow(k - 10) * bmat((k - 11) * ISIZE + j) + _ arow(k - 9) * bmat((k - 10) * ISIZE + j) + arow(k - 8) * bmat((k - 9) * ISIZE + j) Sum = Sum + arow(k - 7) * bmat((k - 8) * ISIZE + j) + arow(k - 6) * bmat((k - 7) * ISIZE + j) + _ arow(k - 5) * bmat((k - 6) * ISIZE + j) + arow(k - 4) * bmat((k - 5) * ISIZE + j) + _ arow(k - 3) * bmat((k - 4) * ISIZE + j) + arow(k - 2) * bmat((k - 3) * ISIZE + j) + _ arow(k - 1) * bmat((k - 2) * ISIZE + j) + arow(k) * bmat((k - 1) * ISIZE + j) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 Sum = Sum + arow(k - 7) * bmat((k - 8) * ISIZE + j) + arow(k - 6) * bmat((k - 7) * ISIZE + j) + _ arow(k - 5) * bmat((k - 6) * ISIZE + j) + arow(k - 4) * bmat((k - 5) * ISIZE + j) + _ arow(k - 3) * bmat((k - 4) * ISIZE + j) + arow(k - 2) * bmat((k - 3) * ISIZE + j) + _ arow(k - 1) * bmat((k - 2) * ISIZE + j) + arow(k) * bmat((k - 1) * ISIZE + j) Next kk End If If ncc4 > 0 Then For kk = 1 To ncc4 k = k + 4 Sum = Sum + arow(k - 3) * bmat((k - 4) * ISIZE + j) + arow(k - 2) * bmat((k - 3) * ISIZE + j) + _ arow(k - 1) * bmat((k - 2) * ISIZE + j) + arow(k) * bmat((k - 1) * ISIZE + j) Next kk End If If ncc4r > 0 Then For kk = 1 To ncc4r k = k + 1 Sum = Sum + arow(k) * bmat((k - 1) * ISIZE + j) Next kk End If cmat((i - 1) * ISIZE + j) = Sum Next j Next i End Sub Sub ascalf(n, epsmch, fvecc(), xjac(), scalef()) '------------------------------------------------------------------------------ ' FEB. 13, 1991 ' ' THIS SUBROUTINE ESTABLISHES SCALING FACTORS FOR THE RESIDUAL VECTOR ' IF FUNCTION ADAPTIVE SCALING IS CHOSEN USING INTEGER VARIABLE ITSCLF. ' ' NOTE: IN QUASI-NEWTON METHODS THE SCALING FACTORS ARE ' UPDATED ONLY WHEN THE JACOBIAN IS EVALUATED EXPLICITLY. ' ' SCALING FACTORS ARE DETERMINED FROM THE INFINITY NORMS OF THE ROWS ' OF THE JACOBIAN AND THE VALUES OF THE CURRENT FUNCTION VECTOR, FVECC. ' ' A MINIMUM TOLERANCE ON THE SCALING FACTOR IS THE SQUARE ' ROOT OF THE MACHINE PRECISION, SQRTEP. '------------------------------------------------------------------------------} sqrtep = Sqr(epsmch) ' I COUNTS THE ROWS For i = 1 To n amax = zero ' FIND MAXIMUM ENTRY IN ROW I For j = 1 To n amax = XMAX(amax, Abs(xjac((i - 1) * ISIZE + j))) Next j amax = XMAX(amax, fvecc(i)) ' SET SCALING FACTOR TO A DEFAULT OF ONE IF ITH ROW IS ZEROS If amax = zero Then amax = one amax = XMAX(amax, sqrtep) scalef(i) = one / amax Next i End Sub Sub ascalx(n, epsmch, xjac(), scalex()) '------------------------------------------------------------------------- ' FEB. 13, 1991 ' ' THIS SUBROUTINE ESTABLISHES SCALING FACTORS FOR THE COMPONENET VECTOR ' IF ADAPTIVE SCALING IS CHOSEN USING INTEGER ITSCLX. ' ' NOTE: IN QUASI-NEWTON METHODS THE SCALING FACTORS ARE ' UPDATED ONLY WHEN THE JACOBIAN IS EVALUATED EXPLICITLY. ' ' SCALING FACTORS ARE DETERMINED FROM THE INFINITY NORMS ' OF THE COLUMNS OF THE JACOBIAN. ' ' A MINIMUM TOLERANCE ON THE SCALING FACTOR IS THE SQUARE ' ROOT OF THE MACHINE PRECISION, SQRTEP. '------------------------------------------------------------------------- sqrtep = Sqr(epsmch) ' J COUNTS COLUMNS For j = 1 To n amax = zero ' FIND MAXIMUM ENTRY IN JTH COLUMN For i = 1 To n amax = XMAX(amax, Abs(xjac((i - 1) * ISIZE + j))) Next i ' IF A COLUMN IS ALL ZEROS SET AMAX TO ONE If amax = zero Then amax = one scalex(j) = XMAX(amax, sqrtep) Next j End Sub Sub atamul(nradec, ncadec, nraact, ncaact, nrbdec, ncbdec, amat(), bmat()) '----------------------------------------------------------------------- ' FEB. 8, 1991 ' ' MATRIX MULTIPLICATION: A^A=B ' ' 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 NUMBER OF GROUPS OF SIZE 32, 16... ncc32 = nraact / 32 ncc32r = nraact - 32 * ncc32 ncc16 = ncc32r / 16 ncc16r = ncc32r - 16 * ncc16 ncc8 = ncc16r / 8 ncc8r = ncc16r - 8 * ncc8 ncc4 = ncc8r / 4 ncc4r = ncc8r - 4 * ncc4 ' FIND ENTRY IN MATRIX B For i = 1 To ncaact 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 Function Log10(x) Log10 = Log(x) / Log(ten) End Function Sub ataov(overfl As Integer, maxexp, n, nunit, output As Integer, _ a(), b(), scalef()) '---------------------------------------------------------------------- ' SEPT. 8, 1991 ' ' THIS SUBROUTINE FINDS THE PRODUCT OF THE TRANSPOSE OF THE MATRIX A ' AND MATRIX A. EACH ENTRY IS CHECKED BEFORE BEING ACCEPTED. ' IF IT WOULD CAUSE AN OVERFLOW 10^MAXEXP IS INSERTED IN ITS PLACE. '---------------------------------------------------------------------- 'Labels 40, 70 eps = one / (ten ^ maxexp) overfl = 0 For i = 1 To n For j = i + 1 To n Sum = zero For k = 1 To n If Log10(Abs(a((k - 1) * ISIZE + i)) + eps) + Log10(Abs(a((k - 1) * ISIZE + j)) _ + eps) + two * Log10(scalef(k)) > maxexp Then overfl = 1 b((i - 1) * ISIZE + j) = Sign(ten ^ maxexp, a((k - 1) * ISIZE + i)) * Sign(one, a((k - 1) * ISIZE + j)) If output > 2 And wrnsup = 0 Then Line0 Print #1, " * WARNING: COMPONENT MATRIX-MATRIX PRODUCT SET TO "; b((i - 1) * ISIZE + j); " " End If GoTo 40 End If Sum = Sum + a((k - 1) * ISIZE + i) * a((k - 1) * ISIZE + j) * scalef(k) * scalef(k) Next k b(i, j) = Sum b(j, i) = Sum 40 Next j Sum = zero For k = 1 To n If two * (Log10(Abs(a((k - 1) * ISIZE + i)) + eps) + Log10(scalef(k))) > maxexp Then overfl = 1 b((i - 1) * ISIZE + i) = ten ^ maxexp If output > 2 And wrnsup = 0 Then Line0 Print #1, " * WARNING: COMPONENT MATRIX-MATRIX PRODUCT SET TO "; b((i - 1) * ISIZE + i); " " End If GoTo 70 End If Sum = Sum + a((k - 1) * ISIZE + i) * a((k - 1) * ISIZE + i) * scalef(k) * scalef(k) Next k b((i - 1) * ISIZE + i) = Sum 70 Next i End Sub Function Sign(a, b) If b >= 0 Then Sign = a Else Sign = -a End If End Function Sub atbmul(ncadec, ncaact, ncbdec, ncbact, ncdec, ncact, amat(), bmat(), _ cmat()) '---------------------------------------------------------------------- ' FEB. 8, 1991 ' ' MATRIX MULTIPLICATION: A^B=C ' ' VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. ' ' NCADEC IS 2ND DIM. OF AMAT; NCAACT IS ACTUAL LIMIT FOR 2ND INDEX ' NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT FOR 2ND INDEX ' NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT ' ' I.E. NCADEC IS NUMBER OF COLUMNS OF A DECLARED ' NCBDEC IS NUMBER OF COLUMNS OF B DECLARED ' NCDEC IS THE NUMBER OF ROWS IN BOTH A AND B DECLARED ' ' MODIFIED VERSION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES ' MACKINNON, QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA '----------------------------------------------------------------------- ' FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32 = ncact / 32 ncc32r = ncact - 32 * ncc32 ncc16 = ncc32r / 16 ncc16r = ncc32r - 16 * ncc16 ncc8 = ncc16r / 8 ncc8r = ncc16r - 8 * ncc8 ncc4 = ncc8r / 4 ncc4r = ncc8r - 4 * ncc4 ' FIND ENTRY IN MATRIX C For i = 1 To ncaact For j = 1 To ncbact Sum = zero k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 Sum = Sum + amat((k - 32) * ISIZE + i) * bmat((k - 32) * ISIZE + j) + amat((k - 31) * ISIZE + i) * bmat((k - 31) * ISIZE + j) + _ amat((k - 30) * ISIZE + i) * bmat((k - 30) * ISIZE + j) + amat((k - 29) * ISIZE + i) * bmat((k - 29) * ISIZE + j) + _ amat((k - 28) * ISIZE + i) * bmat((k - 28) * ISIZE + j) + amat((k - 27) * ISIZE + i) * bmat((k - 27) * ISIZE + j) + _ amat((k - 26) * ISIZE + i) * bmat((k - 26) * ISIZE + j) + amat((k - 25) * ISIZE + i) * bmat((k - 25) * ISIZE + j) Sum = Sum + amat((k - 24) * ISIZE + i) * bmat((k - 24) * ISIZE + j) + amat((k - 23) * ISIZE + i) * bmat((k - 23) * ISIZE + j) + _ amat((k - 22) * ISIZE + i) * bmat((k - 22) * ISIZE + j) + amat((k - 21) * ISIZE + i) * bmat((k - 21) * ISIZE + j) + _ amat((k - 20) * ISIZE + i) * bmat((k - 20) * ISIZE + j) + amat((k - 19) * ISIZE + i) * bmat((k - 19) * ISIZE + j) + _ amat((k - 18) * ISIZE + i) * bmat((k - 18) * ISIZE + j) + amat((k - 17) * ISIZE + i) * bmat((k - 17) * ISIZE + j) Sum = Sum + amat((k - 16) * ISIZE + i) * bmat((k - 16) * ISIZE + j) + amat((k - 15) * ISIZE + i) * bmat((k - 15) * ISIZE + j) + _ amat((k - 14) * ISIZE + i) * bmat((k - 14) * ISIZE + j) + amat((k - 13) * ISIZE + i) * bmat((k - 13) * ISIZE + j) + _ amat((k - 12) * ISIZE + i) * bmat((k - 12) * ISIZE + j) + amat((k - 11) * ISIZE + i) * bmat((k - 11) * ISIZE + j) + _ amat((k - 10) * ISIZE + i) * bmat((k - 10) * ISIZE + j) + amat((k - 9) * ISIZE + i) * bmat((k - 9) * ISIZE + j) Sum = Sum + amat((k - 8) * ISIZE + i) * bmat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * bmat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * bmat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * bmat((k - 5) * ISIZE + j) + _ amat((k - 4) * ISIZE + i) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * bmat((k - 3) * ISIZE + j) + _ amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 16) * ISIZE + j) + amat((k - 15) * ISIZE + i) _ * bmat((k - 15) * ISIZE + j) + amat((k - 14) * ISIZE + i) * bmat((k - 14) * ISIZE + j) + _ amat((k - 13) * ISIZE + i) * bmat((k - 13) * ISIZE + j) + amat((k - 12) * ISIZE + i) * bmat((k - 12) * ISIZE + j) + _ amat((k - 11) * ISIZE + i) * bmat((k - 11) * ISIZE + j) + amat((k - 10) * ISIZE + i) * bmat((k - 10) * ISIZE + j) + _ amat((k - 9) * ISIZE + i) * bmat((k - 9) * ISIZE + j) Sum = Sum + amat((k - 8) * ISIZE + i) * bmat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) _ * bmat((k - 7) * ISIZE + j) + amat((k - 6) * ISIZE + i) * bmat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) _ * bmat((k - 5) * ISIZE + j) + amat((k - 4) * ISIZE + i) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) _ * bmat((k - 3) * ISIZE + j) + amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + _ amat((k - 1) * ISIZE + i) * bmat((k - 1) * ISIZE + j) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 Sum = Sum + amat((k - 7) * ISIZE + i) * bmat((k - 7) * ISIZE + j) + amat((k - 6) * ISIZE + i) _ * bmat((k - 6) * ISIZE + j) + amat((k - 6) * ISIZE + i) * bmat _ ((k - 5) * ISIZE + j) + amat((k - 5) * ISIZE + i) * bmat((k - 5) * ISIZE + j) + _ amat((k - 4) * ISIZE + i) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * bmat((k - 3) * ISIZE + j) + _ amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) _ * bmat((k - 3) * ISIZE + j) + amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + _ amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 1) * ISIZE + j) Next kk End If cmat((i - 1) * ISIZE + j) = Sum Next j Next i End Sub Sub atvov(overfl As Integer, maxexp, n, nunit, output As Integer, amat(), bvec(), _ cvec()) '---------------------------------------------------------------------- ' FEB. 8 ,1991 ' ' THIS SUBROUTINE FINDS THE PRODUCT OF THE TRANSPOSE OF THE MATRIX A ' AND THE VECTOR B WHERE EACH ENTRY IS CHECKED TO PREVENT OVERFLOWS. '---------------------------------------------------------------------- eps = one / (ten ^ maxexp) overfl = 0 For i = 1 To n Sum = zero For j = 1 To n If Log10(Abs(amat((j - 1) * ISIZE + i)) + eps) + Log10(Abs(bvec(j)) + eps) > maxexp Then overfl = 1 cvec(i) = Sign(ten ^ maxexp, amat((j - 1) * ISIZE + i)) * Sign(one, bvec(j)) If output > 2 And wrnsup = 0 Then Line0 Print #1, " * WARNING: COMPONENT MATRIX-VECTOR PRODUCT SET TO "; cvec(i); " *" End If Exit Sub End If Next j Sum = Sum + amat((j - 1) * ISIZE + i) * bvec(j) Next i cvec(i) = Sum End Sub Sub avmul(nradec, nraact, ncdec, ncact, amat(), bvec(), cvec()) '---------------------------------------------------------------------------------- ' FEB. 8, 1991 ' ' MATRIX-VECTOR MULTIPLICATION AB=C ' ' VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 ' EACH ROW OF MATRIX A IS SAVED AS A COLUMN BEFORE USE. ' ' NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX ' NCDEC IS COMMON DIMENSION OF AMAT & BVEC; NCACT IS ACTUAL LIMIT ' ' I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED ' NCDEC IS THE COMMON DECLARED DIMENSION (COLUMNS OF A AND ROWS OF B) ' ' MODIFICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ' QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA '---------------------------------------------------------------------------------- ' FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32 = ncact / 32 ncc32r = ncact - 32 * ncc32 ncc16 = ncc32r / 16 ncc16r = ncc32r - 16 * ncc16 ncc8 = ncc16r / 8 ncc8r = ncc16r - 8 * ncc8 ncc4 = ncc8r / 4 ncc4r = ncc8r - 4 * ncc4 For i = 1 To nraact ' FIND ENTRY FOR VECTOR C Sum = zero k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 Sum = Sum + amat((i - 1) * ISIZE + k - 31) * bvec(k - 31) + amat((i - 1) * ISIZE + k - 30) _ * bvec(k - 30) + amat((i - 1) * ISIZE + k - 29) * bvec(k - 29) + amat((i - 1) * ISIZE + k - 28) * bvec(k - 28) + _ amat((i - 1) * ISIZE + k - 27) * bvec(k - 27) + amat((i - 1) * ISIZE + k - 26) * bvec(k - 26) + _ amat((i - 1) * ISIZE + k - 25) * bvec(k - 25) + amat((i - 1) * ISIZE + k - 24) * bvec(k - 24) Sum = Sum + amat((i - 1) * ISIZE + k - 23) * bvec(k - 23) + amat((i - 1) * ISIZE + k - 22) _ * bvec(k - 22) + amat((i - 1) * ISIZE + k - 21) * bvec(k - 21) + amat((i - 1) * ISIZE + k - 20) * bvec(k - 20) + _ amat((i - 1) * ISIZE + k - 19) * bvec(k - 19) + amat((i - 1) * ISIZE + k - 18) * bvec(k - 18) + _ amat((i - 1) * ISIZE + k - 17) * bvec(k - 17) + amat((i - 1) * ISIZE + k - 16) * bvec(k - 16) Sum = Sum + amat((i - 1) * ISIZE + k - 15) * bvec(k - 15) + amat((i - 1) * ISIZE + k - 14) _ * bvec(k - 14) + amat((i - 1) * ISIZE + k - 13) * bvec(k - 13) + amat((i - 1) * ISIZE + k - 12) * bvec(k - 12) + _ amat((i - 1) * ISIZE + k - 11) * bvec(k - 11) + amat((i - 1) * ISIZE + k - 10) * bvec(k - 10) + _ amat((i - 1) * ISIZE + k - 9) * bvec(k - 9) + amat((i - 1) * ISIZE + k - 8) * bvec(k - 8) Sum = Sum + amat((i - 1) * ISIZE + k - 7) * bvec(k - 7) + amat((i - 1) * ISIZE + k - 6) * bvec(k - 6) + _ amat((i - 1) * ISIZE + k - 5) * bvec(k - 5) + amat((i - 1) * ISIZE + k - 4) * bvec(k - 4) + _ amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + _ amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc16 > 0 Then For kk = 1 To ncc16 k = k + 16 Sum = Sum + amat((i - 1) * ISIZE + k - 15) * bvec(k - 15) + amat((i - 1) * ISIZE + k - 14) * bvec(k - 14) + _ amat((i - 1) * ISIZE + k - 13) * bvec(k - 13) + amat((i - 1) * ISIZE + k - 12) * bvec(k - 12) + _ amat((i - 1) * ISIZE + k - 11) * bvec(k - 11) + amat((i - 1) * ISIZE + k - 10) * bvec(k - 10) + _ amat((i - 1) * ISIZE + k - 9) * bvec(k - 9) + amat((i - 1) * ISIZE + k - 8) * bvec(k - 8) Sum = Sum + amat((i - 1) * ISIZE + k - 7) * bvec(k - 7) + amat((i - 1) * ISIZE + k - 6) * bvec(k - 6) + _ amat((i - 1) * ISIZE + k - 5) * bvec(k - 5) + amat((i - 1) * ISIZE + k - 4) * bvec(k - 4) + _ amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + _ amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 Sum = Sum + amat((i - 1) * ISIZE + k - 7) * bvec(k - 7) + amat((i - 1) * ISIZE + k - 6) * bvec(k - 6) + _ amat((i - 1) * ISIZE + k - 5) * bvec(k - 5) + amat((i - 1) * ISIZE + k - 4) * bvec(k - 4) + _ amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + _ amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc4 > 0 Then For kk = 1 To ncc4 k = k + 4 Sum = Sum + amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + _ amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc4r > 0 Then For kk = 1 To ncc4r k = k + 1 Sum = Sum + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If cvec(i) = Sum Next i End Sub Sub bakdif(overfl As Integer, j, n, deltaj, tempj, fvec(), fvecj1(), _ jacfdm() As Double, xc()) ' FEB. 6, 1991 deltaj = tempj - xc(j) fcn1 overfl, n, fvecj1, xc If overfl = 0 Then For i = 1 To n jacfdm((i - 1) * ISIZE + j) = (fvec(i) - fvecj1(i)) / deltaj Next i End If End Sub Sub bnddif(overfl As Integer, j, n, epsmch, boundl(), boundu(), _ fvecc(), fvecj1(), jacfdm() As Double, xc()) '------------------------------------------------------------------------------- ' FEB. 15, 1991 ' ' FINITE DIFFERENCE CALCULATION WHEN THE BOUNDS FOR COMPONENT J ARE SO CLOSE ' THAT NEITHER A FORWARD NOR BACKWARD DIFFERENCE CAN BE PERFORMED. '------------------------------------------------------------------------------- Dim wv3(ISIZE) eps3q = epsmch ^ 0.75 ' STORE CURRENT For i = 1 To n wv3(i) = fvecc(i) Next i xc(j) = boundu(j) fcn1 overfl, n, fvecj1, xc If overfl = 0 Then xc(j) = boundl(j) fcn1 overfl, n, fvecc, xc If overfl = 0 Then For i = 1 To n ' ENSURE THAT THE JACOBIAN CALCULATION ISN'T JUST NOISE If fvecj1(i) - fvecc(i) > eps3q Then jacfdm((i - 1) * ISIZE + j) = (fvecj1(i) - fvecc(i)) / (boundu(j) - boundl(j)) Else jacfdm((i - 1) * ISIZE + j) = zero End If Next i End If End If For i = 1 To n fvecc(i) = wv3(i) Next i End Sub Sub broyfa(overch As Integer, overfl As Integer, sclfch As Integer, sclxch As Integer, _ maxexp, n, nunit, output As Integer, epsmch, a(), delf(), fvec(), _ fvecc(), jac() As Double, rdiag(), s(), scalef(), scalex(), t(), w(), xc(), xplus()) '-------------------------------------------------------------------------- ' FEB. 23, 1992 ' ' THE BROYDEN QUASI-NEWTON METHOD IS APPLIED TO THE FACTORED ' FORM OF THE JACOBIAN. ' ' NOTE: T AND W ARE TEMPORARY WORKING VECTORS ONLY. ' ' THE UPDATE OCCURS ONLY IF A SIGNIFICANT CHANGE IN THE JACOBIAN ' WOULD RESULT, I.E., NOT ALL THE VALUES IN VECTOR W ARE LESS THAN ' THE THRESHOLD IN MAGNITUDE. IF THE VECTOR W IS ESSENTIALLY ZERO THEN ' THE LOGICAL VARIABLE SKIPUP REMAINS SET AT TRUE. '-------------------------------------------------------------------------- Dim skipup As Integer 'here boolean (0 or 1) Dim tmp(ISIZE) Dim tmp1(ISIZE * ISIZE), tmp2(ISIZE * ISIZE) overfl = 0 eps = one / (ten ^ maxexp) sqrtep = Sqr(epsmch) For i = 1 To n a((i - 1) * ISIZE + i) = rdiag(i) s(i) = xplus(i) - xc(i) Next i ' R IS NOW IN THE UPPER TRIANGLE OF A skipup = 1 ' THE BROYDEN UPDATE IS CONDENSED INTO THE FORM ' A(NEW) = A(OLD) + T S^ ' THE PRODUCT A*S IS FORMED IN TWO STAGES AS R IS IN THE UPPER ' TRIANGLE OF MATRIX A AND Q^ IS IN JAC. ' FIRST MULTIPLY R*S (A IS CONSIDERED UPPER TRIANGULAR) uvmul n, n, n, n, a, s, t ' NOTE: THIS T IS TEMPORARY - IT IS THE T FROM BELOW WHICH ' IS SENT TO SUBROUTINE QRUPDA. For i = 1 To n For k = 1 To n tmp(k) = jac((k - 1) * ISIZE + i) Next k innerp overch, overfl, maxexp, n, n, n, nunit, output, Sum, tmp, t w(i) = scalef(i) * (fvec(i) - fvecc(i)) - Sum ' TEST TO ENSURE VECTOR W IS NONZERO. ANY VALUE GREATER ' THAN THE THRESHOLD WILL SET SKIPUP TO 0=FALSE. If Abs(w(i)) > sqrtep * scalef(i) * (Abs(fvec(i)) + Abs(fvecc(i))) Then skipup = 0 'False Else w(i) = zero End If Next i ' IF W(I)=0 FOR ALL I THEN THE UPDATE IS SKIPPED. If skipup = 0 Then ' T=Q^W; Q^ IS IN JAC avmul n, n, n, n, jac, w, t If sclxch = 1 Then For k = 1 To n w(k) = s(k) * scalex(k) Next k Else For k = 1 To n tmp1((k - 1) * ISIZE + 1) = s(k): tmp2((k - 1) * ISIZE + 1) = w(k) Next k matcop n, n, 1, 1, n, 1, tmp1, tmp2 For k = 1 To n s(k) = tmp1((k - 1) * ISIZE + 1): w(k) = tmp2((k - 1) * ISIZE + 1) Next k End If twonrm overfl, maxexp, n, epsmch, denom, w ' IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN If overfl <> 0 Or Log10(denom + eps) > 0.5 * maxexp Then If output > 3 Then Line0 Msg " WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF" Msg " BROYDEN UPDATE." End If Exit Sub Else denom = denom * denom End If ' IF DENOM IS ZERO AVOID /IDE BY ZERO AND CONTINUE WITH SAME JACOBIAN If denom = zero Then Exit Sub ' THE SCALED VERSION OF S REPLACES THE ORIGINAL BEFORE ' BEING SENT TO QRUPDA. For k = 1 To n s(k) = s(k) * scalex(k) * scalex(k) / denom Next k ' UPDATE THE QR DECOMPOSITION USING A SERIES OF GIVENS ROTATIONS. qrupda overfl, maxexp, n, epsmch, a, jac, t, s ' RESET RDIAG AS DIAGONAL OF CURRENT R WHICH IS IN ' THE UPPER TRIANGE OF A. For i = 1 To n rdiag(i) = a((i - 1) * ISIZE + i) Next i End If ' UPDATE THE GRADIENT VECTOR, DELF. THE NEW Q^ IS IN JAC. ' 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 tmp1((k - 1) * ISIZE + 1) = fvec(k): tmp2((k - 1) * ISIZE + 1) = w(k) Next k matcop n, n, 1, 1, n, 1, tmp1, tmp2 For k = 1 To n fvec(k) = tmp1((k - 1) * ISIZE + 1): w(k) = tmp2((k - 1) * ISIZE + 1) Next k End If avmul n, n, n, n, jac, w, t For k = 1 To n tmp1((k - 1) * ISIZE + 1) = t(k): tmp2((k - 1) * ISIZE + 1) = delf(k) Next k utbmul n, n, 1, 1, n, n, a, tmp1, tmp2 For k = 1 To n t(k) = tmp1((k - 1) * ISIZE + 1): delf(k) = tmp2((k - 1) * ISIZE + 1) Next k End Sub Sub broyun(overfl As Integer, maxexp, n, nunit, output As Integer, _ epsmch, fvec(), fvecc(), jac() As Double, scalex(), xc(), xplus()) '--------------------------------------------------------------------- ' FEB. 23, 1992 ' ' UPDATE THE JACOBIAN USING BROYDEN'S METHOD. '--------------------------------------------------------------------- Dim wv1(ISIZE), tmp1(ISIZE), tmp2(ISIZE) eps = one / (ten ^ maxexp) sqrtep = Sqr(epsmch) For k = 1 To n wv1(k) = (xplus(k) - xc(k)) * scalex(k) Next k twonrm overfl, maxexp, n, epsmch, denom, wv1 ' IF OVERFLOW WOULD OCCUR MAKE NO CHANGE TO JACOBIAN If overfl <> 0 Or Log10(denom + eps) > 0.5 * maxexp Then If output > 3 Then Line0 Msg " WARNING: JACOBIAN NOT UPDATED TO AVOID OVERFLOW IN DENOMINATOR OF" Msg " BROYDEN UPDATE." End If Exit Sub Else denom = denom * denom End If ' IF DENOM IS ZERO, AVOID OVERFLOW, CONTINUE WITH SAME JACOBIAN If denom = zero Then Exit Sub ' UPDATE JACOBIAN BY ROWS For i = 1 To n For k = 1 To n tmp1(k) = jac((i - 1) * ISIZE + k): tmp2(k) = xplus(k) - xc(k) Next k Sum = DotProduct(n, tmp1, tmp2) tempi = fvec(i) - fvecc(i) - Sum ' CHECK TO ENSURE THAT SOME MEANINGFUL CHANGE IS BEING MADE ' TO THE APPROXIMATE JACOBIAN; IF NOT, SKIP UPDATING ROW I. 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 * (xplus(j) - xc(j)) * scalex(j) * scalex(j) Next j End If Next i End Sub Sub cholde(n, maxadd As Double, maxffl As Double, sqrtep, h(), l() As Double) '------------------------------------------------------------------------- ' FEB. 23, 1992 ' ' THIS SUBROUTINE FINDS THE CHOLESKY DECOMPOSITION OF THE MATRIX, H, ' AND RETURNS IT IN THE LOWER TRIANGLE OF MATRIX, L. '------------------------------------------------------------------------- Dim minl As Double Dim minl2 As Double Dim minljj As Double Dim tmp1(ISIZE), tmp2(ISIZE) ' minl=SQRT(sqrtep)*maxffl ' MAXFFL EQUALS 0 WHEN THE MATRIX IS KNOWN TO BE POSITIVE DEFINITE If maxffl = zero Then ' FIND SQUARE ROOT OF LARGEST MAGNITUDE DIAGONAL ELEMENT ' AND SET MINL2. For i = 1 To n maxffl = XMAX(maxffl, Abs(h((i - 1) * ISIZE + i))) Next i maxffl = Sqr(maxffl) minl2 = sqrtep * maxffl End If ' MAXADD CONTAINS THE MAXIMUM AMOUNT WHICH IS IMPLICITLY ADDED ' TO ANY DIAGONAL ELEMENT OF MATRIX H. maxadd = zero For j = 1 To n For k = 1 To j - 1 tmp1(k) = l((j - 1) * ISIZE + k) Next k Sum = DotProduct(j - 1, tmp1, tmp1) l((j - 1) * ISIZE + j) = h((j - 1) * ISIZE + j) - Sum minljj = zero For i = j + 1 To n For k = 1 To j - 1 tmp1(k) = l((i - 1) * ISIZE + k) tmp2(k) = l((j - 1) * ISIZE + k) Next k Sum = DotProduct(j - 1, tmp1, tmp2) l((i - 1) * ISIZE + j) = h((j - 1) * ISIZE + i) - Sum minljj = XMAX(minljj, Abs(l((i - 1) * ISIZE + j))) Next i minljj = XMAX(minljj / maxffl, minl) If l((j - 1) * ISIZE + j) > minljj * minljj Then ' NORMAL CHOLESKY DECOMPOSITION l((j - 1) * ISIZE + j) = Sqr(l((j - 1) * ISIZE + j)) Else ' IMPLICITLY PERTURB DIAGONAL OF H If minljj < minl2 Then minljj = minl2 maxadd = XMAX(maxadd, minljj * minljj - l((j - 1) * ISIZE + j)) l((j - 1) * ISIZE + j) = minljj End If For k = j + 1 To n l((j + k - 1) * ISIZE + j) = l((j + k - 1) * ISIZE + j) / l((j - 1) * ISIZE + j) Next k Next j End Sub Sub chsolv(overch As Integer, overfl As Integer, maxexp, n, nunit, output As Integer, _ l() As Double, rhs(), s()) '--------------------------------------------------------------------- ' ' FEB. 14, 1991 ' ' THIS SUBROUTINE USES FORWARD/BACKWARD SUBSTITUTION TO SOLVE THE ' SYSTEM OF LINEAR EQUATIONS: ' ' (LL^)S=RHS '---------------------------------------------------------------------- Dim wv2(ISIZE) lsolv overch, overfl, maxexp, n, nunit, output, l, wv2, rhs ltsolv overch, overfl, maxexp, n, nunit, output, l, s, wv2 End Sub Sub condno(overch As Integer, overfl As Integer, maxexp, n, nunit, output As Integer, _ connum, a(), rdiag()) '----------------------------------------------------------------------------- ' N.B. Arguments P, PM & Q have been removed. ' ' FEB. 14, 1991 ' ' THIS SUBROUTINE ESTIMATES THE CONDITION NUMBER OF A QR-DECOMPOSED ' MATRIX USING THE METHOD OF CLINE, MOLER, STEWART AND WILKINSON ' (SIAM J. N.A. 16 P368 (1979) ). ' ' IF A POTENTIAL OVERFLOW IS DETECTED AT ANY POINT THEN A CONDITION NUMBER ' EQUIVALENT TO THAT OF A SINGULAR MATRIX IS ASSIGNED BY THE CALLING ' SUBROUTINE. '----------------------------------------------------------------------------- Dim p(ISIZE), pm(ISIZE), q(ISIZE) overfl = 0 eps = one / (ten ^ maxexp) connum = Abs(rdiag(i)) For j = 2 To n temp = zero For i = 1 To j - 1 If overch <> 0 Then If Abs(a((i - 1) * ISIZE + j)) > ten ^ (maxexp - 1) Then overfl = 1 Exit Sub End If End If temp = temp + Abs(a((i - 1) * ISIZE + j)) Next i temp = temp + Abs(rdiag(j)) connum = XMAX(connum, temp) Next j q(1) = one / rdiag(1) For i = 2 To n If overch <> 0 Then If Log10(Abs(q(1)) + eps) + Log10(Abs(a(i)) + eps) > maxexp Then overfl = 1 Exit Sub End If End If p(i) = a(i) * q(1) Next i For j = 2 To n If overch <> 0 Then If Log10(Abs(p(j)) + eps) - Log10(Abs(rdiag(j)) + eps) > maxexp Then overfl = 1 Exit Sub End If End If qp = (one - p(j)) / rdiag(j) qm = (-one - p(j)) / rdiag(j) temp = Abs(qp) tempm = Abs(qm) For i = j + 1 To n If overch <> 0 Then If Log10(Abs(a((j - 1) * ISIZE + i)) + eps) + Log10(Abs(qm) + eps) > maxexp Then overfl = 1 Exit Sub End If End If pm(i) = p(i) + a((j - 1) * ISIZE + i) * qm If overch <> 0 Then If Log10(Abs(pm(i)) + eps) - Log10(Abs(rdiag(i)) + eps) > maxexp Then overfl = 1 Exit Sub End If End If tempm = tempm + (Abs(pm(i)) / Abs(rdiag(i))) If overch <> 0 Then If tempm > ten ^ (maxexp - 1) Then overfl = 1 Exit Sub End If End If If overch <> 0 Then If Log10(Abs(a((j - 1) * ISIZE + i)) + eps) + Log10(Abs(qp) + eps) > maxexp Then overfl = 1 Exit Sub End If End If p(i) = p(i) + a((j - 1) * ISIZE + i) * qp If overch <> 0 Then If Log10(Abs(p(i)) + eps) - Log(Abs(rdiag(i)) + eps) > maxexp Then overfl = 1 Exit Sub End If End If temp = temp + (Abs(p(i)) / Abs(rdiag(i))) If overch <> 0 Then If temp > ten ^ (maxexp - 1) Then overfl = 1 Exit Sub End If End If Next i If temp >= tempm Then q(j) = qp Else q(j) = qm For k = 1 To n p(j + k) = pm(j + k) Next k End If Next j qnorm = zero For j = 1 To n qnorm = qnorm + Abs(q(j)) If overch <> 0 Then If qnorm > ten ^ (maxexp - 1) Then overfl = 1 Exit Sub End If End If Next j If Log10(connum) - Log10(qnorm) > maxexp Then overfl = 1 Exit Sub End If connum = connum / qnorm rsolv overch, overfl, maxexp, n, nunit, output, a, rdiag, q If overfl <> 0 Then Exit Sub qnorm = zero For j = 1 To n qnorm = qnorm + Abs(q(j)) If overch <> 0 Then If qnorm > ten ^ (maxexp - 1) Then overfl = 1 Exit Sub End If End If Next j connum = connum * qnorm End Sub Sub delcau(cauchy As Integer, overch As Integer, overfl As Integer, itnum, _ maxexp, n, nunit, output As Integer, beta, caulen, delta, epsmch, _ maxstp As Double, newlen As Double, sqrtz, a(), delf(), scalex()) '-------------------------------------------------------------------- ' ' FEB. 23, 1992 ' ' THIS SUBROUTINE ESTABLISHES AN INITIAL TRUST REGION, DELTA, ' IF ONE IS NOT SPECIFIED BY THE USER AND FINDS THE LENGTH OF ' THE SCALED CAUCHY STEP, CAULEN, AT EACH STEP IF THE DOUBLE ' DOGLEG OPTION IS BEING USED. ' ' THE USER HAS TWO CHOICES FOR THE INITIAL TRUST REGION: ' 1) BASED ON THE SCALED CAUCHY STEP ' 2) BASED ON THE SCALED NEWTON STEP '-------------------------------------------------------------------- Const three = 3# Dim wv1(ISIZE) Dim s As String overfl = 0 eps = one / (ten ^ maxexp) ' IF DELTA IS NOT GIVEN EVALUATE IT USING EITHER THE CAUCHY ' STEP OR THE NEWTON STEP AS SPECIFIED BY THE USER. ' THE SCALED CAUCHY LENGTH, CAULEN, IS REQUIRED IN TWO CASES. ' 1) WHEN SELECTED AS THE CRITERION FOR THE INITIAL DELTA ' 2) IN THE DOUBLE DOGLEG STEP REGARDLESS OF (1) If output > 3 Then Line0 Msg " DETERMINATION OF SCALED CAUCHY STEP LENGTH, CAULEN" End If ' FIND FACTOR WHICH GIVES CAUCHY POINT WHEN MULTPLYING ' STEEPEST DESCENT DIRECTION, DELF. ' caulen = ZETA ^ 1.5 / beta ' = SQRTZ^3/BETA For k = 1 To n wv1(k) = delf(k) / scalex(k) Next k twonrm overfl, maxexp, n, epsmch, sqrtz, wv1 If overfl <> 0 Then caulen = ten ^ maxexp If output > 4 Then Line0 s = " ZETA SET TO " + Str$(sqrtz) + " TO AVOID OVERFLOW" Msg s s = " SCALED CAUCHY LENGTH, CAULEN SET TO " + Str$(caulen) + " TO AVOID OVERFLOW" Msg s If itnum = 1 Then Line0 Msg " THE PROBLEM SHOULD BE RESCALED OR A NEW STARTING POINT CHOSEN" Msg " EXECUTION CONTINUES WITH SUBSTITUTIONS AS LISTED ABOVE." End If End If Else If output > 4 Then Line0 s = " SQUARE ROOT OF ZETA, SQRTZ: " + Str$(sqrtz): Msg s End If End If ' NOTE: THE LOWER TRIANGLE OF MATRIX A NOW CONTAINS THE ' TRANSPOSE OF R WHERE A=QR. beta = zero For i = 1 To n temp = zero For j = i To n If overch <> 0 Then If Log10(Abs(a((j - 1) * ISIZE + i)) + eps) + Log10(Abs(delf(j)) + eps) > maxexp Then caulen = Sqr(epsmch) GoTo 120 End If End If temp = temp + a((j - 1) * ISIZE + i) * delf(j) / (scalex(j) * scalex(j)) Next j beta = beta + temp * temp Next i If output > 4 Then Line0 s = " BETA: " + Str$(beta) + " NOTE: CAULEN = ZETA^1.5/BETA": Msg s Line0 End If ' AVOID OVERFLOWS IN FINDING CAULEN If three * Log10(sqrtz + eps) - Log10(beta + eps) < maxexp And overfl = 0 And beta <> zero Then ' NORMAL DETERMINATION caulen = sqrtz * sqrtz * sqrtz / beta ' THIS STEP AVOIDS /IDE BY ZERO IN DOGLEG IN THE (RARE) CASE ' WHERE DELF(I)=0 FOR ALL I BUT THE POINT IS NOT YET A SOLUTION - ' MOST LIKELY A BAD STARTING ESTIMATE. caulen = XMAX(caulen, one / (ten ^ maxexp)) Else ' SUBSTITUTION TO AVOID OVERFLOW caulen = ten ^ maxexp End If 120 If output > 3 Then Line0 s = " SCALED CAUCHY LENGTH, CAULEN: " + Str$(caulen): Msg s End If ' ESTABLISH INITIAL TRUST REGION IF NEEDED If delta < zero Then ' USE DISTANCE TO CAUCHY POINT OR LENGTH OF NEWTON STEP If cauchy <> 0 Then delta = XMIN(caulen, maxstp) Else delta = XMIN(newlen, maxstp) End If If output > 3 Then Line0 s = " INITIAL TRUST REGION SIZE, DELTA: " + Str$(delta): Msg s End If End If End Sub Sub deufls(abort As Integer, deuflh As Integer, geoms As Integer, 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, itnum, jupdm, maxexp, maxlin, n, _ nfunc, nunit, output As Integer, qnupdm As Integer, stopcr As Integer, alpha, confac, delfts, _ epsmch, fcnmax, fcnnew, fcnold, lambda As Double, newmax As Double, sbrnrm, sigma, a(), astore(), _ boundl(), boundu(), delf(), fvec(), hhpi(), jac() As Double, rdiag(), rhs(), s(), sbar(), scalef(), _ scalex(), sn(), xc(), xplus()) '---------------------------------------------------------------------------- ' FEB. 23, 1992 ' ' THIS SUBROUTINE CONDUCTS A LINE SEARCH IN THE NEWTON DIRECTION ' IF NO CONSTRAINTS ARE VIOLATED. IF THE FIRST TRIAL IS A FAILURE ' THERE ARE TWO TYPES OF LINE SEARCH AVAILABLE: ' 1) REDUCE THE RELAXATION FACTOR, LAMBDA, TO SIGMA*LAMBDA ' WHERE SIGMA IS USER-SPECIFIED (GEOMETRIC Line0 SEARCH) ' 2) AT THE FIRST STEP MINIMIZE A QUADRATIC THROUGH THE OBJECTIVE ' FUNCTION AT THE CURRENT POINT AND TRIAL ESTIMATE (WHICH MUST BE A ' FAILURE) WITH INITIAL SLOPE DELFTS. AT SUBSEQUENT STEPS MINIMIZE ' A CUBIC THROUGH THE OBJECTIVE FUNCTION AT THE TWO MOST RECENT ' FAILURES AND THE CURRENT POINT, AGAIN USING THE INITIAL SLOPE, ' DELFTS. ' ' CONVIO INDICATES A CONSTRAINT VIOLATION BY ONE OR MORE COMPONENTS ' FRSTST INDICATES THE FIRST STEP IN THE LINE SEARCH. ' ' RATIO RATIO OF PROPOSED STEP LENGTH IN (I)TH DIRECTION ' TO DISTANCE FROM (I)TH COMPONENT TO BOUNDARY VIOLATED ' RATIOM MINIMUM OF RATIOS FOR ALL CONSTAINTS VIOLATED '---------------------------------------------------------------------------- Const point1 = 0.1, point5 = 0.5, three = 3# Dim lampre As Double Dim lamtmp As Double Dim wv2(ISIZE) Dim convio As Integer Dim frstst As Integer Dim ss As String frstst = 1 overfl = 0 abort = 0 eps = one / (ten ^ maxexp) lambda = one 'added by JPM For k = 1 To maxlin ratiom = one convio = 0 ' FIND TRIAL POINT AND CHECK IF CONSTRAINTS VIOLATED (IF ' CONTYP IS NOT EQUAL TO ZERO). For i = 1 To n ' NOTE: WV2 MARKS VIOLATIONS. WV2(I) CHANGES TO 1 FOR LOWER BOUND ' VIOLATIONS AND TO 2 FOR UPPER BOUND VIOLATIONS. ' CONSTRAINT VIOLATIONS CAN ONLY OCCUR AT THE FIRST STEP. wv2(i) = -one xplus(i) = xc(i) + lambda * sn(i) If contyp > 0 And frstst <> 0 Then If xplus(i) < boundl(i) Then convio = 1 wv2(i) = one ElseIf xplus(i) > boundu(i) Then convio = 1 wv2(i) = two Else wv2(i) = -one End If End If Next i ' IF CONSTRAINTS ARE VIOLATED FIRST REDUCE THE STEP SIZES FOR THE ' VIOLATING COMPONENTS TO OBTAIN A FEASIBLE POINT. ' IF THE DIRECTION TO THIS MODIFIED POINT IS NOT A DESCENT DIRECTION ' OR IF THE MODIFIED STEP DOES NOT LEAD TO AN ACCEPTABLE POINT THEN ' RETURN TO THE NEWTON DIRECTION AND START A LINE SEARCH AT A FEASIBLE ' POINT WHERE THE COMPONENT WHICH HAS THE SMALLEST VALUE OF RATIO ' (DEFINED BELOW) IS TAKEN TO "CONFAC" OF THE DISTANCE TO THE BOUNDARY. ' DEFAULT VALUE OF CONFAC IS 0.95. If convio <> 0 Then If output > 3 Then Line0 ss = " LINE SEARCH STEP: " + Str$(k): Msg ss Line0 ss = " LAMBDA FOR ATTEMPTED STEP: " + Str$(lambda): Msg ss Msg " CONSTRAINT VIOLATED" Line0 Msg " TRIAL ESTIMATES (VIOLATIONS MARKED)." Line0 For i = 1 To n If wv2(i) > zero Then ss = " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)) + " *": Msg ss Else ss = " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)): Msg ss End If Next i End If For i = 1 To n If wv2(i) > zero Then ' FIND RATIO FOR THIS VIOLATING COMPONENT If wv2(i) = one Then ratio = -(xc(i) - boundl(i)) / (xplus(i) - xc(i)) ElseIf wv2(i) = two Then ratio = (boundu(i) - xc(i)) / (xplus(i) - xc(i)) End If ' NOTE: THIS LINE IS FOR OUTPUT PURPOSES ONLY wv2(i) = ratio ratiom = XMIN(ratiom, ratio) If ratio > point5 Then s(i) = confac * ratio * lambda * sn(i) Else ' WITHIN BUFFER ZONE - ONLY TAKE 1/2 ' OF THE STEP YOU WOULD TAKE OTHERWISE. s(i) = point5 * confac * ratio * lambda * sn(i) End If ' ESTABLISH MODIFIED TRIAL POINT xplus(i) = xc(i) + s(i) Else ' FOR NONVIOLATORS XPLUS REMAINS UNCHANGED BUT THE COMPONENT ' OF S IS LOADED TO CHECK THE DIRECTIONAL DERIVATIVE. s(i) = lambda * sn(i) End If Next i If output > 3 Then 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 wv2(i) < zero Then ss = " S(" + Str$(i) + ") = " + Str$(s(i)) + " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)) Msg ss Else ss = " S(" + Str$(i) + ") = " + Str$(s(i)) + " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)) + _ " " + Str$(wv2(i)) Msg ss End If Next i Line0 ss = " MINIMUM OF RATIOS, RATIOM: " + Str$(ratiom): Msg ss End If ' CHECK DIRECTIONAL DERIVATIVE FOR MODIFIED POINT, DLFTSM innerp overch, overfl, maxexp, n, n, n, nunit, output, dlftsm, delf, s overfl = 0 If output > 3 Then Line0 ss = " INNER PRODUCT OF DELF AND S FOR MODIFIED S: " + Str$(dlftsm): Msg ss End If ' IF INNER PRODUCT IS POSITIVE RETURN TO NEWTON DIRECTION ' AND CONDUCT A LINE SEARCH WITHIN THE FEASIBLE REGION. If dlftsm > zero Then If output > 3 Then Msg " DELFTS > 0 START LINE SEARCH AT LAMBDA = CONFAC*LAMBDA*RATIOM" Msg " NOTE: NO TRIAL POINT WAS EVALUATED AT THIS STEP OF LINE SEARCH." End If ' THE STARTING POINT IS SET AT JUST INSIDE THE MOST ' VIOLATED BOUNDARY. lambda = confac * ratiom * lambda ' LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE End If End If 'if convio ' NO CONSTRAINTS VIOLATED - EVALUATE RESIDUAL VECTOR AT NEW POINT fcn1 overfl, n, fvec, xplus nfunc = nfunc + 1 ' CHECK FOR OVERFLOW IN FUNCTION VECTOR EVALUATION. ' IF SO, REDUCE STEP LENGTH AND CONTINUE LINE SEARCH. If overfl <> 0 Then If output > 3 Then _ Msg " OVERFLOW IN FUNCTION VECTOR - STEP LENGTH REDUCED." ' FORCE STEP TO BE WITHIN CONSTRAINTS - DON'T CALL ' THIS THE FIRST STEP, I.E. FRSTST STAYS AT 1 (FOR TRUE). If convio <> 0 Then lambda = ratiom * confac * lambda Else lambda = sigma * lambda End If ' LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE End If ' EVALUATE (POSSIBLY SCALED) OBJECTIVE FUNCTION AT NEW POINT fcnevl overfl, maxexp, n, nunit, output, epsmch, fcnnew, fvec, scalef If overfl <> 0 Then If output > 3 Then Msg " OVERFLOW IN OBJECTIVE FUNCTION - STEP LENGTH REDUCED." End If ' FORCE STEP TO BE WITHIN CONSTRAINTS - DON'T CALL ' THIS THE FIRST STEP, I.E. FRSTST STAYS AT TRUE. If convio <> 0 Then lambda = ratiom * confac * lambda Else lambda = sigma * lambda End If End If ' IF DEUFLHARD'S METHOD IS BEING USED FOR EITHER RELAXATION FACTOR ' INITIALIZATION OR THE SECOND ACCEPTANCE CRITERION THEN EVALUATE SBAR. ' EVALUATION METHOD DEPENDS UPON WHETHER THE JACOBIAN WAS PERTURBED ' IN THE SOLUTION OF THE LINEAR SYSTEM. ' LOGICAL VARIABLE QRSING IS TRUE IF PERTURBATION TOOK PLACE. If deuflh <> 0 Or acptcr = 12 Then If qrsing <> 0 Then ' FORM -J^F AS RIGHT HAND SIDE - METHOD DEPENDS ON WHETHER ' QNUPDM EQUALS 0 OR 1 IF A QUASI-NEWTON UPDATE IS BEING USED. ' IF JUPDM IS 0 THEN THE NEWTON STEP HAS BEEN FOUND IN SUBROUTINE ' NSTPUN. If jupdm = 0 Or qnupdm = 0 Then ' UNSCALED JACOBIAN IN MATRIX JAC For i = 1 To n If sclfch <> 0 Then wv2(i) = -fvec(i) * scalef(i) * scalef(i) Else wv2(i) = -fvec(i) End If Next i mmulv n, jac, wv2, 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 wv2(i) = zero For j = 1 To n wv2(i) = wv2(i) - jac((i - 1) * ISIZE + j) * fvec(j) * scalef(j) Next j Next i rhs(1) = rdiag(1) * wv2(1) For j = 2 To n rhs(j) = zero For i = 1 To j - 1 rhs(j) = rhs(j) + a((i - 1) * ISIZE + j) * wv2(i) Next i rhs(j) = rhs(j) + rdiag(j) * wv2(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 ii = 1 To n sbar(ii) = -fvec(ii) * scalef(ii) Next ii qrsolv overch, overfl, maxexp, n, nunit, output, astore, hhpi, rdiag, sbar Else ' SET UP RIGHT HAND SIDE - MULTIPLY -FVEC BY Q (STORED IN JAC). For ii = 1 To n wv2(ii) = -fvec(ii) * scalef(ii) Next ii avmul n, n, n, n, jac, wv2, 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 If acptcr = 12 Then For ii = 1 To n wv2(ii) = scalex(ii) * sbar(ii) Next ii twonrm overfl, maxexp, n, epsmch, sbrnrm, wv2 End If End If If output > 3 Then Line0 ss = " LINE SEARCH STEP: " + Str$(k): Msg ss Line0 ss = " LAMBDA FOR ATTEMPTED STEP: " + Str$(lambda): Msg ss Line0 If convio = 0 Then Msg " NEW COMPONENT/FCN VECTORS (XPLUS(I) = XC(I) + LAMBDA*SN(I)" Else Msg " NEW FUNCTION VECTORS AT MODIFIED POINT" End If Line0 For i = 1 To n ss = " XPLUS(" + Str$(i) + ") = " + Str$(xplus(i)) + " FVEC(" + Str$(i) + ") = " + Str$(fvec(i)) Msg ss Next i Line0 If sclfch = 0 Then ss = " OBJECTIVE FUNCTION VALUE AT XPLUS: " + Str$(fcnnew): Msg ss Else ss = " SCALED OBJECTIVE FUNCTION VALUE AT XPLUS: " + Str$(fcnnew): Msg ss End If ss = " FCNMAX + ALPHA*LAMBDA*DELFTS: " + Str$(fcnmax + alpha * lambda * delfts): Msg ss If deuflh <> 0 Or acptcr = 12 Then If itnum > 0 Then If sclxch = 0 Then Line0 Msg " DEUFLHARD SBAR VECTOR:" Line0 For i = 1 To n ss = " SBAR(" + Str$(i) + ") = " + Str$(sbar(i)): Msg ss Next i Else Line0 Msg " DEUFLHARD SBAR VECTOR IN SCALE DX UNITS:" Line0 For i = 1 To n ss = " SBAR(" + Str$(i) + ") = " + Str$(scalex(i) * sbar(i)): Msg ss Next i End If End If End If If acptcr = 12 Then Line0 If sclxch = 0 Then ss = " VALUE OF SBRNRM AT XPLUS: ..............." + Str$(sbrnrm): Msg ss Else ss = " VALUE OF SCALED SBRNRM AT XPLUS: ........" + Str$(sbrnrm): Msg ss End If ss = " NEWMAX: ................................." + Str$(newmax): Msg ss End If End If ' *** CHECK FOR ACCEPTABLE STEP *** If fcnnew < fcnmax + alpha * lambda * delfts Then acpcod = 1 ' NOTE: STEP WILL BE ACCEPTED REGARDLESS OF NEXT TEST. ' THIS SECTION IS FOR BOOKKEEPING ONLY. If acptcr = 12 Then If sbrnrm < newmax Then acpcod = 12 End If Exit Sub End If If acptcr = 12 And sbrnrm < newmax Then acpcod = 2 Exit Sub End If ' FAILURE OF STEP ACCEPTANCE TEST If convio <> 0 Then lambda = confac * ratiom * lambda End If ' LAMBDA IS ALREADY SET - SKIP NORMAL PROCEDURE If lambda = zero Then If output > 0 Then Line0 Msg " LAMBDA IS 0.0: NO PROGRESS POSSIBLE - CHECK BOUNDS OR START." End If abort = 1 Exit Sub End If If geoms <> 0 Then ' GEOMETRIC LINE SEARCH lambda = sigma * lambda Else If frstst <> 0 Then frstst = 0 ' FIND MINIMUM OF QUADRATIC AT FIRST STEP lamtmp = -lambda * lambda * delfts / (two * (fcnnew - fcnold - lambda * delfts)) If output > 4 Then Line0 ss = " TEMPORARY LAMBDA FROM QUADRATIC MODEL: " + Str$(lamtmp): Msg ss End If Else ' FIND MINIMUM OF CUBIC AT SUBSEQUENT STEPS FACTOR = one / (lambda - lampre) If lambda * lambda = zero Then lambda = sigma * lambda ' NOTE: IF THIS LAMBDA^2 WAS ZERO ANY SUBSEQUENT ' LAMBDA^2 WILL ALSO BE ZERO. acubic = FACTOR * ((one / lambda * lambda) * (fcnnew - fcnold - lambda * delfts) - ((one / lampre * lambda) _ * (fplpre - fcnold - lampre * delfts))) bcubic = FACTOR * ((-lampre / lambda * lambda) * (fcnnew - fcnold - lambda * delfts) + ((lambda / lampre * lambda) _ * (fplpre - fcnold - lampre * delfts))) If two * Log10(Abs(bcubic) + eps) > maxexp Then lamtmp = sigma * lambda If output > 4 Then Line0 Msg " POTENTIAL OVERFLOW IN CALCULATING TRIAL LAMBDA FROM" Msg " CUBIC MODEL - LAMBDA SET TO SIGMA*LAMBDA." End If Else disc = bcubic * bcubic - three * acubic * delfts If Abs(acubic) <= epsmch Then lamtmp = -delfts / (two * bcubic) Else If disc < zero Then lamtmp = sigma * lambda Else lamtmp = (-bcubic + Sqr(disc)) / (three * acubic) End If End If If output > 4 Then Line0 ss = " TEMPORARY LAMBDA FROM CUBIC MODEL : ....." + Str$(lamtmp): Msg ss End If If lamtmp > sigma * lambda Then lamtmp = sigma * lambda If output > 4 Then Line0 Msg " LAMTMP TOO LARGE - REDUCED TO SIGMA*LAMBDA." End If End If End If '3rd else End If '2nd else End If '1st else lampre = lambda fplpre = fcnnew If lamtmp < point1 * lambda Then If output > 4 Then Line0 Msg " LAMTMP TOO SMALL - INCREASED TO 0.1*LAMBDA." End If lambda = point1 * lambda Else If output > 4 And lamtmp <> sigma * lambda Then Line0 Msg " LAMTMP WITHIN LIMITS - LAMBDA SET TO LAMTMP." End If lambda = lamtmp End If Next k 'end of very long k loop ' FAILURE IN LINE SEARCH acpcod = 0 ' IF A QUASI-NEWTON STEP HAS FAILED IN THE LINE SEARCH THEN SET QNFAIL ' TO TRUE ANS RETURN TO SUBROUTINE NNES. THIS WILL CAUSE THE JACOBIAN ' TO BE RE-EVALUATED EXPLICITLY AND A LINE SEARCH IN THE NEW DIRECTION ' CONDUCTED. If restrt = 0 Then qnfail = 1 Exit Sub End If ' FALL THROUGH MAIN LOOP - WARNING GIVEN If output > 2 And wrnsup = 0 Then Line0 ss = " WARNING: " + Str$(maxlin) + " CYCLES COMPLETED IN LINE SEARCH WITHOUT SUCCESS." Msg ss End If If stopcr = 2 Or stopcr = 3 Then stopcr = 12 Line0 Msg " STOPPING CRITERION RESET FROM 2 TO 12 TO AVOID HANGING." End If End Sub Sub fcnevl(overfl As Integer, maxexp, n, nunit, output As Integer, _ epsmch, fcnnew, fvec(), scalef()) '---------------------------------------------------------------------- ' FEB. 23, 1992 ' ' THE OBJECTIVE FUNCTION, FCNNEW, DEFINED BY: ' ' FCNNEW = 1/2(SCALEF*FVEC^SCALEF*FVEC) ' ' IS EVALUATED BY THIS SUBROUTINE. '---------------------------------------------------------------------- Dim wv1(ISIZE) Dim s As String overfl = 0 eps = one / (ten ^ maxexp) For k = 1 To n wv1(k) = fvec(k) * scalef(k) Next k twonrm overfl, maxexp, n, epsmch, fcnnew, wv1 ' IF AN OVERFLOW WOULD OCCUR SUBSTITUTE A LARGE VALUE FOR FCNNEW If overfl <> 0 Or two * Log10(fcnnew + eps) > maxexp Then overfl = 1 fcnnew = ten ^ maxexp If output > 2 And wrnsup = 0 Then Line0 s = " WARNING: TO AVOID OVERFLOW, OBJECTIVE FUNCTION SET TO: " + Str$(fcnnew): Msg s End If Exit Sub End If fcnnew = fcnnew * fcnnew / two End Sub Sub fordif(overfl As Integer, j, n, deltaj, fvec(), fvecj1(), jacfdm() As Double, xc()) '---------------------------------------------- ' FEB. 6, 1991 '---------------------------------------------- fcn1 overfl, n, fvecj1, xc If overfl = 0 Then For k = 1 To n jacfdm((k - 1) * ISIZE + j) = (fvecj1(k) - fvec(k)) / deltaj Next k End If End Sub Sub innerp(overch As Integer, overfl As Integer, maxexp, ldima, ldimb, n, nunit, _ output As Integer, dtpro, a(), b()) '-------------------------------------------------------------------- ' FEB. 14, 1991 ' ' THIS SUBROUTINE FINDS THE INNER PRODUCT OF TWO VECTORS, A AND B. ' IF OVERCH IS FALSE, UNROLLED LOOPS ARE USED. ' ' LDIMA IS THE DIMENSION OF A ' LDIMB IS THE DIMENSION OF B ' N IS THE DEPTH INTO A AND B THE INNER PRODUCT IS DESIRED. ' ^[USUALLY LDIMA=LDIMB=N) '-------------------------------------------------------------------} Dim s As String eps = one / (ten ^ maxexp) overfl = 0 dtpro = zero If overch <> 0 Then For i = 1 To n If Log10(Abs(a(i)) + eps) + Log10(Abs(b(i)) + eps) > maxexp Then overfl = 1 dtpro = Sign(ten ^ maxexp, a(i)) * Sign(one, b(i)) If output > 2 And wrnsup = 0 Then Line0 s = " WARNING: TO AVOID OVERFLOW, INNER PRODUCT SET TO " + Str$(dtpro): Msg s End If Exit Sub End If dtpro = dtpro + a(i) * b(i) Next i Else ' SET NUMBER OF GROUPS OF EACH SIZE ng32 = n / 32 ng32r = n - 32 * ng32 ng16 = ng32r / 16 ng16r = ng32r - 16 * ng16 ng8 = ng16r / 8 ng8r = ng16r - 8 * ng8 ng4 = ng8r / 4 ng4r = ng8r - 4 * ng4 ' FIND INNER PRODUCT k = 0 If ng32 > 0 Then For kk = 1 To ng32 k = k + 32 dtpro = dtpro + a(k - 31) * b(k - 31) + a(k - 30) * b(k - 30) + _ a(k - 29) * b(k - 29) + a(k - 28) * b(k - 28) + _ a(k - 27) * b(k - 27) + a(k - 26) * b(k - 26) + a(k - 25) * b(k - 25) + a(k - 24) * b(k - 24) dtpro = dtpro + a(k - 23) * b(k - 23) + a(k - 22) * b(k - 22) + _ a(k - 21) * b(k - 21) + a(k - 20) * b(k - 20) + _ a(k - 19) * b(k - 19) + a(k - 18) * b(k - 18) + a(k - 17) * b(k - 17) + a(k - 16) * b(k - 16) dtpro = dtpro + a(k - 15) * b(k - 15) + a(k - 14) * b(k - 14) + _ a(k - 13) * b(k - 13) + a(k - 12) * b(k - 12) + _ a(k - 11) * b(k - 11) + a(k - 10) * b(k - 10) + a(k - 9) * b(k - 9) + a(k - 8) * b(k - 8) dtpro = dtpro + a(k - 7) * b(k - 7) + a(k - 6) * b(k - 6) + a(k - 5) * b(k - 5) + _ a(k - 4) * b(k - 4) + a(k - 3) * b(k - 3) + a(k - 2) * b(k - 2) + a(k - 1) * b(k - 1) + a(k) * b(k) Next kk End If If ng16 > 0 Then For kk = 1 To ng16 k = k + 16 dtpro = dtpro + a(k - 15) * b(k - 15) + a(k - 14) * b(k - 14) + _ a(k - 13) * b(k - 13) + a(k - 12) * b(k - 12) + _ a(k - 11) * b(k - 11) + a(k - 10) * b(k - 10) + a(k - 9) * b(k - 9) + a(k - 8) * b(k - 8) dtpro = dtpro + a(k - 7) * b(k - 7) + a(k - 6) * b(k - 6) + _ a(k - 5) * b(k - 5) + a(k - 4) * b(k - 4) + a(k - 3) * b(k - 3) + _ a(k - 2) * b(k - 2) + a(k - 1) * b(k - 1) + a(k) * b(k) Next kk End If If ng8 > 0 Then For kk = 1 To ng8 k = k + 8 dtpro = dtpro + a(k - 7) * b(k - 7) + a(k - 6) * b(k - 6) + _ a(k - 5) * b(k - 5) + a(k - 4) * b(k - 4) + a(k - 3) * b(k - 3) + _ a(k - 2) * b(k - 2) + a(k - 1) * b(k - 1) + a(k) * b(k) Next kk End If If ng4 > 0 Then For kk = 1 To ng4 k = k + 4 dtpro = dtpro + a(k - 3) * b(k - 3) + a(k - 2) * b(k - 2) + a(k - 1) * b(k - 1) + a(k) * b(k) Next kk End If If ng4r > 0 Then For kk = 1 To ng4r k = k + 1 dtpro = dtpro + a(k) * b(k) Next kk End If End If End Sub Sub jacrot(overfl As Integer, i, maxexp, n, arot, brot, epsmch, a(), jac() As Double) '---------------------------------------------------------- ' FEB. 11, 1991 ' ' JACOBI (OR GIVENS) ROTATION. '---------------------------------------------------------- Dim hold(2) If arot = zero Then c = zero s = -Sign(one, brot) Else hold(1) = arot hold(2) = brot ldhold = 2 twonrm overfl, maxexp, ldhold, epsmch, denom, hold c = arot / denom s = -brot / denom End If For j = i To n y = a((i - 1) * ISIZE + j) w = a(i * ISIZE + j) a((i - 1) * ISIZE + j) = c * y - s * w a(i * ISIZE + j) = s * y + c * w Next j For j = 1 To n y = jac((i - 1) * ISIZE + j) w = jac(i * ISIZE + j) jac((i - 1) * ISIZE + j) = c * y - s * w jac(i * ISIZE + j) = s * y + c * w Next j End Sub Sub lsolv(overch As Integer, overfl As Integer, maxexp, n, nunit, output As Integer, _ l() As Double, b(), rhs()) '----------------------------------------------------------------------------- ' FEB. 14, 1991 ' ' THIS SUBROUTINE SOLVES: ' ' LB = RHS ' ' WHERE L IS TAKEN FROM THE CHOLESKY DECOMPOSITION ' RHS IS A GIVEN RIGHT HAND SIDE WHICH IS NOT OVERWRITTEN ' B IS THE SOLUTION VECTOR ' ' FRSTER IS USED FOR OUTPUT PURPOSES ONLY. '----------------------------------------------------------------------------- 'Label: 30 Dim maxlog As Double Dim frster As Integer Dim s As String frster = 1 overfl = 0 eps = one / (ten ^ maxexp) If overch <> 0 Then If Log10(Abs(rhs(1)) + eps) - Log10(Abs(l(1)) + eps) > maxexp Then overfl = 1 b(i) = Sign(ten ^ maxexp, rhs(1)) * Sign(one, l(1)) If output > 3 Then Line0 s = " WARNING: COMPONENT 1 SET TO " + Str$(b(1)): Msg s End If GoTo 30 End If End If b(1) = rhs(1) / l(1) 30 For i = 2 To n If overch <> 0 Then ' CHECK TO FIND IF ANY TERMS IN THE EVALUATION WOULD OVERFLOW maxlog = Log10(Abs(rhs(i)) + eps) - Log10(Abs(l((i - 1) * ISIZE + i)) + eps) jstar = 0 For j = 1 To i - 1 tmplog = Log10(Abs(l((i - 1) * ISIZE + j)) + eps) + Log10(Abs(b(j)) + eps) - Log10(Abs(l((i - 1) * ISIZE + i)) + eps) If tmplog > maxlog Then jstar = j maxlog = tmplog End If Next j ' IF AN OVERFLOW WOULD OCCUR ASSIGN A VALUE FOR THE ' TERM WITH CORRECT SIGN. If maxlog > maxexp Then overfl = 1 If jstar = 0 Then b(i) = Sign(ten ^ maxexp, rhs(i)) * Sign(one, l((i - 1) * ISIZE + i)) Else b(i) = -Sign(ten ^ maxexp, l((i - 1) * ISIZE + jstar)) * Sign(one, b(jstar)) * Sign(one, l((i - 1) * ISIZE + i)) End If If frster <> 0 Then frster = 0 Line0 End If If output > 3 Then s = " WARNING: COMPONENT " + Str$(i) + " SET TO " + Str$(b(1)): Msg s End If End If End If ' SUM FOR EACH TERM, ORDERING OPERATIONS TO MINIMIZE ' POSSIBILITY OF OVERFLOW. Sum = zero For j = 1 To i - 1 Sum = Sum + (XMIN(Abs(l((i - 1) * ISIZE + j)), Abs(b(j))) / l((i - 1) * ISIZE + i)) * (XMAX(Abs(l((i - 1) * ISIZE + j)), Abs(b(j)))) _ * Sign(one, l((i - 1) * ISIZE + j)) * Sign(one, b(j)) Next j b(i) = rhs(i) / l((i - 1) * ISIZE + i) - Sum Next i End Sub Sub ltsolv(overch As Integer, overfl As Integer, maxexp, n, nunit, output As Integer, _ l() As Double, y(), b()) '-------------------------------------------------------------------------- ' FEB. 14, 1991 ' ' THIS SUBROUTINE SOLVES: ' ' LY = B ' ' WHERE L IS TAKEN FROM THE CHOLESKY DECOMPOSITION ' B IS A GIVEN RIGHT HAND SIDE WHICH IS NOT OVERWRITTEN ' Y IS THE SOLUTION VECTOR ' ' FRSTER IS USED FOR OUTPUT PURPOSES ONLY. '-------------------------------------------------------------------------- 'Label: 30 Dim maxlog As Double Dim frster As Integer Dim s As String frster = 1 overfl = 0 eps = one / (ten ^ maxexp) If overch <> 0 Then If Log10(Abs(b(n)) + eps) - Log10(Abs(l((n - 1) * ISIZE + n)) + eps) > maxexp Then overfl = 1 y(n) = Sign(ten ^ maxexp, b(n)) * Sign(one, l((n - 1) * ISIZE + n)) If output > 3 Then frster = 0 Line0 s = " WARNING: COMPONENT " + Str$(n) + " SET TO " + Str$(y(n)): Msg s End If GoTo 30 End If End If y(n) = b(n) / l((n - 1) * ISIZE + n) 30 For i = n - 1 To 1 Step -1 If overch <> 0 Then ' CHECK TO FIND IF ANY TERMS IN THE EVALUATION WOULD OVERFLOW maxlog = Log10(Abs(b(i)) + eps) - Log10(Abs(l((i - 1) * ISIZE + i)) + eps) jstar = 0 For j = i + 1 To n tmplog = Log10(Abs(l((j - 1) * ISIZE + i)) + eps) + Log10(Abs(y(j)) + eps) - Log10(Abs(l((i - 1) * ISIZE + i)) + eps) If tmplog > maxlog Then jstar = j maxlog = tmplog End If Next j ' IF AN OVERFLOW WOULD OCCUR ASSIGN A VALUE FOR THE ' TERM WITH CORRECT SIGN. If maxlog > maxexp Then overfl = 1 If jstar = 0 Then y(i) = Sign(ten ^ maxexp, b(i)) * Sign(one, l((i - 1) * ISIZE + i)) Else y(i) = -Sign(ten ^ maxexp, l((jstar - 1) * ISIZE + i)) * Sign(one, y(jstar)) * Sign(one, l((i - 1) * ISIZE + i)) End If If frster <> 0 Then frster = 0 Line0 End If If output > 3 Then s = " WARNING: COMPONENT " + Str$(i) + " SET TO " + Str$(y(i)): Msg s End If End If End If ' SUM FOR EACH TERM ORDERING OPERATIONS TO MINIMIZE ' POSSIBILITY OF OVERFLOW. Sum = zero For j = i + 1 To n Sum = Sum + (XMIN(Abs(l((j - 1) * ISIZE + i)), Abs(y(j))) / l((i - 1) * ISIZE + i)) * (XMAX(Abs(l((j - 1) * ISIZE + i)), _ Abs(y(j)))) * Sign(one, l((j - 1) * ISIZE + i)) * Sign(one, y(j)) Next j y(i) = b(i) / l((i - 1) * ISIZE + i) - Sum Next i End Sub Sub matcop(nradec, nraact, ncadec, ncaact, nrbdec, ncbdec, amat(), bmat()) '------------------------------------------------------------------ ' SEPT. 15, 1991 ' ' COPY A CONTINGUOUS RECTANGULAR PORTION OF ONE MATRIX ' INTO ANOTHER (ELEMENT (1,1) MUST BE INCLUDED). ' ' NRADEC IS 1ST DIMENSION OF AMAT, NRAACT IS LIMIT OF 1ST INDEX ' NCADEC IS 2ND DIMENSION OF AMAT, NCAACT IS LIMIT OF 2ND INDEX ' NRBDEC IS 1ST DIMENSION OF BMAT ' NCBDEC IS 2ND DIMENSION OF BMAT '------------------------------------------------------------------ ' FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32 = nraact / 32 ncc32r = nraact - 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 = 1 To ncaact ' COPY ENTRIES INTO MATRIX B BY COLUMN k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 bmat((k - 32) * ISIZE + j) = amat((k - 32) * ISIZE + j) bmat((k - 31) * ISIZE + j) = amat((k - 31) * ISIZE + j) bmat((k - 30) * ISIZE + j) = amat((k - 30) * ISIZE + j) bmat((k - 29) * ISIZE + j) = amat((k - 29) * ISIZE + j) bmat((k - 28) * ISIZE + j) = amat((k - 28) * ISIZE + j) bmat((k - 27) * ISIZE + j) = amat((k - 27) * ISIZE + j) bmat((k - 26) * ISIZE + j) = amat((k - 26) * ISIZE + j) bmat((k - 25) * ISIZE + j) = amat((k - 25) * ISIZE + j) bmat((k - 24) * ISIZE + j) = amat((k - 24) * ISIZE + j) bmat((k - 23) * ISIZE + j) = amat((k - 23) * ISIZE + j) bmat((k - 22) * ISIZE + j) = amat((k - 22) * ISIZE + j) bmat((k - 21) * ISIZE + j) = amat((k - 21) * ISIZE + j) bmat((k - 20) * ISIZE + j) = amat((k - 20) * ISIZE + j) bmat((k - 19) * ISIZE + j) = amat((k - 19) * ISIZE + j) bmat((k - 18) * ISIZE + j) = amat((k - 18) * ISIZE + j) bmat((k - 17) * ISIZE + j) = amat((k - 17) * ISIZE + j) bmat((k - 16) * ISIZE + j) = amat((k - 16) * ISIZE + j) bmat((k - 15) * ISIZE + j) = amat((k - 15) * ISIZE + j) bmat((k - 14) * ISIZE + j) = amat((k - 14) * ISIZE + j) bmat((k - 13) * ISIZE + j) = amat((k - 13) * ISIZE + j) bmat((k - 12) * ISIZE + j) = amat((k - 12) * ISIZE + j) bmat((k - 11) * ISIZE + j) = amat((k - 11) * ISIZE + j) bmat((k - 10) * ISIZE + j) = amat((k - 10) * ISIZE + j) bmat((k - 9) * ISIZE + j) = amat((k - 9) * ISIZE + j) bmat((k - 8) * ISIZE + j) = amat((k - 8) * ISIZE + j) bmat((k - 7) * ISIZE + j) = amat((k - 7) * ISIZE + j) bmat((k - 6) * ISIZE + j) = amat((k - 6) * ISIZE + j) bmat((k - 5) * ISIZE + j) = amat((k - 5) * ISIZE + j) bmat((k - 4) * ISIZE + j) = amat((k - 4) * ISIZE + j) bmat((k - 3) * ISIZE + j) = amat((k - 3) * ISIZE + j) bmat((k - 2) * ISIZE + j) = amat((k - 2) * ISIZE + j) bmat((k - 1) * ISIZE + j) = amat((k - 1) * ISIZE + j) Next kk End If If ncc16 > 0 Then For kk = 1 To ncc16 k = k + 16 bmat((k - 16) * ISIZE + j) = amat((k - 16) * ISIZE + j) bmat((k - 15) * ISIZE + j) = amat((k - 15) * ISIZE + j) bmat((k - 14) * ISIZE + j) = amat((k - 14) * ISIZE + j) bmat((k - 13) * ISIZE + j) = amat((k - 13) * ISIZE + j) bmat((k - 12) * ISIZE + j) = amat((k - 12) * ISIZE + j) bmat((k - 11) * ISIZE + j) = amat((k - 11) * ISIZE + j) bmat((k - 10) * ISIZE + j) = amat((k - 10) * ISIZE + j) bmat((k - 9) * ISIZE + j) = amat((k - 9) * ISIZE + j) bmat((k - 8) * ISIZE + j) = amat((k - 8) * ISIZE + j) bmat((k - 7) * ISIZE + j) = amat((k - 7) * ISIZE + j) bmat((k - 6) * ISIZE + j) = amat((k - 6) * ISIZE + j) bmat((k - 5) * ISIZE + j) = amat((k - 5) * ISIZE + j) bmat((k - 4) * ISIZE + j) = amat((k - 4) * ISIZE + j) bmat((k - 3) * ISIZE + j) = amat((k - 3) * ISIZE + j) bmat((k - 2) * ISIZE + j) = amat((k - 2) * ISIZE + j) bmat((k - 1) * ISIZE + j) = amat((k - 1) * ISIZE + j) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 bmat((k - 8) * ISIZE + j) = amat((k - 8) * ISIZE + j) bmat((k - 7) * ISIZE + j) = amat((k - 7) * ISIZE + j) bmat((k - 6) * ISIZE + j) = amat((k - 6) * ISIZE + j) bmat((k - 5) * ISIZE + j) = amat((k - 5) * ISIZE + j) bmat((k - 4) * ISIZE + j) = amat((k - 4) * ISIZE + j) bmat((k - 3) * ISIZE + j) = amat((k - 3) * ISIZE + j) bmat((k - 2) * ISIZE + j) = amat((k - 2) * ISIZE + j) bmat((k - 1) * ISIZE + j) = amat((k - 1) * ISIZE + j) Next kk End If If ncc4 > 0 Then For kk = 1 To ncc4 k = k + 4 bmat((k - 4) * ISIZE + j) = amat((k - 4) * ISIZE + j) bmat((k - 3) * ISIZE + j) = amat((k - 3) * ISIZE + j) bmat((k - 2) * ISIZE + j) = amat((k - 2) * ISIZE + j) bmat((k - 1) * ISIZE + j) = amat((k - 1) * ISIZE + j) Next kk End If If ncc4r > 0 Then For kk = 1 To ncc4r k = k + 1 bmat((k - 1) * ISIZE + j) = amat((k - 1) * ISIZE + j) Next kk End If Next j End Sub Sub qrsolv(overch As Integer, overfl As Integer, maxexp, n, nunit, output As Integer, _ a(), hhpi(), rdiag(), b()) '--------------------------------------------------------------------------- ' FEB. 2, 1991 ' ' THIS SUBROUTINE SOLVES ' ' ^[QR]X=B ' ' WHERE Q AND R ARE OBTAINED FROM THE QR DECOMPOSITION ' B IS A GIVEN RIGHT HAND SIDE WHICH IS OVERWRITTEN ' ' R IS CONTAINED IN THE STRICT UPPER TRIANGLE OF ' MATRIX A AND THE VECTOR RDIAG ' Q IS "CONTAINED" IN THE LOWER TRIANGLE OF MATRIX A ' FRSTOV INDICATES FIRST OVERFLOW - USED ONLY TO SET BORDER FOR OUTPUT '--------------------------------------------------------------------------- Dim frstov As Integer Dim s As String eps = one / (ten ^ maxexp) frstov = 1 overfl = 0 ' MULTIPLY RIGHT HAND SIDE BY Q THEN SOLVE USING R STORED IN MATRIX A For j = 1 To n - 1 tau = zero For i = j To n If overch <> 0 Then If Log10(Abs(a((i - 1) * ISIZE + j)) + eps) + Log10(Abs(b(i)) + eps) - Log10(hhpi(j) + eps) > maxexp Then overfl = 1 tau = Sign(ten ^ maxexp, a((i - 1) * ISIZE + j)) * Sign(one, b(j)) End If End If tau = tau + a((i - 1) * ISIZE + j) * b(i) / hhpi(j) Next i For i = j To n If overch <> 0 Then If Log10(Abs(tau) + eps) + Log10(Abs(a((i - 1) * ISIZE + j)) + eps) > maxexp Then overfl = 1 b(i) = -Sign(ten ^ maxexp, tau) * Sign(one, a((i - 1) * ISIZE + j)) If output > 2 And wrnsup = 0 Then If frstov <> 0 Then Line0 s = " WARNING: COMPONENT " + Str$(i) + " SET TO " + Str$(b(i)) + " IN QRSOLV BEFORE RSOLV": Msg s End If End If End If b(i) = b(i) - tau * a((i - 1) * ISIZE + j) Next i Next j rsolv overch, overfl, maxexp, n, nunit, output, a, rdiag, b End Sub Sub qrupda(overfl As Integer, maxexp, n, epsmch, a(), jac() As Double, u(), V()) '----------------------------------------------------------------- ' FEB. 12, 1991 ' ' UPDATE QR DECOMPOSITION USING A SERIES OF GIVENS ROTATIONS. '----------------------------------------------------------------- Dim hold(2) ' REPLACE SUBDIAGONAL WITH ZEROS SO THAT WHEN R IS MULTIPLIED BY GIVENS ' (JACOBI) ROTATIONS THE SUBDIAGONAL ELEMENTS DO NOT AFFECT THE OUTCOME For i = 2 To n a((i - 1) * ISIZE + i - 1) = zero Next i ' FIND LARGEST K FOR WHICH U(K) DOES NOT EQUAL ZERO k = n For l = 1 To n If u(k) = zero Then If k > 1 Then k = k - 1 Else Exit Sub End If Else Exit Sub End If Next l ' MULTIPLY UV^ BY A SERIES OF ROTATIONS SO THAT ALL BUT THE ' TOP ROW IS MADE ZERO (THEORETICALLY THIS IS WHAT HAPPENS ' ALTHOUGH THIS MATRIX ISN'T ACTUALLY FORMED). For i = k - 1 To 1 Step -1 jacrot overfl, i, maxexp, n, u(i), u(i + 1), epsmch, a, jac If u(i) = zero Then ' THIS STEP JUST AVOIDS ADDING ZERO u(i) = Abs(u(i + 1)) Else hold(1) = u(i) hold(2) = u(i + 1) ldhold = 2 twonrm overfl, maxexp, ldhold, epsmch, eucnrm, hold u(i) = eucnrm End If Next i ' ADD THE TOP ROW TO THE TOP ROW OF A - THIS FORMS THE ' UPPER HESSENBERG MATRIX. For i = 1 To n a(i) = a(i) + u(1) * V(i) Next i ' FORM THE UPPER TRIANGULAR R MATRIX BY A SERIES OF ROTATIONS ' TO ZERO OUT THE SUBDIAGONALS. For i = 1 To k - 1 jacrot overfl, i, maxexp, n, a((i - 1) * ISIZE + i), a(i * ISIZE + i), epsmch, a, jac Next i End Sub Sub rsolv(overch As Integer, overfl As Integer, maxexp, n, nunit, output As Integer, _ a(), rdiag(), b()) '-------------------------------------------------------------------------- ' FEB. 14, 1991 ' ' THIS SUBROUTINE SOLVES, BY BACKWARDS SUBSTITUTION, ' ' RX = B ' ' WHERE R IS TAKEN FROM THE QR DECOMPOSITION AND IS STORED IN ' THE STRICT UPPER TRIANGLE OF MATRIX A AND THE VECTOR, ' RDIAG ' B IS A GIVEN RIGHT HAND SIDE WHICH IS OVERWRITTEN ' ' FRSTOV INDICATES FIRST OVERFLOW - USED ONLY TO SET BORDERS FOR OUTPUT '-------------------------------------------------------------------------- 'Label: 30 Dim maxlog As Double Dim frstov As Integer Dim s As String frstov = 1 overfl = 0 eps = one / (ten ^ maxexp) If overch <> 0 Then If Log10(Abs(b(n)) + eps) - Log10(Abs(rdiag(n)) + eps) > maxexp Then overfl = 1 b(n) = Sign(ten ^ maxexp, b(n)) * Sign(one, rdiag(n)) If output > 2 And wrnsup = 0 Then frstov = 0 Line0 s = " WARNING: COMPONENT " + Str$(n) + " SET TO " + Str$(b(n)): Msg s End If GoTo 30 End If End If b(n) = b(n) / rdiag(n) 30: For i = n - 1 To 1 Step -1 If overch <> 0 Then ' CHECK TO FIND IF ANY TERMS IN THE EVALUATION WOULD OVERFLOW maxlog = Log10(Abs(b(i)) + eps) - Log10(Abs(rdiag(i)) + eps) jstar = 0 For j = i + 1 To n tmplog = Log10(Abs(a((i - 1) * ISIZE + j)) + eps) + Log10(Abs(b(j)) + eps) - Log10(Abs(rdiag(i)) + eps) If tmplog > maxlog Then jstar = j maxlog = tmplog End If Next j ' IF AN OVERFLOW WOULD OCCUR ASSIGN A VALUE FOR THE ' TERM WITH CORRECT SIGN. If maxlog > maxexp Then overfl = 1 If jstar = 0 Then b(i) = Sign(ten ^ maxexp, b(i)) * Sign(one, rdiag(i)) Else b(i) = -Sign(ten ^ maxexp, a((i - 1) * ISIZE + jstar)) * Sign(one, b(jstar)) * Sign(one, rdiag(i)) End If If frstov <> 0 Then frstov = 0 Line0 End If If output > 2 And wrnsup = 0 Then s = " WARNING: COMPONENT " + Str$(i) + " SET TO " + Str$(b(i)): Msg s End If End If End If ' SUM FOR EACH TERM ORDERING OPERATIONS TO MINIMIZE ' POSSIBILITY OF OVERFLOW. Sum = zero For j = i + 1 To n Sum = Sum + (XMIN(Abs(a((i - 1) * ISIZE + j)), Abs(b(j))) / rdiag(i)) * (XMAX(Abs(a((i - 1) * ISIZE + j)), Abs(b(j)))) _ * Sign(one, a((i - 1) * ISIZE + j)) * Sign(one, b(j)) Next j b(i) = b(i) / rdiag(i) - Sum Next i End Sub Sub setup(absnew As Integer, cauchy As Integer, deuflh As Integer, geoms As Integer, linesr, newton, _ overch As Integer, acptcr As Integer, itsclf, itsclx, jactyp, jupdm, maxexp, maxit, maxns, _ maxqns, minqns, n, narmij, niejev, njacch, output As Integer, qnupdm As Integer, stopcr _ As Integer, supprs As Integer, trupdm As Integer, alpha, confac, delta, delfac, epsmch, etafac, _ fdtolj, ftol, lam0 As Double, mstpf As Double, nsttol As Double, omega, ratiof, sigma, stptol, _ boundl(), boundu(), scalef(), scalex(), help As String) '------------------------------------------------------------------------ ' DEC. 7, 1991 ' ' SUBROUTINE SETUP ASSIGNS DEFAULT VALUES TO ALL REQUISITE PARAMETERS. '------------------------------------------------------------------------ ' LOGICAL VALUES absnew = 0 bypass = 0 cauchy = 0 deuflh = 1 geoms = 1 linesr = 1 matsup = 0 newton = 0 overch = 0 wrnsup = 0 ' INTEGER VALUES acptcr = 12 itsclf = 0 itsclx = 0 jactyp = 1 jupdm = 0 maxit = 100 'max. number of iterations maxns = 50 maxqns = 10 minqns = 7 narmij = 1 nfetot = 0 niejev = 1 njacch = 1 output = 2 'set to 4 to have more outputs qnupdm = 1 stopcr = 12 supprs = 0 trupdm = 0 ' REAL VALUES alpha = 0.0001 confac = 0.95 delta = -1# delfac = 2# etafac = 0.2 lam0 = 1# mstpf = 1000# omega = 0.1 ratiof = 0.7 sigma = 0.5 ' CHARACTER VARIABLE help = "NONE" ' NOTE: NOTATIONAL CHANGES IN CALLING PROGRAM FROM MACHAR ' 1) EPSMCH DENOTES MACHINE EPSILON ' 2) MINEBB DENOTES MINIMUM EXPONENT BASE BETA ' 3) MAXEBB DENOTES MAXIMUM EXPONENT BASE BETA 'Values extracted from Fortan 90 version IT = 53 ibeta = 2 minebb = -1021 maxebb = 1024 epsmch = 2.22044604925031E-16 XMAX1 = 1.7976931348E+300 maxexp = 308 ' VALUES FOR TWO-NORM CALCULATIONS smallb = (one * ibeta) ^ (minebb + 1) / 2 bigb = (one * ibeta) ^ (maxebb - IT + 1) / 2 smalls = (one * ibeta) ^ (minebb - 1) / 2 'bigs = (one * ibeta) ^ (maxebb + IT - 1) / 2 bigs = 8.997827589E+160 bigr = XMAX1 ' SET STOPPING CRITERIA PARAMETERS fdtolj = 0.000001 ftol = epsmch ^ 0.333 nsttol = ftol * ftol stptol = nsttol ' VECTOR VALUES temp = -(ten ^ maxexp) For i = 1 To n boundl(i) = temp boundu(i) = -temp scalef(i) = one scalex(i) = one Next i End Sub Sub twonrm(overfl As Integer, maxexp, n, epsmch, eucnrm, V()) '---------------------------------------------------------------------- ' FEB. 23 ,1992 ' ' THIS SUBROUTINE EVALUATES THE EUCLIDEAN NORM OF A VECTOR, V. ' IT FOLLOWS THE ALGORITHM OF J.L. BLUE IN ACM TOMS V4 15 (1978) ' BUT USES SLIGHTLY DIFFERENT CUTS. THE CONSTANTS IN COMMON BLOCK ' NNES_5 ARE CALCULATED IN THE SUBROUTINE MACHAR OR ARE PROVIDED ' BY THE USER IN THE DRIVER. '---------------------------------------------------------------------- overfl = 0 sqrtep = Sqr(epsmch) asmall = zero amed = zero abig = zero For i = 1 To n ' ACCUMULATE SUMS OF SQUARES IN ONE OF THREE ACCULULATORS, ' ABIG, AMED AND ASMALL, DEPENDING ON THEIR SIZES. absvi = Abs(V(i)) ' THIS COMPARISON RESTRICTS THE MAXIMUM VALUE OF AMED TO BE ' B/N => CANNOT SUM SO THAT AMED OVERFLOWS. If absvi > bigb / (one * n) Then ' THIS /ISOR OF 10N RESTRICTS ABIG FROM (PATHALOGICALLY) ' OVERFLOWING FROM SUMMATION. abig = abig + (V(i) / (ten * n * bigs) ^ 2) ElseIf absvi < smalls Then asmall = asmall + (V(i) / smalls) ^ 2 Else amed = amed + V(i) * V(i) End If Next i If abig > zero Then ' IF OVERFLOW WOULD OCCUR ASSIGN BIGR AS NORM AND SIGNAL TO ' CALLING SUBROUTINE VIA OVERFL. If Log10(abig) / two + Log10(bigs) + one + Log10(one * n) > maxexp Then eucnrm = bigr overfl = True Exit Sub End If ' IF AMED IS POSITIVE IT COULD CONTRIBUTE TO THE NORM - ' DETERMINATION IS DELAYED UNTIL LATER TO SAVE REPEATING CODE. If amed > zero Then ymin = XMIN(Sqr(amed), ten * n * bigs * Sqr(abig)) ymax = XMAX(Sqr(amed), ten * n * bigs * Sqr(abig)) Else ' AMED DOESN'T CONTRIBUTE AND ASMALL WON'T MATTER IF ' ABIG IS NONZERO - FIND NORM USING ABIG AND RETURN. eucnrm = ten * n * bigs * Sqr(abig) Exit Sub End If ElseIf (asmall > zero) Then If (amed > zero) Then ymin = XMIN(Sqr(amed), smalls * Sqr(asmall)) ymax = XMAX(Sqr(amed), smalls * Sqr(asmall)) Else ' ABIG AND AMED ARE ZERO SO USE ASMALL ONLY. eucnrm = smalls * Sqr(asmall) Exit Sub End If Else eucnrm = Sqr(amed) Exit Sub End If If ymin < sqrtep * ymax Then ' SMALLER PORTION DOES NOT CONTRIBUTE TO NORM. eucnrm = ymax Else ' SMALLER PORTION CONTRIBUTES TO NORM. eucnrm = ymax * Sqr((one + ymin / ymax) * (one + ymin / ymax)) End If End Sub Sub update(mnew, mold, n, trmcod As Integer, fcnnew, fcnold, fvec(), fvecc(), xc(), xplus()) '------------------------------------------------------------------------ ' FEB. 9, 1991 ' ' THIS SUBROUTINE RESETS CURRENT ESTIMATES OF SOLUTION AND UPDATES THE ' OBJECTIVE FUNCTION VALUE, M (USED TO SET HOW MANY PREVIOUS VALUES TO ' LOOK AT IN THE NON- MONOTONIC COMPARISONS) AND THE TERMINATION CODE, ' TRMCOD. '------------------------------------------------------------------------ fcnold = fcnnew mold = mnew trmcod = 0 For i = 1 To n fvecc(i) = fvec(i) xc(i) = xplus(i) Next i End Sub Sub utbmul(ncadec, ncaact, ncbdec, ncbact, ncdec, ncact, amat(), bmat(), cmat()) '--------------------------------------------------------------------------- ' FEB. 8, 1991 ' ' MATRIX MULTIPLICATION: A^B=C WHERE A IS UPPER TRIANGULAR ' ' VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4. ' ' NCADEC IS 2ND DIM. OF AMAT; NCAACT IS ACTUAL LIMIT FOR 2ND INDEX ' NCBDEC IS 2ND DIM. OF BMAT; NCBACT IS ACTUAL LIMIT FOR 2ND INDEX ' NCDEC IS COMMON DIMENSION OF AMAT & BMAT; NCACT IS ACTUAL LIMIT ' ' I.E. NCADEC IS NUMBER OF COLUMNS OF A DECLARED ' NCBDEC IS NUMBER OF COLUMNS OF B DECLARED ' NCDEC IS THE NUMBER OF ROWS IN BOTH A AND B DECLARED ' ' MODIFICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ' QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA '--------------------------------------------------------------------------- For i = 1 To ncaact ' FIND NUMBER OF GROUPS OF SIZE 32, 16... nend = IMIN(i, ncact) ' THIS ADJUSTMENT IS REQUIRED WHEN NCACT IS LESS THAN NCAACT 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 ' FIND ENTRY IN MATRIX C For j = 1 To ncbact Sum = zero k = 0 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 Sum = Sum + amat((k - 32) * ISIZE + i) * bmat((k - 32) * ISIZE + j) + amat((k - 31) * ISIZE + i) * bmat((k - 31) * ISIZE + j) + _ amat((k - 30) * ISIZE + i) * bmat((k - 30) * ISIZE + j) + amat((k - 29) * ISIZE + i) * bmat((k - 29) * ISIZE + j) + _ amat((k - 28) * ISIZE + i) * bmat((k - 28) * ISIZE + j) + amat((k - 27) * ISIZE + i) * bmat((k - 27) * ISIZE + j) + _ amat((k - 26) * ISIZE + i) * bmat((k - 26) * ISIZE + j) + amat((k - 25) * ISIZE + i) * bmat((k - 25) * ISIZE + j) Sum = Sum + amat((k - 24) * ISIZE + i) * bmat(k - 23, j) + amat(k - 22, i) * bmat(k - 22, j) _ + amat((k - 22) * ISIZE + i) * bmat((k - 22) * ISIZE + j) + amat((k - 21) * ISIZE + i) * bmat((k - 21) * ISIZE + j) + _ amat((k - 20) * ISIZE + i) * bmat((k - 20) * ISIZE + j) + amat((k - 19) * ISIZE + i) * bmat((k - 19) * ISIZE + j) + _ amat((k - 18) * ISIZE + i) * bmat((k - 18) * ISIZE + j) + amat((k - 17) * ISIZE + i) * bmat((k - 17) * ISIZE + j) Sum = Sum + amat((k - 16) * ISIZE + i) * bmat((k - 16) * ISIZE + j) + amat((k - 15) * ISIZE + i) * bmat((k - 15) * ISIZE + j) + _ amat((k - 14) * ISIZE + i) * bmat((k - 14) * ISIZE + j) + amat((k - 13) * ISIZE + i) * bmat((k - 13) * ISIZE + j) + _ amat((k - 12) * ISIZE + i) * bmat((k - 12) * ISIZE + j) + amat((k - 11) * ISIZE + i) * bmat((k - 11) * ISIZE + j) + _ amat((k - 10) * ISIZE + i) * bmat((k - 10) * ISIZE + j) + amat((k - 9) * ISIZE + i) * bmat((k - 9) * ISIZE + j) Sum = Sum + amat((k - 8) * ISIZE + i) * bmat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * bmat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * bmat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * bmat((k - 5) * ISIZE + j) + _ amat((k - 4) * ISIZE + i) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * bmat((k - 3) * ISIZE + j) + _ amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 16) * ISIZE + j) + amat((k - 15) * ISIZE + i) * bmat((k - 15) * ISIZE + j) + _ amat((k - 14) * ISIZE + i) * bmat((k - 14) * ISIZE + j) + amat((k - 13) * ISIZE + i) * bmat((k - 13) * ISIZE + j) + _ amat((k - 12) * ISIZE + i) * bmat((k - 12) * ISIZE + j) + amat((k - 11) * ISIZE + i) * bmat((k - 11) * ISIZE + j) + _ amat((k - 10) * ISIZE + i) * bmat((k - 10) * ISIZE + j) + amat((k - 9) * ISIZE + i) * bmat((k - 9) * ISIZE + j) Sum = Sum + amat((k - 8) * ISIZE + i) * bmat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * bmat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * bmat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * bmat((k - 5) * ISIZE + j) + _ amat((k - 4) * ISIZE + i) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * bmat((k - 3) * ISIZE + j) + _ amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 8) * ISIZE + j) + amat((k - 7) * ISIZE + i) * bmat((k - 7) * ISIZE + j) + _ amat((k - 6) * ISIZE + i) * bmat((k - 6) * ISIZE + j) + amat((k - 5) * ISIZE + i) * bmat((k - 5) * ISIZE + j) + _ amat((k - 4) * ISIZE + i) * bmat((k - 4) * ISIZE + j) + amat((k - 3) * ISIZE + i) * bmat((k - 3) * ISIZE + j) + _ amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 4) * ISIZE + j) + _ amat((k - 3) * ISIZE + i) * bmat((k - 3) * ISIZE + j) + _ amat((k - 2) * ISIZE + i) * bmat((k - 2) * ISIZE + j) + _ amat((k - 1) * ISIZE + i) * bmat((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) * bmat((k - 1) * ISIZE + j) Next kk End If cmat((i - 1) * ISIZE + j) = Sum Next j Next i End Sub Sub uvmul(nradec, nraact, ncdec, ncact, amat(), bvec(), cvec()) '------------------------------------------------------------------------------- ' FEB. 8, 1991 ' ' MATRIX-VECTOR MULTIPLICATION: AB=C WHERE A IS UPPER TRIANGULAR ' ' VERSION WITH INNER LOOP UNROLLED TO DEPTHS 32, 16, 8 AND 4 ' EACH ROW OF MATRIX A IS SAVED AS A COLUMN BEFORE USE. ' ' NRADEC IS 1ST DIM. OF AMAT; NRAACT IS ACTUAL LIMIT FOR 1ST INDEX ' NCDEC IS COMMON DIMENSION OF AMAT & BVEC; NCACT IS ACTUAL LIMIT ' ' I.E. NRADEC IS THE NUMBER OF ROWS OF A DECLARED ' NCDEC IS THE COMMON DECLARED DIMENSION (COLUMNS OF A AND ROWS OF B) ' ' MODIFICATION OF THE MATRIX MULTIPLIER DONATED BY PROF. JAMES MACKINNON, ' QUEEN'S UNIVERSITY, KINGSTON, ONTARIO, CANADA '------------------------------------------------------------------------------- For i = 1 To nraact ' FIND NUMBER OF GROUPS OF SIZE 32, 16... ncc32 = (ncact - (i - 1)) / 32 ncc32r = (ncact - (i - 1)) - 32 * ncc32 ncc16 = ncc32r / 16 ncc16r = ncc32r - 16 * ncc16 ncc8 = ncc16r / 8 ncc8r = ncc16r - 8 * ncc8 ncc4 = ncc8r / 4 ncc4r = ncc8r - 4 * ncc4 ' FIND ENTRY FOR VECTOR C Sum = zero k = i - 1 If ncc32 > 0 Then For kk = 1 To ncc32 k = k + 32 Sum = Sum + amat((i - 1) * ISIZE + k - 31) * bvec(k - 31) + amat((i - 1) * ISIZE + k - 30) * bvec(k - 30) + _ amat((i - 1) * ISIZE + k - 29) * bvec(k - 29) + amat((i - 1) * ISIZE + k - 28) * bvec(k - 28) + _ amat((i - 1) * ISIZE + k - 27) * bvec(k - 27) + amat((i - 1) * ISIZE + k - 26) * bvec(k - 26) + _ amat((i - 1) * ISIZE + k - 25) * bvec(k - 25) + amat((i - 1) * ISIZE + k - 24) * bvec(k - 24) Sum = Sum + amat((i - 1) * ISIZE + k - 23) * bvec(k - 23) + amat((i - 1) * ISIZE + k - 22) * bvec(k - 22) + _ amat((i - 1) * ISIZE + k - 21) * bvec(k - 21) + amat((i - 1) * ISIZE + k - 20) * bvec(k - 20) + _ amat((i - 1) * ISIZE + k - 19) * bvec(k - 19) + amat((i - 1) * ISIZE + k - 18) * bvec(k - 18) + _ amat((i - 1) * ISIZE + k - 17) * bvec(k - 17) + amat((i - 1) * ISIZE + k - 16) * bvec(k - 16) Sum = Sum + amat((i - 1) * ISIZE + k - 15) * bvec(k - 15) + amat((i - 1) * ISIZE + k - 14) * bvec(k - 14) + _ amat((i - 1) * ISIZE + k - 13) * bvec(k - 13) + amat((i - 1) * ISIZE + k - 12) * bvec(k - 12) + _ amat((i - 1) * ISIZE + k - 11) * bvec(k - 11) + amat((i - 1) * ISIZE + k - 10) * bvec(k - 10) + _ amat((i - 1) * ISIZE + k - 9) * bvec(k - 9) + amat((i - 1) * ISIZE + k - 8) * bvec(k - 8) Sum = Sum + amat((i - 1) * ISIZE + k - 7) * bvec(k - 7) + amat((i - 1) * ISIZE + k - 6) * bvec(k - 6) + _ amat((i - 1) * ISIZE + k - 5) * bvec(k - 5) + amat((i - 1) * ISIZE + k - 4) * bvec(k - 4) + amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + _ amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc16 > 0 Then For kk = 1 To ncc16 k = k + 16 Sum = Sum + amat((i - 1) * ISIZE + k - 15) * bvec(k - 15) + amat((i - 1) * ISIZE + k - 14) * bvec(k - 14) + _ amat((i - 1) * ISIZE + k - 13) * bvec(k - 13) + amat((i - 1) * ISIZE + k - 12) * bvec(k - 12) + _ amat((i - 1) * ISIZE + k - 11) * bvec(k - 11) + amat((i - 1) * ISIZE + k - 10) * bvec(k - 10) + _ amat((i - 1) * ISIZE + k - 9) * bvec(k - 9) + amat((i - 1) * ISIZE + k - 8) * bvec(k - 8) Sum = Sum + amat((i - 1) * ISIZE + k - 7) * bvec(k - 7) + amat((i - 1) * ISIZE + k - 6) * bvec(k - 6) + _ amat((i - 1) * ISIZE + k - 5) * bvec(k - 5) + amat((i - 1) * ISIZE + k - 4) * bvec(k - 4) + amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + _ amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc8 > 0 Then For kk = 1 To ncc8 k = k + 8 Sum = Sum + amat((i - 1) * ISIZE + k - 7) * bvec(k - 7) + amat((i - 1) * ISIZE + k - 6) * bvec(k - 6) + _ amat((i - 1) * ISIZE + k - 5) * bvec(k - 5) + amat((i - 1) * ISIZE + k - 4) * bvec(k - 4) + amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + _ amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc4 > 0 Then For kk = 1 To ncc4 k = k + 4 Sum = Sum + amat((i - 1) * ISIZE + k - 3) * bvec(k - 3) + amat((i - 1) * ISIZE + k - 2) * bvec(k - 2) + _ amat((i - 1) * ISIZE + k - 1) * bvec(k - 1) + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If If ncc4r > 0 Then For kk = 1 To ncc4r k = k + 1 Sum = Sum + amat((i - 1) * ISIZE + k) * bvec(k) Next kk End If cvec(i) = Sum Next i End Sub