‘——————————————————————————–
Sub SolveAxb(Amatrix() As Double, bvec() As Double, ByRef xvec() As Double, _
N As Integer, iptr As Integer, jptr As Integer, kptr As Integer)
‘
‘ Solve Amatrix(iptr:iptr+N-1,jptr:jptr+N-1)*xvec(kptr:kptr+N-1) = bvec(kptr:kptr+N-1)
‘
Dim wAmatrix() As Double: ReDim wAmatrix(1 To N, 1 To N)
Dim row As Integer, column As Integer
For row = 1 To N
For column = 1 To N: wAmatrix(row, column) = Amatrix(iptr + row – 1, jptr + column – 1): Next column
Next row
Dim wbvec() As Double: ReDim wbvec(1 To N)
For row = 1 To N: wbvec(row) = bvec(kptr + row – 1): Next row
Dim wxvec() As Double: ReDim wxvec(1 To N)
Dim Vmatrix() As Double: ReDim Vmatrix(1 To N, 1 To N)
Dim Wvec() As Double: ReDim Wvec(1 To N)
Call DSVDCMP(wAmatrix, N, N, Wvec, Vmatrix)
Call DSVBKSB(wAmatrix, Wvec, Vmatrix, N, N, wbvec, wxvec)
Dim i As Integer
For i = kptr To kptr + N – 1: xvec(i) = wxvec(i – kptr + 1): Next i
End Sub
‘
‘
‘——————————————————–
‘This is Module xlSVD from SVD02.xls
‘
‘Purpose:
‘ Provide SVD subroutines IAW Numerical Recipies with examples
‘
‘Date: 13 Sep 2003
‘
‘Revisions:
‘13 Sep 2003 Added DMRQMIN, MRQFNC, and SVDFNC
‘12 Jan 2001 New
‘
‘Contents:
‘DSVDCMP(A#(), M%, N%, W#(), V#())
‘DSVBKSB(U#(), W#(), V#(), M%, N%, B#(), X#())
‘Dpythag#(A#, B#)
‘
‘Author(s):
‘T. McCloskey, Clifton Park, New York, 12065
‘
‘**********************************************************
‘
‘———————————————————-
Public Sub DSVDCMP(A() As Double, M As Integer, N As Integer, _
W() As Double, V() As Double)
‘———————————————————-
‘Purpose:
‘ Given a matrix A(1:M, 1:N), this routine computes its
‘ Singular Value Decomposition (SVD),
‘
‘A = U*W*TRANS(V).
‘
‘ Notes:
‘ The matrix U replaces A on output.
‘
‘ The diagonal matrix of singular values, W, is output as a vector W(1:N).
‘ They are unordered
‘
‘ The matrix V (not the transpose TRANS(V)) is output as V(1:N, 1:N).
‘
‘Dates/Revisions:
‘ 7 Jan 2001 – Remove line label dependencies
‘24 Dec 2000 – New
‘
‘USES IMIN, DMAX, DSIGN, Dpythag
‘
‘Externals:
‘From Module ‘SVD02.xForSpt’
‘IMIN%(I%,J%)
‘DMAX#(A#,B#)
‘DSIGN#(A#,B#)
‘
‘From this module
‘Dpythag#(A#, B#)
‘
‘Method:
‘ Householder bidiagonalization and a variant of the QR algoritm
‘ are used.See references.
‘
‘Reference:
‘This Visual/BASIC adaptation is an adaptation of the FORTRAN 77 reference [1],
‘which appears to be a FORTRAN 77 adaptation of the FORTRAN 66 code of
‘refernce [2], which, is an adaptation of the ALGO procedure of reference [3].
‘
‘[1] NUMERICAL RECIPIES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING,
‘Cambridge University Press, 1986-1992
‘
‘[2] Forsythe, G., Malcom, M., Moler, C., COMPUTER METHODS FOR
‘MATHEMATICAL COMPUTATIONS, 1977, Prentice-Hall, Inc., NJ 07632
‘
‘[3] Golub, G.H., and C. Reinsch, “Singular value decomposition and least
‘squares solutions.” in J. H. Wilkinson and C. Reinsch, HANDBOOK FOR
‘AUTOMATIC COMPUTERS, vol II: “Linear Algebra”, Heidelburg: Springer
‘
Dim i As Integer, Its As Integer
Dim j As Integer, JJ As Integer
Dim k As Integer, L As Integer, Nm As Integer
Dim Anorm As Double, C As Double, F As Double, G As Double
Dim H As Double, S As Double, Dscale As Double
Dim X As Double, Y As Double, Z As Double
Dim Rv1() As Double
Dim Atmp As Double, DsclTmp As Double
Dim Iflg1 As Boolean
ReDim Rv1(1 To N)
‘Householder reduction to bidiagonal form
G = 0#
Dscale = 0#
Anorm = 0#
For i = 1 To N
L = i + 1
Rv1(i) = Dscale * G
G = 0#
S = 0#
Dscale = 0#
If i <= M ThenFor k = i To MDscale = Dscale + Abs(A(k, i))Next kIf Dscale <> 0# Then
DsclTmp = 1# / Dscale
For k = i To M
A(k, i) = A(k, i) * DsclTmp
S = S + A(k, i) ^ 2
Next k
F = A(i, i)
G = -DSIGN(Sqr(S), F)
H = F * G – S
A(i, i) = F – G
For j = L To N
S = 0#
For k = i To M
S = S + A(k, i) * A(k, j)
Next k
F = S / H
For k = i To M
A(k, j) = A(k, j) + F * A(k, i)
Next k
Next j
For k = i To M
A(k, i) = Dscale * A(k, i)
Next k
End If
End If
W(i) = Dscale * G
G = 0#
S = 0#
Dscale = 0#
If (i <= M) And (i <> N) Then
For k = L To N
Dscale = Dscale + Abs(A(i, k))
Next k
If Dscale <> 0# Then
DsclTmp = 1# / Dscale
For k = L To N
A(i, k) = A(i, k) * DsclTmp
S = S + A(i, k) ^ 2
Next k
F = A(i, L)
G = -DSIGN(Sqr(S), F)
H = F * G – S
A(i, L) = F – G
For k = L To N
Rv1(k) = A(i, k) / H
Next k
For j = L To M
S = 0#
For k = L To N
S = S + A(j, k) * A(i, k)
Next k
For k = L To N
A(j, k) = A(j, k) + S * Rv1(k)
Next k
Next j
For k = L To N
A(i, k) = Dscale * A(i, k)
Next k
End If
End If
Anorm = DMAX(Anorm, (Abs(W(i)) + Abs(Rv1(i))))
Next i
‘
‘ Accumulation of right-hand transformations.
‘
For i = N To 1 Step -1
If i < N ThenIf G <> 0# Then
For j = L To N‘ Double Division to avoid possible underflow
V(j, i) = (A(i, j) / A(i, L)) / G
Next j
For j = L To N
S = 0#
For k = L To N
S = S + A(i, k) * V(k, j)
Next k
For k = L To N
V(k, j) = V(k, j) + S * V(k, i)
Next k
Next j
End If
For j = L To N
V(i, j) = 0#
V(j, i) = 0#
Next j
End If
V(i, i) = 1#
G = Rv1(i)
L = i
Next i
‘
‘ Accumulation of Left-hand Transformation
‘
For i = IMIN(M, N) To 1 Step -1
L = i + 1
G = W(i)
For j = L To N
A(i, j) = 0#
Next j
If G <> 0# Then
G = 1# / G
For j = L To N
S = 0#
For k = L To M
S = S + A(k, i) * A(k, j)
Next k
F = (S / A(i, i)) * G
For k = i To M
A(k, j) = A(k, j) + F * A(k, i)
Next k
Next j
For j = i To M
A(j, i) = A(j, i) * G
Next j
Else
For j = i To M
A(j, i) = 0#
Next j
End If
A(i, i) = A(i, i) + 1#
Next i
‘
‘ Diagonalization of the Bidiagonal form:Loop over
‘ singular values, and over allowed iterations
‘
For k = N To 1 Step -1
For Its = 1 To 30
Iflg1 = True
For L = k To 1 Step -1
Nm = L – 1
If ((Abs(Rv1(L)) + Anorm) = Anorm) Then
Iflg1 = False
Exit For‘ GoTo SVDCMP002
End If
If ((Abs(W(Nm)) + Anorm) = Anorm) Then
Exit For‘ GoTo SVDCMP001
End If
Next L
‘
If Iflg1 Then
‘SVDCMP001:
C = 0#‘Cancellation of RV1(1), if L > 1
S = 1#
For i = L To k
F = S * Rv1(i)
If ((Abs(F) + Anorm) = Anorm) Then Exit For ‘ GoTo SVDCMP002
G = W(i)
H = Dpythag(F, G)
W(i) = H
H = 1# / H
C = G * H
S = -F * H
For j = 1 To M
Y = A(j, Nm)
Z = A(j, i)
A(j, Nm) = (Y * C) + (Z * S)
A(j, i) = (Z * C) – (Y * S)
Next j
Next i
End If
‘SVDCMP002:
Z = W(k)
If L = k Then
If Z < 0# ThenW(k) = -ZFor j = 1 To NV(j, k) = -V(j, k)Next jEnd IfExit For’ GoTo SVDCMP003End IfIf Its = 30 ThenCall MsgBox(“No Convergence in SVDCMP – aborting”)Erase Rv1()Exit SubEnd IfX = W(L)Nm = k – 1Y = W(Nm)G = Rv1(Nm)H = Rv1(k)F = ((Y – Z) * (Y + Z) + (G – H) * (G + H)) / (2# * H * Y)G = Dpythag(F, 1#)F = ((X – Z) * (X + Z) + H * ((Y / (F + DSIGN(G, F))) – H)) / X” Next QR transformation.’C = 1#S = 1#For j = L To Nmi = j + 1G = Rv1(i)Y = W(i)H = S * GG = C * GZ = Dpythag(F, H)Rv1(j) = ZC = F / ZS = H / ZF = (X * C) + (G * S)G = G * C – X * SH = Y * SY = Y * CFor JJ = 1 To NX = V(JJ, j)Z = V(JJ, i)V(JJ, j) = (Z * S) + (X * C)V(JJ, i) = (Z * C) – (X * S)Next JJZ = Dpythag(F, H)W(j) = Z’Rotation can be arbitrary if z = 0If Z <> 0# Then
Z = 1# / Z
C = F * Z
S = H * Z
End If
F = (C * G) + (S * Y)
X = C * Y – S * G
For JJ = 1 To M
Y = A(JJ, j)
Z = A(JJ, i)
A(JJ, j) = (Z * S) + (Y * C)
A(JJ, i) = (Z * C) – (Y * S)
Next JJ
Next j
Rv1(L) = 0#
Rv1(k) = F
W(k) = X
Next Its
‘SVDCMP003:
Next k
Erase Rv1()
End Sub
‘———————————————————-
Public Sub DSVBKSB(U() As Double, W() As Double, V() As Double, _
M As Integer, N As Integer, B() As Double, X() As Double)
‘———————————————————-
‘Purpose:
‘Solves A*X = B for a vector X, where A is specified by the
‘arrays U, W, V as returned by DSVDCMP.M and N are the
‘logical dimensions of A, and will be equal for square matrices.
‘B(1:M) is the input right-hand side. X(1:n) is the output
‘solution vector. No input quantitie are destroyed, so the
‘routine may be called sequentially with different loads B().
‘
‘Date/Revisions:
‘24 Dec 2000
‘
‘Method:
‘Standard back-substitution. Requires user to “zero” small
‘entries of W(j), ie, after DSVDCMP call, but prior to DSVBKSB
‘call, replace small values of W(j) for zero.One example
‘coding to perform this edit automatically is as follows:
‘
‘Wtol = 2E-16
‘Wmax = 0#
‘For J = 1 To N
‘If W(J) > Wmax Then Wmax = W(J)
‘Next J
‘Wmin = Wmax * WTol
‘For J = 1 To N
‘If W(J) < Wmin Then W(J) = 0#’Next J”References:’NUMERICAL RECIPIES IN FORTRAN 77: THE ART OF SCIENTIFIC’COMPUTING, Cambridge University Press, 1986-1992′———————————————————-Dim i As Integer, j As Integer, JJ As IntegerDim S As Double, tmp() As Double”Dimension Temporary vectorReDim tmp(N)’Calculate TRANS(U)*BFor j = 1 To NS = 0#If W(j) <> 0# Then‘Non-zero result only if w(j) <> 0
For i = 1 To M
S = S + U(i, j) * B(i)
Next i
S = S / W(j)‘This is the divide by w(j)
End If
tmp(j) = S
Next j
‘Matrix multiply by V to get answer.
For j = 1 To N
S = 0#
For JJ = 1 To N
S = S + V(j, JJ) * tmp(JJ)
Next JJ
X(j) = S
Next j
‘Release Dynamic memory
Erase tmp
End Sub
‘———————————————————-
Public Function Dpythag(A As Double, B As Double) As Double
‘———————————————————-
‘Purpose:
‘ Computes SQR(a^2 + b^2) without destructive underflow or overflow.
‘
‘Date/Revisions:
‘24 Dec 2000
‘
‘References:
‘NUMERICAL RECIPIES IN FORTRAN 77: THE ART OF SCIENTIFIC
‘COMPUTING, Cambridge University Press, 1986-1992
Dim AbsA As Double, AbsB As Double
‘
AbsA = Abs(A)
AbsB = Abs(B)
If AbsA > AbsB Then
Dpythag = AbsA * Sqr(1# + (AbsB / AbsA) ^ 2)
Else
If AbsB = 0# Then
Dpythag = 0#
Else
Dpythag = AbsB * Sqr(1# + (AbsA / AbsB) ^ 2)
End If
End If
End Function
‘Attribute VB_Name = “xlForSpt”
‘**********************************************************
‘This is Module xlForSpt from SVD02.xls
‘
‘Purpose:
‘ Provide some FORTRAN like Intrinsics
‘
‘Date: 26 Oct 2000
‘
‘Author(s):
‘T. McCloskey, Sillwater, New York, 12170
‘
‘**********************************************************
‘
‘———————————————————-
Public Function DMAX(x1 As Variant, x2 As Variant) As Variant
‘———————————————————-
‘DMAX = CDbl(WorksheetFunction.Max(x1, x2))
If x1 < x2 Then DMAX = x2 Else DMAX = x1End Function’———————————————————-Public Function DMIN(x1 As Variant, x2 As Variant) As Variant’———————————————————-‘DMIN = CDbl(WorksheetFunction.Min(x1, x2))If x1 < x2 Then DMIN = x1 Else DMIN = x2End Function’———————————————————-Public Function IMAX(iX1 As Integer, iX2 As Integer) As Integer’———————————————————- ‘IMAX = CInt(WorksheetFunction.Max(x1, x2))If iX1 < iX2 Then IMAX = iX2 Else IMAX = iX1End Function’———————————————————-Public Function IMIN(iX1 As Integer, iX2 As Integer) As Integer’———————————————————-‘IMIN = CInt(WorksheetFunction.Min(x1, x2))If iX1 < iX2 Then IMIN = iX1 Else IMIN = iX2End Function’———————————————————-Public Function JMIN(iX1 As Long, iX2 As Long) As Long’———————————————————-‘JMIN = CInt(WorksheetFunction.Min(x1, x2))If iX1 < iX2 Then JMIN = iX1 Else JMIN = iX2End Function’———————————————————-Public Function DSIGN(x1 As Double, x2 As Double) As Double’———————————————————-Dim x3 As Double”Transfer of sign’x3 = Abs(x1)If x2 < 0# ThenDSIGN = -x3ElseDSIGN = x3End IfEnd Function’============================================================Function SciFmt(Mantisa As Integer) As String’============================================================’ Purpose:’To allow construction of a formating string’for use with the Basic FORMAT function.” Date: 1998 Jan 24” Method:’The format string is specific for Scientific’format with one leading digit and a user specified’mantisa length.Non-Negative values are preceded’with a space, negative values are preceded with a’negative sign.For example,”Mantisa% = 3’sFmt$ = SciFmt(Mantisa%)”would yield sFmt = ” 0.000E+000;-0.000E+000″.”This allows uniform formated listings for’positive and negative values.’————————————————-Dim sFmt As StringIf Mantisa > 15 Then Mantisa = 15
sFmt = String$(Mantisa, “0”) & “E+000″
SciFmt = ” 0.” & sFmt & “;” & “-0.” & sFmt
End Function
‘============================================================
Public Function IsEven2(k As Integer) As Boolean
‘============================================================
‘ Purpose:
‘Returns TRUE if integer2 variable is EVEN
‘
‘ Date: 11 May 2002
‘————————————————-
IsEven2 = ((Abs(k) And 1) = 0)
End Function
‘============================================================
Public Function IsOdd2(k As Integer) As Boolean
‘============================================================
‘ Purpose:
‘Returns TRUE if integer2 variable is ODD
‘
‘ Date: 11 May 2002
‘————————————————-
IsOdd2 = ((Abs(k) And 1) = 1)
End Function
Reviews
There are no reviews yet.