Source Codes

Private Sub Matrix_Inversion(ByVal A(,) As Single, ByVal Bound As Integer, ByRef A_inv(,) As Single)

        '---- Gauss Jordan Inverse Matrix |A I| = |I B| where B = A^-1

        Dim _triang(Bound, (Bound * 2) + 1) As Single

        Dim p, q As Integer

        '----Forming |A I| matrix

        For p = 0 To Bound

            For q = 0 To Bound

                _triang(p, q) = A(p, q)

                If p = q Then

                    _triang(p, q + (Bound + 1)) = 1

                End If

            Next

        Next

        '----Gauss Jordan for reducing to |B I|

        Dim line_1 As Integer

        Dim temp, pivot As Single

        For k = 0 To Bound

            'Bring a non-zero element first by changes lines if necessary

            If _triang(k, k) = 0 Then

                For n = k To Bound

                    If _triang(n, k) <> 0 Then line_1 = n : Exit For 'Finds line_1 with non-zero element

                Next n

                'Change line k with line_1

                For m = k To ((Bound * 2) + 1)

                    temp = _triang(k, m)

                    _triang(k, m) = _triang(line_1, m)

                    _triang(line_1, m) = temp

                Next m

            End If

            pivot = _triang(k, k)

            For n = k To (Bound * 2) + 1

                _triang(k, n) = _triang(k, n) / pivot

            Next n

            'For other lines, make a zero element by using:

            'Ai1=Aij-A11*(Aij/A11)

            'and change all the line using the same formula for other elements

            Dim mult_1 As Single

            For n = 0 To Bound

                If n = k And n = Bound Then Exit For

                If n = k And n < Bound Then n = n + 1 'Do not change that element (already equals to 1), go for next one

                If _triang(n, k) <> 0 Then 'if it is zero, stays as it is

                    mult_1 = _triang(n, k) / _triang(k, k)

                    For m = k To (Bound * 2) + 1

                        _triang(n, m) = _triang(n, m) - (_triang(k, m) * mult_1)

                    Next m

                End If

            Next n

        Next k

        For p = 0 To Bound

            For q = 0 To Bound

                A_inv(p, q) = _triang(p, q + (Bound + 1))

            Next

        Next

    End Sub

    Private Sub Matrix_Transpose(ByVal A(,) As Single, ByVal Bound As Integer, ByRef AT(,) As Single)

        Dim i, j As Integer

        For i = 0 To Bound

            For j = 0 To Bound

                AT(i, j) = A(j, i)

            Next

        Next

    End Sub

Private Sub Choleski_Decomposition(ByVal A(,) As Single, ByVal n As Integer, ByRef L(,) As Single)

        Dim Sum1, Sum2, Sum3 As Single

        ReDim L(n, n)

        L(0, 0) = Math.Sqrt(A(0, 0))

        For j = 1 To n

            L(j, 0) = A(j, 0) / L(0, 0)

        Next

        For i = 1 To n - 1

            For k = 0 To i

                Sum1 = Sum1 + Math.Pow(L(i, k), 2)

            Next

            L(i, i) = Math.Sqrt(A(i, i) - Sum1)

            For j = i + 1 To n

                For k = 0 To i

                    Sum2 = Sum2 + (L(j, k) * L(i, k))

                Next

                L(j, i) = ((A(j, i) - Sum2) / (L(i, i)))

            Next

        Next

        For k = 0 To n - 1

            Sum3 = Sum3 + Math.Pow(L(n, k), 2)

        Next

        L(n, n) = Math.Sqrt(A(n, n) - Sum3)

    End Sub

