Attribute VB_Name = "Module2"
'*****************************************************************
'*         Gauss algorithm for solving linear equations          *
'* ------------------------------------------------------------- *
'* REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges *
'*        and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11].    *
'*                                                               *
'*                    Visual Basic Release By J-P Moreau, Paris. *
'*                                (www.jpmoreau.fr)              *
'*****************************************************************
DefDbl A-H, O-Z
DefInt I-N

Option Base 0

Dim d()
Dim iperm(), xlu()

Dim perm() As Integer
Dim lu() As Double, x()
Dim fs(), fsg()


Sub ISWAP(ia, ib)
Dim tnp As Integer
  tmp = ib: ib = ia: ia = tmp
End Sub

Sub SWAP(a, b)
  tmp = b: b = a: a = tmp
End Sub

' Gauss algorithm for solving linear equations .......................
Sub gauss(mode, n, xmat(), xlumat(), iperm(), b(), x(), signd As Integer, code As Integer)
' integer  mode                      Modus: 0, 1, 2, 3 ...............
' integer  n                         Dimension of matrix .............
' double   xmat(0:n-1,0:n-1),        Input matrix ....................
'          xlumat(0:n-1,0:n-1)       LU decomposition ................
' integer  iperm(0:n-1,0:n-1)        row remutation vector ...........
' double   b(0:n-1),                 right hand side .................
'          x(0:n-1)                  solution of the system ..........
' integer  signd,                    sign of the permutation .........
'          code                      return error code ...............
'Note: in matrices xmat,xlumat, iperm, element [i,j] is replaced by
'      vector element [n*i+j].
 '*====================================================================*
 '*                                                                    *
 '*  The procedure gauss solves a linear system :  mat * x = b.        *
 '*  Here mat is the nonsingular system matrix, b the right hand side  *
 '*  of the system and x the solution vector.                          *
 '*                                                                    *
 '*  gauss uses the Gauss algorithm and computes a triangular factori- *
 '*  zation of mat and scaled column pivot search.  (Crout method with *
 '*  row swaps).                                                       *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   Application:                                                     *
 '*   ============                                                     *
 '*      Solve general linear system with a nonsingular coefficient    *
 '*      matrix.                                                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Control parameter:                                               *
 '*   ==================                                               *
 '*      mode     integer;                                             *
 '*               calling modus for gauss:                             *
 '*       = 0     Find factorization and solve linear system           *
 '*       = 1     Find factorization only.                             *
 '*       = 2     Solve linear system only; the factorization is       *
 '*               already available in lumat. This saves work when     *
 '*               solving a linear system repeatedly for several right *
 '*               hand sides and the same system matrix such as when   *
 '*               inverting the matrix.                                *
 '*       = 3     as under 2, additionally we improve the solution     *
 '*               via iterative refinement (not available here).       *
 '*                                                                    *
 '*   Input parameters:                                                *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of mat and lumat,                          *
 '*               size of the vector b, the right hand side, the       *
 '*               solution x and the permutation vector perm.          *
 '*      xmat     matrix of the linear system. It is stored in vector  *
 '*               form.                                                *
 '*      xlumat   (for mode = 2, 3)                                    *
 '*               LU factors of mat                                    *
 '*               xlumat can be stored in the space of xmat.           *
 '*      iperm    (for mode = 2, 3)                                    *
 '*               Permutation vector, of the row interchangfes in      *
 '*               xlumat.                                              *
 '*      b        Right hand side of the system.                       *
 '*      signd    (for mode = 2, 3)                                    *
 '*               sign of the permutation in perm; the determinant of  *
 '*               mat can be computed as the product of the diagonal   *
 '*               entries of lumat times signd.                        *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   ==================                                               *
 '*      xlumat   (for mode = 0, 1)                                    *
 '*               LU factorization of xmat.                            *
 '*      iperm    (for mode = 0, 1)                                    *
 '*               row ermutation vector                                *
 '*      x        (for mode = 0, 2, 3)                                 *
 '*               solution vector(0:n-1).                              *
 '*      signd    (for mode = 0, 1)                                    *
 '*               sign of perm.                                        *
 '*                                                                    *
 '*   Return value (code):                                             *
 '*   ===================                                              *
 '*      =-1      Max. number (MAXITER) of iterative refinements       *
 '*               reached (MAXITER) while mode = 3                     *
 '*      = 0      all ok                                               *
 '*      = 1      n < 1 or other invalid input                         *
 '*      = 2      lack of memory                                       *
 '*      = 3      Matrix singular                                      *
 '*      = 4      Matrix numerically singular                          *
 '*      = 5      incorrect call                                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   subroutines used:                                                *
 '*   ================                                                 *
 '*                                                                    *
 '*      gaudec: determines  LU decomposition                          *
 '*      gausol: solves the linear system                              *
 '*                                                                    *
 '*====================================================================*
