{**************************************************************** * Gauss algorithm for solving linear equations * * ------------------------------------------------------------- * * REF.: "Numerical Algorithms with C, By Gisela Engeln-Muellges * * and Frank Uhlig, Springer-Verlag, 1996" [BIBLI 11]. * * * * TPW Release By J-P Moreau, Paris. * * (www.jpmoreau.fr) * ****************************************************************} UNIT FGAUSS; INTERFACE Const MACH_EPS = 1.2e-16; MAXROOT = 1e16; NMAX = 10; ONE = 1.0; ZERO = 0.0; Type pMAT = ^MATRIX; MATRIX = Array[0..NMAX,0..NMAX] of Double; pVEC = ^VECTOR; VECTOR = Array[0..NMAX] of Double; pIVEC = ^IVEC; IVEC = Array[0..NMAX] of Integer; Procedure gauss(mode,n: integer; mat: pMAT; lumat: pMAT; perm: pIVEC; b:pVEC; x: pVEC; var signd: integer; var code: integer); Procedure mgauss {Gauss for multiple right hand sides} ( n, {Dimension of system .............} k: Integer; {number of right hand sides ......} mat: pMAT; {original matrix .................} rmat: pMAT; {Right hand sides/solutions ......} var code: Integer {Error code ......................} ); IMPLEMENTATION {headers of functions used by gauss} Procedure gaudec {Gauss decomposition .......................} ( n: integer; { size of matrix ..................} mat: pMAT; { Input matrix ....................} lumat: pMAT; { matrix decomposition ............} perm: pIVEC; { row interchanges ................} var signd: integer; { sign of perm ....................} var rc: integer { error code ......................} ); Forward; Procedure gausol { Gauss solution ............................} ( n: integer; { size of matrix ..................} lumat: pMAT; { decomposed matrix (LU) ..........} perm: pIVEC; { row permutation vector ..........} b: pVEC; { Right hand side .................} x: pVEC; { solution ........................} var rc: integer { error code ......................} ); Forward; Procedure ISWAP(var a,b: integer); Var tmp:integer; Begin tmp:=b; b:=a; a:=tmp End; Procedure SWAP(var a,b: double); Var tmp:double; Begin tmp:=b; b:=a; a:=tmp End; Procedure gauss {Gauss algorithm for solving linear equations} ( mode, { Modus: 0, 1, 2, 3 ...............} n: integer; { Dimension of matrix .............} mat: pMAT; { Input matrix ....................} lumat: pMAT; { LU decomposition ................} perm: pIVEC; { row remutation vector ...........} b: pVEC; { right hand side .................} x: pVEC; { solution of the system ..........} var signd: integer; { sign of the permutation .........} var code: integer ); {====================================================================* * * * 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. * * mat type pMAT; pointer to matrix of the linear system. * * lumat type pMAT; (for mode = 2, 3) * * LU factors of mat * * lumat can be stored in the space of mat. * * perm type pIVEC; (for mode = 2, 3) * * Permutation vector, of the row interchangfes in lumat* * b type pVEC; (for mode = 0, 2, 3) * * Right hand side of the system. * * signd type Integer; (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: * * ================== * * lumat type pMAT; (for mode = 0, 1) * * LU factorization of mat. * * perm type pIVEC; (for mode = 0, 1) * * row ermutation vektor * * x type pVEC; (for mode = 0, 2, 3) * * pointer to solution vector. * * signd type Integer; (for mode = 0, 1) * * sign of perm. * * * * Return value (code): * * =================== * * =-1 Max. number (MAXITER) of iterative refinements * * reached (MAXITER) while mod = 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 * * * *====================================================================* * * * Procedures used: * * ================ * * * * gaudec: determines LU decomposition * * gausol: solves the linear system * * * *====================================================================} Label fin; {to exit procedure} Var rc: integer; Begin if n < 1 then begin code := 1; goto fin end; Case mode of 0: begin {Find factorization and solve system ...................} gaudec (n, mat, lumat, perm, signd,rc); if rc = 0 then begin gausol(n, lumat, perm, b, x,rc); code:=rc end else code:=rc; goto fin end; 1: begin {Find factorization only ...............................} gaudec(n, mat, lumat, perm, signd,rc); code:=rc; goto fin end; 2: begin {Solve only ............................................} gausol(n, lumat, perm, b, x,rc); code:=rc; goto fin end; 3: begin {solve and then use iterative refinement ...............} writeln(' fgauss: gausoli not implemented.'); code:=5; goto fin end End; code := 5; {Wrong call} fin:End; Procedure gaudec {Gauss decomposition .......................} ( n: integer; { size of matrix ..................} mat: pMAT; { Input matrix ....................} lumat: pMAT; { matrix decomposition ............} perm: pIVEC; { row interchanges ................} var signd: integer; { sign of perm ....................} var rc: integer { error code ......................} ); {====================================================================* * * * 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. * * mat type pMAT; * * pointer to original system matrix in vector form. * * * * Output parameters: * * ================== * * lumat type pMAT; * * pointer to LU factorization * * perm type pIVEC; * * pointer to row permutation vector for lumat * * signd TYPE integer; * * sign of perm. The determinant of mat can be computed * * as the product of the diagonal entreis of lumat 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 * * * *====================================================================} Label fin; {to exit procedure} Var m, j, i, j0: Integer; piv, tmp, zmax: Double; d: VECTOR; {scaling vector for pivoting} Begin if n < 1 then begin rc:=1; {Invalid parameters} goto fin end; {copy mat to lumat} for i:=0 to n-1 do for j:=0 to n-1 do lumat^[i,j] := mat^[i,j]; for i := 0 to n-1 do begin perm^[i] := i; { Initialize perm} zmax := ZERO; for j := 0 to n-1 do {find row maxima} begin tmp := ABS(lumat^[i,j]); if tmp > zmax then zmax := tmp end; if zmax = ZERO then {mat is singular} begin rc := 3; goto fin end; d[i] := ONE / zmax end; signd := 1; {initialize sign of perm} for i := 0 to n-1 do begin piv := ABS(lumat^[i,i]) * d[i]; j0 := i; {Search for pivot element} for j := i + 1 to n-1 do begin tmp := ABS(lumat^[j,i]) * d[j]; if piv < tmp then begin piv := tmp; {Mark pivot element and} j0 := j {its location} end end; if piv < MACH_EPS then {If piv is small, mat is} begin {nearly singular} signd := 0; rc :=4; goto fin end; if j0 <> i then begin signd := - signd; {update signd} ISWAP(perm^[j0], perm^[i]); {swap pivot entries} SWAP(d[j0], d[i]); {swap scaling vector} for j:=0 to n-1 do SWAP(lumat^[j0,j], lumat^[i,j]) {swap j0-th and i-th rows of lumat} end; for j := i + 1 to n-1 do {Gauss elimination} begin if lumat^[j,i] <> ZERO then begin lumat^[j,i] := lumat^[j,i] / lumat^[i,i]; tmp := lumat^[j,i]; for m := i + 1 to n-1 do lumat^[j,m] := lumat^[j,m] - tmp * lumat^[i,m] end end end; {i loop} rc:=0; {all ok} fin:End; Procedure gausol { Gauss solution ............................} ( n: integer; { size of matrix ..................} lumat: pMAT; { decomposed matrix (LU) ..........} perm: pIVEC; { row permutation vector ..........} b: pVEC; { Right hand side .................} x: pVEC; { solution ........................} var rc: integer { 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 lumat, * * lumat type pMAT; * * pointer to LU factorization, as produced from gaudec * * perm type pIVEC; * * pointer to row permutation vector for lumat * * b type pVEC; * * pointer to Right hand side of the system. * * * * Output parameter: * * ================ * * x type pVEC; * * pointer to solution vector * * * * Return value (rc): * * ================= * * = 0 all ok * * = 1 n < 1 or other invalid input parameter * * = 3 improper LU decomposition ( zero diagonal entry) * * * *====================================================================} Label fin; {to exit procedure} Var j,k: integer; sum: double; Begin if n < 1 then begin rc:=1; {Invalid input parameter} goto fin end; for k := 0 to n-1 do {update b} begin x^[k] := b^[perm^[k]]; for j := 0 to k-1 do x^[k] := x^[k] - lumat^[k,j] * x^[j] end; for k := n - 1 Downto 0 do {back substitute} begin sum := ZERO; for j := k + 1 to n-1 do sum := sum + lumat^[k,j] * x^[j]; if ABS(lumat^[k,k]) < MACH_EPS then begin rc:=3; goto fin end; x^[k] := (x^[k] - sum) / lumat^[k,k] end; rc := 0; {all ok} fin:End; Function det {Determinant ..............................} ( n: integer; {Dimension of the matrix .........} mat:pMAT {matrix ..........................} ): Double; {====================================================================* * * * det computes the determinant of an n x n real matrix mat * * * *====================================================================* * * * Input parameter: * * ================ * * n integer; (n > 0) * * Dimension of mat * * mat pointer to (type pMAT) * * n x n matrix * * * * Return value: * * ============= * * REAL Determinant of mat. * * If the return value = 0, then the matrix is singular * * or the storage is insufficient * * * *====================================================================* * * * Functions in use : * * =================== * * * * int gaudec (): LU decomposition of mat * * * *====================================================================} Label fin; {to exit function} Var i, rc, signd: Integer; perm: pIVEC; lu: pMAT; tmpdet: Double; Begin New(perm); New(lu); if n < 1 then begin det:=ZERO; goto fin end; gaudec(n, mat, lu, perm, signd,rc); {decompose} if (rc <> 0) or (signd = 0) then begin det:=ZERO; goto fin end; tmpdet := 1.0*signd; for i := 0 to n-1 do if ABS(tmpdet) < MACH_EPS then begin det:=ZERO; goto fin end else if (ABS(tmpdet) > MAXROOT) or (ABS(lu^[i,i]) > MAXROOT) then begin det:=MAXROOT; goto fin end else tmpdet := tmpdet * lu^[i,i]; {compute det} det:=tmpdet; fin: Dispose(perm); Dispose(lu) End; Procedure mgauss {Gauss for multiple right hand sides} ( n, {Dimension of system .............} k: Integer; {number of right hand sides ......} mat: pMAT; {original matrix .................} rmat: pMAT; {Right hand sides/solutions ......} var code: Integer {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 mat. * * k integer k; (k > 0) * * number of right hand sides * * mat type pMAT; pointer to a * * n x n original system matrix * * rmat type pMAT; * * pointer to matrix of right hand sides * * * * Output parameter: * * ================ * * rmat type pMAT; * * 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. * * * * ------------------------------------------------------------------ * * * * procedures used: * * ================ * * * * gaudec: LU decomposition of mat. * * * *====================================================================} Label fin; {to exit procedure} Var i, j: integer; m, signd, rc: integer; perm: pIVEC; lu: pMAT; x: pVEC; sum: Double; Begin New(perm); New(lu); New(x); if (n < 1) or (k < 1) then {Invalid parameter} begin code:=1; goto fin end; gaudec (n, mat, lu, perm, signd, rc); {compute factorization} {in matrix lu } if (rc <> 0) or (signd = 0) then {if not possible } begin code:=3; {exit with code=3 } goto fin end; for m := 0 to k-1 do {Loop over the right hand sides} begin for i := 0 to n-1 do {Updating the b's} begin x^[i] := rmat^[perm^[i],m]; for j := 0 to i-1 do x^[i] := x^[i] - lu^[i,j] * x^[j] end; for i := n - 1 Downto 0 do {back substitution} begin sum := ZERO; for j := i + 1 to n-1 do sum := sum + lu^[i,j] * x^[j]; if ABS(lu^[i,i]) < MACH_EPS then {invalid LU decomposition} begin code := 2; goto fin end; x^[i] := (x^[i] - sum) / lu^[i,i]; end; for j := 0 to n-1 do {Save result} rmat^[j,m] := x^[j] end; code := 0; fin: Dispose(perm); Dispose(lu); Dispose(x) End; END. {of unit} { ------------------------ END fgauss.pas -------------------------- }