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 --------------------------