Dim rc As Integer

  If (n < 1) Then
    code = 1
    Return
  End If

  'Select mode

    If mode = 0 Then 'Find factorization and solve system ...................
       Call gaudec(n, xmat, xlumat, iperm, signd, rc)
       If rc = 0 Then
         Call gausol(n, xlumat, iperm, b, x, rc)
         code = rc
       Else
         code = rc
         Exit Sub
       End If

    ElseIf mode = 1 Then 'Find factorization only ...........................
       Call gaudec(n, xmat, xlumat, iperm, signd, rc)
       code = rc
       Exit Sub

    ElseIf mode = 2 Then 'Solve only ........................................
       Call gausol(n, xlumat, iperm, b, x, rc)
       code = rc
       Exit Sub

    ElseIf mode = 3 Then 'Solve and then use iterative refinement ...........
       Form1.Print " fgauss: gausoli not implemented."
       code = 5
       Exit Sub
    End If
  
  code = 5                                               'Wrong call
End Sub

' Gauss decomposition ................................................
Sub gaudec(n, xmat(), xlumat(), iperm(), signd As Integer, rc As Integer)
' integer   n                        size of matrix ..................
' double    xmat(0:n-1,0:n-1),       Input matrix ....................
'           xlumat(0:n-1,0:n-1)      matrix decomposition ............
' integer   iperm(0:n-1),            row interchanges ................
'           signd,                   sign of perm ....................
'           rc                       error code ......................
'Note: in matrices xmat,xlumat, iperm, element [i,j] is replaced by
'      vector element [n*i+j].
 '*====================================================================*
 '*                                                                    *
 '*  gaudec decomposes a nonsingular n x n matrix into a product of a  *
 '*  unit lower and an upper triangular matrix. Both triangular factors*
 '*  are stored in lumat (minus the unit diagonal, of course).         *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   Input parameter:                                                 *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of  mat and lumat,                         *
 '*               size of  b , x and perm.                             *
 '*      xmat     pointer to original system matrix in vector form.    *
 '*                                                                    *
 '*   Output parameters:                                               *
 '*   ==================                                               *
 '*      xlumat   pointer to LU factorization                          *
 '*      iperm    pointer to row permutation vector for xlumat         *
 '*      signd    sign of perm. The determinant of xmat can be computed*
 '*               as the product of the diagonal entries of xlumat     *
 '*               times signd.                                         *
 '*                                                                    *
 '*   Return value (rc):                                               *
 '*   =================                                                *
 '*      = 0      all ok                                               *
 '*      = 1      n < 1 or invalid input                               *
 '*      = 2      lack of memory                                       *
 '*      = 3      Matrix is singular                                   *
 '*      = 4      Matrix numerically singular                          *
 '*                                                                    *
 '*====================================================================*
  Dim MACH_EPS As Double
  ReDim d(n) As Double                  'scaling vector for pivoting

  MACH_EPS = 1.2E-16

  If n < 1 Then
    rc = 1                              'Invalid parameters
    Exit Sub
  End If

  'copy xmat to xlumat
  For i = 0 To n - 1
    For j = 0 To n - 1
      xlumat(n * i + j) = xmat(n * i + j)
    Next j
  Next i

  For i = 0 To n - 1
    iperm(i) = i                         'Initialize iperm
    zmax = 0#
    For j = 0 To n - 1                   'find row maxima
      tmp = Abs(xlumat(n * i + j))
      If tmp > zmax Then zmax = tmp
    Next j

    If zmax = 0# Then                    'xmat is singular
      rc = 3
      Exit Sub
    End If
    d(i) = 1# / zmax
  Next i

  signd = 1                              'initialize sign of iperm

  For i = 0 To n - 1
    piv = Abs(xlumat(n * i + i)) * d(i)
    j0 = i                               'Search for pivot element
    For j = i + 1 To n - 1
      tmp = Abs(xlumat(n * j + i)) * d(j)
      If piv < tmp Then
        piv = tmp                        'Mark pivot element and
        j0 = j                           'its location
      End If
    Next j

    If piv < MACH_EPS Then               'If piv is small, xmat is
      signd = 0                          'nearly singular
      rc = 4
      Exit Sub
    End If

    If j0 <> i Then
      signd = -signd                     'update signd
      Call ISWAP(iperm(j0), iperm(i))    'swap pivot entries
      Call SWAP(d(j0), d(i))             'swap scaling vector
      For j = 0 To n - 1
        'swap j0-th and i-th rows of xlumat
        Call SWAP(xlumat(n * j0 + j), xlumat(n * i + j))
      Next j
    End If

    For j = i + 1 To n - 1               'Gauss elimination
      If xlumat(n * j + i) <> 0# Then
        xlumat(n * j + i) = xlumat(n * j + i) / xlumat(n * i + i)
        tmp = xlumat(n * j + i)
        For m = i + 1 To n - 1
          xlumat(n * j + m) = xlumat(n * j + m) - tmp * xlumat(n * i + m)
        Next m
      End If
    Next j
  Next i

  rc = 0  'all ok

