Attribute VB_Name = "modMatrix"
' vbMatrix.bas
'Visual Basic module for matrix operations. Written by Kelly Brumbelow, Dept. of Civil Engineering, Texas A&M University.
Option Explicit
Option Base 1
'NOTE: For ALL routines...
'It is assumed that all arrays are meant to be dimensioned with lower bounds equal to 1: e.g., Dim A(1 to 3, 1 to 4).
' The arrays do not need to have been declared this way, but the zero index is ignored.
Function MatMult(A As Variant, B As Variant) As Variant()
'Multiplies A x B
'Ar is number of rows in A. Ac is number of columns in A.
'Br is number of rows in B. Bc is number of columns in B.
'P is the product of A x B, copied to MatMult
'A check is performed to make sure that the number of columns in A is
' equal to the number of rows in B. If not, sub is left immediately.
Dim Ar, Ac, Br, Bc, X, Y, i, P(), StatusString As String
On Error GoTo ErrHandler
StatusString = "Finding Ar"
Ar = UBound(A, 1)
StatusString = "Found Ar"
Ac = UBound(A, 2)
Br = UBound(B, 1)
StatusString = "Finding Bc"
Bc = UBound(B, 2)
StatusString = "Found Bc"
If Ac <> Br Then Exit Function
If Ar > 0 And Bc > 0 Then ' Both A and B are 2-D Matrices; Product is 2-D matrix
ReDim P(1 To Ar, 1 To Bc)
For X = 1 To Ar
For Y = 1 To Bc
P(X, Y) = 0
Next Y
Next X
For X = 1 To Ar
For Y = 1 To Bc
For i = 1 To Ac
P(X, Y) = P(X, Y) + (A(X, i) * B(i, Y))
Next i
Next Y
Next X
ElseIf Ar = 0 Then ' A is a row vector, B is a 2-D matrix; Product is a row vector
ReDim P(1 To Bc)
For Y = 1 To Bc
P(Y) = 0
Next Y
For Y = 1 To Bc
For i = 1 To Ac
P(Y) = P(Y) + (A(i) * B(i, Y))
Next i
Next Y
ElseIf Bc = 0 Then ' A is a 2-D matrix, B is a column vector; Product is a column vector
ReDim P(1 To Ar)
For X = 1 To Ar
P(X) = 0
Next X
For X = 1 To Ar
For i = 1 To Ac
P(X) = P(X) + (A(X, i) * B(i))
Next i
Next X
End If
MatMult = P
Exit Function
ErrHandler:
Select Case StatusString
Case "Finding Ar"
Ar = 0
Case "Finding Bc"
Bc = 0
End Select
Resume Next
End Function
Function MatInv(A As Variant) As Variant()
'A is the square matrix to be inverted
'r is the dimension of the matrix A (r by r)
'Ainv is the inverse matrix
'Dim A0(1000, 2000), I0(1000, 1000), Ix(1000, 1000), Ax(1000, 2000), Axp1(1000, 2000)
Dim A0(), I0(), Ix(), Ax(), Axp1(), Ainv()
Dim i, j, X, r
r = UBound(A, 1)
j = UBound(A, 2)
If r <> j Then Exit Function
ReDim A0(r, 2 * r), I0(r, r), Ix(r, r), Ax(r, 2 * r), Axp1(r, 2 * r), Ainv(r, r)
'Build A0 matrix (A plus Ident Matrix appended on right side)
For i = 1 To r
For j = 1 To r
A0(i, j) = A(i, j)
Next j
Next i
For i = 1 To r
For j = r + 1 To 2 * r '1 To R
If i = j - r Then
A0(i, j) = 1
Else
A0(i, j) = 0
End If
Next j
Next i
'Define Identity Matrix
For i = 1 To r
For j = 1 To r
If i = j Then
I0(i, j) = 1
Else
I0(i, j) = 0
End If
Next j
Next i
'Perform inversion (Gauss-Jacobs(?) method)
Ax = MatCopy(A0)
For X = 1 To r
Ix = MatCopy(I0)
For i = 1 To r
If i = X Then
Ix(i, X) = 1 / Ax(X, X)
Else
Ix(i, X) = -Ax(i, X) / Ax(X, X)
End If
Next i
Axp1 = MatMult(Ix, Ax)
Ax = MatCopy(Axp1)
Next X
For i = 1 To r
For j = 1 To r
Ainv(i, j) = Ax(i, j + r)
Next j
Next i
MatInv = Ainv
End Function
Function MatCopy(A As Variant) As Variant()
'A is the original matrix
'B is the copied matrix
'r is the number of rows in A
'c is the number of columns in A
Dim r, c, i, j, B(), StatusString As String
On Error GoTo ErrHandler
r = UBound(A, 1)
StatusString = "Finding c"
c = UBound(A, 2)
StatusString = "Found c"
If c > 0 Then
ReDim B(1 To r, 1 To c)
For i = 1 To r
For j = 1 To c
B(i, j) = A(i, j)
Next j
Next i
Else
ReDim B(1 To r)
For i = 1 To r
B(i) = A(i)
Next i
End If
MatCopy = B
Exit Function
ErrHandler:
Select Case StatusString
Case "Finding c"
c = 0
End Select
Resume Next
End Function
Function MatTran(A As Variant) As Variant()
'A is the matrix to be transposed
'AT is the transposed matrix
'Ar is the number of rows in A
'Ac is the number of columns in A
Dim Ar, Ac, i, j, AT(), StatusString As String
On Error GoTo ErrHandler
Ar = UBound(A, 1)
StatusString = "Finding Ac"
Ac = UBound(A, 2)
StatusString = "Found Ac"
If Ac > 0 Then
ReDim AT(1 To Ac, 1 To Ar)
For i = 1 To Ar
For j = 1 To Ac
AT(j, i) = A(i, j)
Next j
Next i
Else
ReDim AT(1 To Ar)
For i = 1 To Ar
AT(i) = A(i)
Next i
End If
MatTran = AT
Exit Function
ErrHandler:
Select Case StatusString
Case "Finding Ac"
Ac = 0
End Select
Resume Next
End Function
Function MatSqRt(A As Variant) As Variant()
'Determines the square root of a symmetric positive definite matrix, A().
'Coded in VB by Kelly Brumbelow. Algorithm by H.Spaeth, 1967 (originally in ALGOL).
'Theta and Eps are input into the subroutine in Spaeth's code, but
' knowing these values ahead of time is difficult.
Dim Theta, Eps
Dim i, j, k
Dim Delta, s, c
Dim BB(1000) 'only R elements are used
Dim Iters
Dim r
Dim Asqrt()
r = UBound(A, 1)
j = UBound(A, 2)
If r <> j Then Exit Function
ReDim Asqrt(1 To r, 1 To r)
c = 0
Theta = 0.999
Eps = 0.0000000000001
Delta = 1
Iters = 0
For i = 1 To r
s = 0
For j = 1 To r
s = s + Abs(A(i, j))
Next j
If c < s Then c = s
Next i
c = 0.5 * Theta / c ^ 0.5
For i = 1 To r
For j = 1 To r
Asqrt(i, j) = 2 * c * A(i, j)
Next j
Next i
Do Until Delta < Eps
Delta = 0
For i = 1 To r
For j = 1 To r
s = 0
For k = 1 To r
s = s - Asqrt(i, k) * Asqrt(k, j)
Next k
BB(j) = Asqrt(i, j) + c * (A(i, j) + s)
Next j
For j = 1 To r
s = Abs(Asqrt(i, j) - BB(j))
If Delta < s Then Delta = s
Asqrt(i, j) = BB(j)
Next j
Next i
For i = 1 To r - 1
For j = i + 1 To r
Asqrt(j, i) = Asqrt(i, j)
Next j
Next i
Iters = Iters + 1: If Iters > 2000 Then Exit Do
Loop
MatSqRt = Asqrt
End Function
Function MatAdd(A As Variant, B As Variant) As Variant()
Dim Ar, Ac, Br, Bc, r, c, Sum(), StatusString As String
On Error GoTo ErrHandler
StatusString = "Finding Ar"
Ar = UBound(A, 1)
StatusString = "Found Ar"
StatusString = "Finding Ac"
Ac = UBound(A, 2)
StatusString = "Found Ac"
StatusString = "Finding Br"
Br = UBound(B, 1)
StatusString = "Found Br"
StatusString = "Finding Bc"
Bc = UBound(B, 2)
StatusString = "Found Bc"
If Ar <> Br Or Ac <> Bc Then Exit Function
If Ar > 0 And Ac > 0 Then
ReDim Sum(1 To Ar, 1 To Ac)
For r = 1 To Ar
For c = 1 To Ac
Sum(r, c) = A(r, c) + B(r, c)
Next c
Next r
ElseIf Ar = 0 Then
ReDim Sum(1 To Ac)
For c = 1 To Ac
Sum(c) = A(c) + B(c)
Next c
ElseIf Ac = 0 Then
ReDim Sum(1 To Ar)
For r = 1 To Ar
Sum(r) = A(r) + B(r)
Next r
End If
MatAdd = Sum
Exit Function
ErrHandler:
Select Case StatusString
Case "Finding Ar"
Ar = 0
Case "Finding Ac"
Ac = 0
Case "Finding Br"
Br = 0
Case "Finding Bc"
Bc = 0
End Select
Resume Next
End Function
Function MatNeg(A As Variant) As Variant()
Dim Ar, Ac, r, c, Neg(), StatusString
On Error GoTo ErrHandler
StatusString = "Finding Ar"
Ar = UBound(A, 1)
StatusString = "Found Ar"
StatusString = "Finding Ac"
Ac = UBound(A, 2)
StatusString = "Found Ac"
If Ac > 0 Then
ReDim Neg(1 To Ar, 1 To Ac)
For r = 1 To Ar
For c = 1 To Ac
Neg(r, c) = A(r, c) * (-1)
Next c
Next r
Else
ReDim Neg(1 To Ar)
For r = 1 To Ar
Neg(r) = A(r) * (-1)
Next r
End If
MatNeg = Neg
Exit Function
ErrHandler:
Select Case StatusString
Case "Finding Ac"
Ac = 0
End Select
Resume Next
End Function