'*******************************************************************************************
'*             Inversion of a real square matrix by Householder's method                   *
'* --------------------------------------------------------------------------------------- *
'* SAMPLE RUN:                                                                             *
'*                                                                                         *
'* Input file householder.dat contains:                                                    *
'*                                                                                         *
'* 4                                                                                       *
'* 1.0       0.5        0.333333  0.25                                                     *
'* 0.5       0.333333   0.25      0.2                                                      *
'* 0.333333  0.25       0.2       0.166667                                                 *
'* 0.25      0.2        0.166667  0.142857                                                 *
'*                                                                                         *
'* Output file householder.txt contains:                                                   *
'*                                                                                         *
'* Input square matrix:                                                                    *
'*  1   0.5   0.333333   0.25                                                              *
'*  0.5   0.333333   0.25   0.2                                                            *
'*  0.333333   0.25   0.2   0.166667                                                       *
'*  0.25   0.2   0.166667   0.142857                                                       *
'*                                                                                         *
'* Transformation # 1                                                                      *
'* -1.19315166214903  -0.670492926740837  -0.474932372787676  -0.369835516397145           *
'*  0   0.0664812024380633   0.0657297129201566   0.0586884055729692                       * 
'*  0   0.0720989795265739   0.0771532647936291   0.0724593645897091                       *
'*  0   0.0665741012190317   0.0745318564600783   0.0722012027864846                       *
'* Transformation # 2                                                                      *
'* -1.19315166214903  -0.670492926740837  -0.474932372787676  -0.369835516397145           * 
'*  0  -0.118533219307947  -0.125655520474182  -0.117542173236451                          *
'*  0   0   0.00257161924423549   0.0037833945024712                                       *
'*  0   0   0.00566533362590668   0.00878780895033712                                      *
'* Transformation # 3                                                                      *
'* -1.19315166214903  -0.670492926740837  -0.474932372787676  -0.369835516397145           * 
'*  0  -0.118533219307947  -0.125655520474182  -0.117542173236451                          *
'*  0   0  -0.00622167426262025  -0.00956580449944891                                      * 
'*  0   0   0   0.000187201461739776                                                       *
'*                                                                                         *
'* Inverted matrix:                                                                        *
'*  16.0326813541712  -120.367368862522   240.887472621948  -140.578963337606              *
'* -120.367368862524   1204.1386409279  -2710.00936710887   1686.5344030603                *
'*  240.887472621951  -2710.00936710888   6504.2212289674  -4215.8174559316                * 
'* -140.578963337609   1686.53440306031  -4215.81745593161   2810.33136738203              *
'*                                                                                         *
'* Determinant: -1.64722238571534E-07                                                      *
'*                                                                                         *
'* AP=A*A^-1 matrix:                                                                       *
'*  0.999999999999972   3.41060513164848E-13  -6.82121026329696E-13   3.41060513164848E-13 *
'*  7.105427357601E-15   1   0   0                                                         *
'* -7.105427357601E-15   1.13686837721616E-13   1   1.13686837721616E-13                   *
'* -3.5527136788005E-15   2.8421709430404E-14   0   0.999999999999943                      *
'*                                                                                         * 
'* --------------------------------------------------------------------------------------- *
'* Ref.: "Méthodes de calcul numérique - Tome 2 - By Claude Nowakowski,                    *
'*        PSI Editions, 1984".                                                             *
'*                                                                                         *
'*                                        Visual Basic 2008 Version By J-P Moreau, Paris.  *
'*******************************************************************************************
Module Householder
    Dim n As Integer  'size of matrix A
    Dim A(25, 50), A0(25, 25), V(25) As Double
    Dim D, S, XA, XB, XG, XL As Double
    Dim B(25), AI(25, 25), AP(25, 25), X(25) As Double

    Sub Main()
        'read matrix to be inverted
        'Open input/output text files
        FileOpen(1, "householder.dat", OpenMode.Input)
        FileOpen(2, "householder.txt", OpenMode.Output)
        PrintLine(2)
        PrintLine(2, " Input square matrix:")
        Input(1, n)
        For i = 1 To n
            For j = 1 To n
                Input(1, A(i, j))
                If j = n Then
                    PrintLine(2, "  " + Str$(A(i, j)))
                Else
                    Print(2, "  " + Str$(A(i, j)))
                End If
            Next j
        Next i
        FileClose(1)
        For i = 1 To n
            For j = 1 To n
                A(i, j + n) = 0
                A0(i, j) = A(i, j)
            Next j
            A(i, i + n) = 1
        Next i
        'Transform A into triangular matrix
        S1000()
        'N linear systems to solve
        For K = 1 To n
            For i = 1 To n
                B(i) = A(i, K + n)
            Next i
            'Solve triangular system
            S2000()
            For i = 1 To n
                AI(i, K) = X(i)
            Next i
        Next K
        PrintLine(2)
        PrintLine(2, " Inverted matrix:")
        For i = 1 To n
            For j = 1 To n
                If j = n Then
                    PrintLine(2, "  " + Str$(AI(i, j)))
                Else
                    Print(2, "  " + Str$(AI(i, j)))
                End If
            Next j
        Next i
        'Calculate determinant
        S4000()
        PrintLine(2)
        PrintLine(2, " Determinant: " + Str$(D))
        'Check AP=A*A^-1=Unity Matrix
        PrintLine(2)
        PrintLine(2, " AP=A*A^-1 matrix:")
        S3000()
        For i = 1 To n
            For j = 1 To n
                If j = n Then
                    PrintLine(2, "  " + Str$(AP(i, j)))
                Else
                    Print(2, "  " + Str$(AP(i, j)))
                End If
            Next j
        Next i

        FileClose(2)
    End Sub

    Sub S1000()
        PrintLine(2)
        For K = 1 To n - 1
            S = 0
            For i = K To n
                S = S + A(i, K) * A(i, K)
            Next i
            XL = -Math.Sign(A(K, K)) * Math.Sqrt(S)
            XA = XL * (XL - A(K, K))
            V(K) = A(K, K) - XL
            For i = K + 1 To n
                V(i) = A(i, K)
            Next i
            For j = K + 1 To 2 * n
                XB = 0
                For i = K To n
                    XB = XB + V(i) * A(i, j)
                Next i
                XG = XB / XA
                For i = K To n
                    A(i, j) = A(i, j) - XG * V(i)
                Next i
            Next j
            A(K, K) = XL
            For i = K + 1 To n
                A(i, K) = 0
            Next i
            PrintLine(2, " Transformation #" + Str$(K))
            For i = 1 To n
                For j = 1 To n
                    If j = n Then
                        PrintLine(2, "  " + Str$(A(i, j)))
                    Else
                        Print(2, "  " + Str$(A(i, j)))
                    End If
                Next j
            Next i
        Next K
    End Sub

    Sub S2000()
        Dim j As Integer
        For i = n To 1 Step -1
            j = n : S = 0
2020:       If i = j Then GoTo 2040
            S = S - A(i, j) * X(j)
            j = j - 1
            GoTo 2020
2040:       X(i) = (B(i) + S) / A(i, i)
        Next i
    End Sub

    Sub S3000()
        For i = 1 To n
            For j = 1 To n
                S = 0
                For k = 1 To n
                    S = S + A0(i, k) * AI(k, j)
                Next k
                AP(i, j) = S
            Next j
        Next i
    End Sub

    Sub S4000()
        D = 1
        For i = 1 To n
            D = D * A(i, i)
        Next i
    End Sub

End Module