#Region "Gauss Elimination Method"

    '-----Redefined Gauss Elimination Procedure

    Private Sub Gauss(ByVal A(,) As Single, ByVal B() As Single, ByRef X() As Single, ByVal Bound As Integer)

        '___________________________________________________________________________________________

        '---------------- HouseHolder Transformation developed by Samson Mano ----------------------

        '-------------------- based on the following notes -----------------------------------------

        '----------------- 1. https://en.wikipedia.org/wiki/Gaussian_elimination--------------------

        '___________________________________________________________________________________________

        Dim Triangular_A(Bound, Bound + 1) As Single

        Dim soln(Bound) As Single 'Solution matrix

        For m = 0 To Bound

            For n = 0 To Bound

                Triangular_A(m, n) = A(m, n)

            Next

        Next

        '.... substituting the force to triangularmatrics....

        For n = 0 To Bound

            Triangular_A(n, Bound + 1) = B(n)

        Next

        ForwardSubstitution(Triangular_A, Bound)

        ReverseElimination(Triangular_A, X, Bound)

    End Sub

    Private Sub ForwardSubstitution(ByRef _triang(,) As Single, ByVal bound As Integer)

        'Forward Elimination

        'Dim _fraction As Single

        For k = 0 To bound - 1

            For i = k + 1 To bound

                If _triang(k, k) = 0 Then

                    Continue For

                End If

                _triang(i, k) = _triang(i, k) / _triang(k, k)

                For j = k + 1 To bound + 1

                    _triang(i, j) = _triang(i, j) - (_triang(i, k) * _triang(k, j))

                Next

            Next

        Next

    End Sub

    Private Sub ReverseElimination(ByRef _triang(,) As Single, ByRef X() As Single, ByVal bound As Integer)

        'Back Substitution

        For i = 0 To bound

            X(i) = _triang(i, bound + 1)

        Next

        For i = bound To 0 Step -1

            For j = i + 1 To bound

                X(i) = X(i) - (_triang(i, j) * X(j))

            Next

            X(i) = X(i) / _triang(i, i)

        Next

    End Sub

#End Region

 #Region "EigenValues and EigenVectors Calculator Using HouseHolder Transformation and QR Decomposition"

    '___________________________________________________________________________________________________________________________________________________________________________________________

    '___________________________________________________________________________________________________________________________________________________________________________________________

    Private Sub EigenValues_and_EigenVectors_HouseHolder_Transfromation(ByVal A_Matrix(,) As Double, ByVal n As Integer, ByRef EigenValues() As Double, ByRef Eigen_Vector(,) As Double)

        '___________________________________________________________________________________________

        '---------------- EigenValues and EigenVector calculator developed by Samson Mano ----------

        '------------ based on HouseHolder Transfromation and QR Decomposition ---------------------

        '--------------  Eigen values and Eigen Vectors of A_Matrix is calculated ------------------

        '--------- for 4x4 matrix n = 3 , 7x7 matrix n =4 and ExE matrix n = E-1 !!!! Important Note

        '___________________________________________________________________________________________

        Dim TPPk_Matrix(n, n) As Double

        HouseHolder_Transformation(A_Matrix, n, TPPk_Matrix)

        Dim Possib_EigenVector(n, n) As Double

        Dim lv1, lv2, lv3 As Integer

        '----- keep unit matrix for possib_Eigenvector

        For lv1 = 0 To n Step +1

            Possib_EigenVector(lv1, lv1) = 1

        Next

        Do

            Dim TQMatrix(n, n) As Double

            Dim TRMatrix(n, n) As Double

            QRDecomposition(A_Matrix, n, TQMatrix, TRMatrix)

            '---- Update A matrix to A1 which is A = QR => Q^-1A = R => Q^-1AQ = RQ = A1 and Q^-1 = Q^T

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    A_Matrix(lv1, lv2) = 0

                    For lv3 = 0 To n Step +1

                        A_Matrix(lv1, lv2) = A_Matrix(lv1, lv2) + (TRMatrix(lv1, lv3) * TQMatrix(lv3, lv2))

                    Next

                Next

            Next

            '---- Update EigenVector here which is Q1*Q2*Q3 .... Qk

            Dim Temp_Q(n, n) As Double

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    Temp_Q(lv1, lv2) = 0

                    For lv3 = 0 To n Step +1

                        Temp_Q(lv1, lv2) = Temp_Q(lv1, lv2) + (Possib_EigenVector(lv1, lv3) * TQMatrix(lv3, lv2))

                    Next

                Next

            Next

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    Possib_EigenVector(lv1, lv2) = Temp_Q(lv1, lv2)

                Next

            Next

            '---- Exit if lower triangle of A_Matrix becomes Zero

            If Check_IfLowerTriangleisZero(A_Matrix, n) = True Then

                Exit Do

            End If

        Loop

        '-------------- Fixing the EigenValues here

        For lv1 = 0 To n

            EigenValues(lv1) = A_Matrix(lv1, lv1)

        Next

        '------------ Fixing the EigenVectors here

        For lv1 = 0 To n Step +1

            For lv2 = 0 To n Step +1

                Eigen_Vector(lv1, lv2) = 0

                For lv3 = 0 To n Step +1

                    Eigen_Vector(lv1, lv2) = Eigen_Vector(lv1, lv2) + (TPPk_Matrix(lv1, lv3) * Possib_EigenVector(lv3, lv2))

                Next

            Next

        Next

        '********** Sorting EigenValues / EigenVectors

        For lv1 = 0 To n - 1 Step +1

            For lv2 = lv1 + 1 To n Step +1

                'Sorting Ascending Minimum to Maximum

                If EigenValues(lv1) > EigenValues(lv2) Then

                    SortSwap(EigenValues(lv1), EigenValues(lv2))

                    For lv3 = 0 To n Step +1

                        SortSwap(Eigen_Vector(lv3, lv1), Eigen_Vector(lv3, lv2))

                    Next

                End If

            Next

        Next

    End Sub

    Private Sub SortSwap(ByRef a As Double, ByRef b As Double)

        Dim temp As Double

        temp = a

        a = b

        b = temp

    End Sub

    Private Function Check_IfLowerTriangleisZero(ByRef A_Matrix(,) As Double, ByVal n As Integer) As Boolean

        Dim lv1, lv2 As Integer

        Dim Chk_Bool As Boolean

        Chk_Bool = True

        For lv1 = 1 To n Step +1

            For lv2 = 0 To (lv1 - 1) Step +1

                If Math.Round(A_Matrix(lv1, lv2), 5) <> 0 Then

                    Chk_Bool = False

                    GoTo Ex1

                End If

            Next

        Next