End Sub


' Gauss solution .....................................................
Sub gausol(n, xlumat(), iperm(), b(), x(), rc As Integer)
'   integer n                        size of matrix ..................
'   real*8  xlumat(0:n-1,0:n-1)      decomposed matrix (LU) ..........
'   integer perm(0:n-1)              row permutation vector ..........
'   real*8  b(0:n-1),                Right hand side .................
'           x(0:n-1)                 solution ........................
'   integer rc                       error code ......................
 '====================================================================*
 '                                                                    *
 '  gausol  finds the solution x of the linear system  lumat * x = b  *
 '  for the product matrix lumat, that describes an LU decomposition, *
 '  as produced by gaudec.                                            *
 '                                                                    *
 '====================================================================*
 '                                                                    *
 '   Input parameters:                                                *
 '   ================                                                 *
 '      n        integer (n > 0)                                      *
 '               Dimension of xlumat.                                 *
 '      xlumat   Matrix(0:n-1,0:n-1) stored in a vector(n*n).         *
 '               LU factorization, as produced from gaudec.           *
 '      iperm    integer vector(0:n-1)                                *
 '               row permutation vector for xlumat                    *
 '      b        vector(0:n-1)                                        *
 '               Right hand side of the system.                       *
 '                                                                    *
 '   Output parameter:                                                *
 '   ================                                                 *
 '      x        vector(0:n-1)                                        *
 '               Solution vector                                      *
 '                                                                    *
 '   Return value (rc):                                               *
 '   =================                                                *
 '      = 0      all ok                                               *
 '      = 1      n < 1 or other invalid input parameter               *
 '      = 3      improper LU decomposition ( zero diagonal entry)     *
 '                                                                    *
 '====================================================================*
  Dim MACH_EPS As Double
  
  MACH_EPS = 1.2E-16

  If n < 1 Then
     rc = 1                                   'Invalid input parameter
     Exit Sub
  End If

  MACH_EPS = 1.2E-16

  For k = 0 To n - 1                                         'update b
    x(k) = b(iperm(k))
    For j = 0 To k - 1
      x(k) = x(k) - xlumat(n * k + j) * x(j)
    Next j
  Next k

  For k = n - 1 To 0 Step -1                          'back substitute
    Sum = 0#
    For j = k + 1 To n - 1
      Sum = Sum + xlumat(n * k + j) * x(j)
    Next j

    If Abs(xlumat(n * k + k)) < MACH_EPS Then
      rc = 3
      Exit Sub
    End If
    x(k) = (x(k) - Sum) / xlumat(n * k + k)
  Next k

  rc = 0  'all ok

End Sub


' Determinant  .......................................................
Function det(n, xmat()) As Double
'      integer n                     Dimension of the matrix .........
'      double xmat(0:n-1,0:n-1)       matrix ..........................
 '*====================================================================*
 '*                                                                    *
 '*  det computes the determinant of an n x n real matrix xmat          *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   Input parameter:                                                 *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of xmat                                    *
 '*      xmat     matrix(0:n-1,0:n-1)                                  *
 '*               stored in a vector(n*n).                             *
 '*                                                                    *
 '*   Return value:                                                    *
 '*   =============                                                    *
 '*      REAL     Determinant of mat.                                  *
 '*               If the return value = 0, then the matrix is singular *
 '*               or the storage is insufficient                       *
 '*                                                                    *
 '*====================================================================*
 '*                                                                    *
 '*   subroutine used:                                                 *
 '*   ================                                                 *
 '*                                                                    *
 '*          gaudec ():    LU decomposition of mat                     *
 '*                                                                    *
 '*====================================================================*
 Dim rc As Integer, signd As Integer
 ReDim iperm(n), xlu(n * n)
 Dim MACH_EPS As Double
 Dim MAXROOT As Double

  If n < 1 Then   'n not valid
    det = 0#
    Exit Function
  End If
 
  MACH_EPS = 1.2E-16
  MAXROOT = 1E+16

  Call gaudec(n, xmat, xlu, iperm, signd, rc)              'decompose

  If rc <> 0 Or signd = 0 Then
    det = 0#
    Exit Function
  End If

  tmpdet = 1# * signd

  For i = 0 To n - 1
    If Abs(tmpdet) < MACH_EPS Then
      det = 0#
      Exit Function
    ElseIf Abs(tmpdet) > MAXROOT Or Abs(xlu(n * i + i)) > MAXROOT Then
      det = MAXROOT
      Exit Function
    Else
      tmpdet = tmpdet * xlu(n * i + i)                    'compute det
    End If
  Next i
  
  det = tmpdet

