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