Ex1:

        Return Chk_Bool

    End Function

    Private Sub HouseHolder_Transformation(ByRef A_Matrix(,) As Double, ByVal n As Integer, ByRef PPk_Matrix(,) As Double)

        '___________________________________________________________________________________________

        '---------------- HouseHolder Transformation developed by Samson Mano ----------------------

        '-------------------- based on the following notes -----------------------------------------

        '----------------- 1. https://en.wikipedia.org/wiki/Householder_transformation--------------

        '___________________________________________________________________________________________

        Dim k As Integer

        '--- Fix Unit Matrix for P0 so that we can multiply p1p2..pk

        For k = 0 To n Step +1

            PPk_Matrix(k, k) = 1

        Next

        For k = 0 To n - 2 Step +1

            '----- Find the alpha value alpha = -1*(sgn(a(k+1,k))*sqrt(a(j,k)^2) ->j from k+1 to n

            Dim alpha As Double = 0

            Dim lv1, lv2, lv3 As Integer

            For lv1 = k + 1 To n Step +1

                alpha = alpha + (A_Matrix(lv1, k) ^ 2)

            Next

            alpha = (-1) * (Math.Sign(A_Matrix(k + 1, k))) * Math.Sqrt(alpha)

            '----- Find the r value r_value = math.sqrt((1/2)*(alpha^2 - (A_matrix(k+1,k)*alpha))

            Dim r_value As Double

            r_value = Math.Sqrt((1 / 2) * ((alpha ^ 2) - (A_Matrix(k + 1, k) * alpha)))

            '----- Find the v matrix v(0)=v(1)=v(2)=...=v(k)=0 and v(k+1)=(a(k+1,k)-alpha)/(2*r) and v(j)=a(j,k)/(2*r) -> j = k+2,k+3,...,n

            Dim v_matrix(n) As Double

            For lv1 = 0 To k Step +1

                v_matrix(lv1) = 0

            Next

            v_matrix(k + 1) = (A_Matrix(k + 1, k) - alpha) / (2 * r_value)

            For lv1 = k + 2 To n Step +1

                v_matrix(lv1) = A_Matrix(lv1, k) / (2 * r_value)

            Next

            '----- Find the vvT matrix and P_Matrix

            Dim vvT(n, n) As Double

            Dim I_Matix(n, n) As Double

            'Multiply v matrix with vT matrix

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    vvT(lv1, lv2) = 0

                    I_Matix(lv1, lv2) = 0

                    vvT(lv1, lv2) = v_matrix(lv1) * v_matrix(lv2)

                Next

                I_Matix(lv1, lv1) = 1 '--- All the diagonals are 1

            Next

            'Multiply vvT matrix with 2

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    vvT(lv1, lv2) = 2 * vvT(lv1, lv2)

                Next

            Next

            'Forming P_Matrix here

            Dim P_Matrix(n, n) As Double

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    P_Matrix(lv1, lv2) = I_Matix(lv1, lv2) - vvT(lv1, lv2)

                Next

            Next

            '----- Store P_matrix as mutiples to use in transformation PPk = p1p2p3...pk

            Dim Temp_P(n, n) As Double

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    Temp_P(lv1, lv2) = 0

                    For lv3 = 0 To n Step +1

                        Temp_P(lv1, lv2) = Temp_P(lv1, lv2) + (PPk_Matrix(lv1, lv3) * P_Matrix(lv3, lv2))

                    Next

                Next

            Next

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    PPk_Matrix(lv1, lv2) = Temp_P(lv1, lv2)

                Next

            Next

            '----- Find P_Matrix*A_Matrix*P_Matrix which is the A_Matrix HouseHolder Transformation

            Dim AP_Matrix(n, n) As Double

            'Multiply A_Matrix with P_Matrix to find intermittent AP_Matrix

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    AP_Matrix(lv1, lv2) = 0

                    For lv3 = 0 To n Step +1

                        AP_Matrix(lv1, lv2) = AP_Matrix(lv1, lv2) + (A_Matrix(lv1, lv3) * P_Matrix(lv3, lv2))

                    Next

                Next

            Next

            'Multiply P_Matrix with AP_Matrix

            For lv1 = 0 To n Step +1

                For lv2 = 0 To n Step +1

                    A_Matrix(lv1, lv2) = 0

                    For lv3 = 0 To n Step +1

                        A_Matrix(lv1, lv2) = A_Matrix(lv1, lv2) + (P_Matrix(lv1, lv3) * AP_Matrix(lv3, lv2))

                    Next

                Next

            Next

        Next

    End Sub

    Private Sub QRDecomposition(ByRef A_Matrix(,) As Double, ByRef n As Integer, ByRef Q_Matrix(,) As Double, ByRef R_Matrix(,) As Double)

        '___________________________________________________________________________________________

        '---------------- QR Algorithgm developed by Samson Mano -----------------------------------

        '-------------------- based on the following notes -----------------------------------------

        '----------------- 1. http://ranger.uta.edu/~weems/NOTES4351/hansen.pdf---------------------

        '----------------- 2. http://web.engr.illinois.edu/~heath/scicomp/notes/chap03.pdf----------

        '----------------- 3. http://www.math.usm.edu/lambers/mat610/sum10/lecture9.pdf-------------

        '___________________________________________________________________________________________

        '--------------------Initialization A- Matrix

        Dim A(n, n) As Double

        Dim R(n - 1)(,) As Double

        Dim i, j, k As Integer

        For i = 0 To n

            For j = 0 To n

                A(i, j) = A_Matrix(i, j)

            Next

        Next

        '___________________________________________________________________________

        For k = 0 To (n - 1) Step 1

            Dim A_Norm As Double = 0

            For i = k To n

                A_Norm = A_Norm + (A(i, k) ^ 2) '---- Finding the length of column vector

            Next

            A_Norm = Math.Sqrt(A_Norm)

            Dim Alpha_A As Double = If(A_Norm < 0, A_Norm, -A_Norm)

            '----- Forming V - Matrix v = a - alpha * (e)

            Dim V(n) As Double

            For i = k To n Step 1

                V(i) = A(i, k)

            Next

            V(k) = V(k) - Alpha_A

            '-------------------------

            '--- H = I - (2/V^T*V)* V *V^T

            V_VT_Mul(V, k, n, R)

            For i = k To n Step 1

                Dim V_Ratio As Double

                Dim A_Col(n) As Double

                For j = k To n Step 1

                    A_Col(j) = A(j, i) '---- Transformation for every column

                Next

                V_Ratio = VT_A_Mul(V, A_Col, k, n) * VT_V_Mul(V, k, n)

                For j = k To n Step 1

                    A(j, i) = A_Col(j) - (V_Ratio * V(j)) '---- Transformation for every column

                Next

            Next

        Next

        '------------- Multiplying all the H - Matrix to form Q_Matrix

        Dim Q(0)(,) As Double

        ReDim Q(0)(n, n)

        For i = 0 To n

            Q(0)(i, i) = 1 '----- Creating I - matrix

        Next

        For k = 0 To (n - 1) Step 1

            H_Multiplication(Q, R, k, n, Q)

        Next

        '_______________________________________________________________

        'Transferring the result

        For i = 0 To n

            For j = 0 To n

                R_Matrix(i, j) = A(i, j)

            Next

        Next

        For i = 0 To n

            For j = 0 To n

                Q_Matrix(i, j) = Q(0)(i, j)

            Next

        Next

    End Sub

    Private Sub H_Multiplication(ByRef A()(,) As Double, _

                                 ByRef B()(,) As Double, _

                                 ByRef k As Integer, _

                                 ByRef n As Integer, _

                                 ByRef C()(,) As Double)

        '------- Ordinary Matrix multiplication

        Dim T(n, n) As Double

        For i = 0 To n

            For j = 0 To n

                T(i, j) = 0

                For d = 0 To n

                    T(i, j) = T(i, j) + (A(0)(i, d) * B(k)(d, j))

                Next

            Next

        Next

        ReDim C(0)(n, n)

        For i = 0 To n

            For j = 0 To n

                C(0)(i, j) = T(i, j)

            Next

        Next

    End Sub

    Private Sub V_VT_Mul(ByRef V() As Double, ByRef k As Integer, ByRef n As Integer, ByRef H()(,) As Double)

        ReDim H(k)(n, n)

        Dim c As Double = VT_V_Mul(V, k, n)

        '------- Initilaizing

        For i = 0 To n

            For j = 0 To n

                If i = j Then

                    H(k)(i, j) = 1

                Else

                    H(k)(i, j) = 0

                End If

            Next

        Next

        '---- Doing the operation H = I - (2/vT*V)*(V*VT)

        For i = k To n

            For j = k To n

                If i = j Then

                    H(k)(i, j) = 1 - (c * V(i) * V(j))

                Else

                    H(k)(i, j) = 0 - (c * V(i) * V(j))

                End If

            Next

        Next

    End Sub

    Private Function VT_A_Mul(ByRef V() As Double, ByRef A() As Double, ByRef k As Integer, ByRef n As Integer)

        '-------- Function calculates ratio of |V^T*A/V^T*V|

        Dim VtV As Double = 0 ' Holds the value of V^T*V

        Dim VtA As Double = 0 ' Holds the value of V^T*A

        For i = k To n

            VtA = VtA + (V(i) * A(i))

        Next

        Return (VtA)

    End Function

    Private Function VT_V_Mul(ByRef V() As Double, ByRef k As Integer, ByRef n As Integer)

        '-------- Function calculates ratio of |V^T*A/V^T*V|

        Dim VtV As Double = 0 ' Holds the value of V^T*V

        Dim VtA As Double = 0 ' Holds the value of V^T*A

        For i = k To n

            VtV = VtV + (V(i) ^ 2)

        Next

        Return (2 / VtV)

    End Function

    '___________________________________________________________________________________________________________________________________________________________________________________________

    '___________________________________________________________________________________________________________________________________________________________________________________________

#End Region