End Function

' Gauss for multiple right hand sides ................................
Sub mgauss(n, k, xmat(), rmat(), code As Integer)
'   integer n,                       Dimension of system .............
'           k                        number of right hand sides ......
'   real*8  mat(0:n-1,0:n-1),        original matrix .................
'           rmat(0:n-1,0:n-1)        Right hand sides/solutions ......
'   integer code                     Error code ......................
 '*====================================================================*
 '*                                                                    *
 '*  mgauss  finds the solution matrix x for the linear system         *
 '*  mat * x = rmat with an  n x n coefficient matrix mat and a        *
 '*  n x k matrix rmat of right hand sides. Here mat must be           *
 '*  nonsingular.                                                      *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   Input parameters:                                                *
 '*   ================                                                 *
 '*      n        integer;  (n > 0)                                    *
 '*               Dimension of xmat.                                   *
 '*      k        integer k;  (k > 0)                                  *
 '*               number of right hand sides                           *
 '*      xmat     matrix(0:n-1,0:n-1) stored in a vector(n*n)          *
 '*               n x n original system matrix                         *
 '*      rmat     matrix(0:n-1,0:n-1) stored in a vector(n*n)          *
 '*               Right hand sides                                     *
 '*                                                                    *
 '*   Output parameter:                                                *
 '*   ================                                                 *
 '*      rmat     solution matrix for the system.                      *
 '*               The input right hand sides are lost.                 *
 '*                                                                    *
 '*   Return value (code):                                             *
 '*   ===================                                              *
 '*      = 0      all ok                                               *
 '*      = 1      n < 1 or k < 1 or invalid input parameter            *
 '*      = 2      lack of memory (not used here)                       *
 '*      = 3      mat is numerically singular.                         *
 '*                                                                    *
 '* ------------------------------------------------------------------ *
 '*                                                                    *
 '*   subroutine used:                                                 *
 '*   ================                                                 *
 '*                                                                    *
 '*      gaudec:  LU decomposition of mat.                             *
 '*                                                                    *
 '*====================================================================*
 Dim signd As Integer, rc As Integer
 ReDim perm(n)
 ReDim lu(n * n), x(n)
 Dim MACH_EPS As Double
 
  If n < 1 Or k < 1 Then                        'Invalid parameter
    code = 1
    Exit Sub
  End If

  MACH_EPS = 1.2E-16

  Call gaudec(n, xmat, lu, perm, signd, rc)     'compute factorization
                                                'in matrix lu
  If rc <> 0 Or signd = 0 Then                  'if not possible
    code = 3                                    'exit with code=3
    Exit Sub
  End If

  For m = 0 To k - 1                   'Loop over the right hand sides
    For i = 0 To n - 1                               'Updating the b's
      x(i) = rmat(n * perm(i) + m)
      For j = 0 To i - 1
        x(i) = x(i) - lu(n * i + j) * x(j)
      Next j
    Next i

    For i = n - 1 To 0 Step -1                      'back substitution
      Sum = 0#
      For j = i + 1 To n - 1
        Sum = Sum + lu(n * i + j) * x(j)
      Next j

      If Abs(xlu(n * i + i)) < MACH_EPS Then 'invalid LU decomposition
        code = 2
        Exit Sub
      End If
      x(i) = (x(i) - Sum) / lu(n * i + i)
    Next i

    For j = 0 To n - 1                                    'Save result
      rmat(n * j + m) = x(j)
    Next j
  Next m

  code = 0

End Sub

' ------------------------ END fgauss.bas --------------------------