\ Linear equation solver using Gaussian elimination with row pivoting \ as described in Noble -- "Scientific Forth: A Modern Language for \ Scientific Computing" (Mechum Banks Publishing, Ivy, VA, 1992). \ \ version of October 4th, 1999 \ --------------------------------------------------- \ (c) Copyright 1999 Julian V. Noble. \ \ Permission is granted by the author to \ \ use this software for any application pro- \ \ vided this copyright notice is preserved. \ \ --------------------------------------------------- \ This is an ANS Forth program requiring the \ FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets. \ \ Environmental dependences: \ Assumes independent floating point stack \ Uses a FORmula TRANslator for clarity FALSE [IF] Algorithm: : }}solve N 0 DO find pivot element among rows with J > I SWAP rows J & I divide row I by pivot element: N I DO A(J,I) = A(J,I) / A(J,J) N I+1 DO A(K,J) = A(K,J) - A(K,I)*A(I,J) LOOP LOOP LOOP back-solve ; [THEN] MARKER -solve \ conditional compilation INCLUDE ftran201.f [undefined] frame| [IF] INCLUDE flocals.f [THEN] [undefined] } [IF] INCLUDE arrays.f [THEN] [undefined] zdup [IF] : zdup FOVER FOVER ; [THEN] [undefined] 1/f [IF] : 1/f 1.e0 FSWAP F/ ; [THEN] \ data structures 0 VALUE Nmax FVARIABLE Det \ determinant 1.0e-20 FCONSTANT tiny \ condition criterion 1000 long 1 CELLS 1array Iperm{ \ array of permuted row labels \ locate next pivot row : LessThan F< ; : MoreThan F> ; : }}get_pivot ( A{{ col# -- Ipiv) ( f: --) 0 \ dummy argument LOCALS| Iperm Col mat{{ | \ local names Iperm{ Col } @ TO Iperm f" ABS( mat{{ Iperm Col }} ) " ( f: -- |a[col,col]| ) Col \ 1st pivot value on stack Nmax Col 1+ \ loop limits ?DO Iperm{ I } @ TO Iperm \ begin loop f" LessThan( zdup( ABS( mat{{ Iperm Col }} ) ) )" \ f" ABS( mat{{ Iperm Col }} ) " ( f: -- |a| |a'| ) \ zdup ( f: -- |a| |a'| |a| |a'|) \ F< \ new.elt > old.elt ? IF FSWAP DROP I THEN FDROP LOOP \ end loop FDROP ; \ multiply a row by a constant : }}row*x ( M{{ row --) ( f: x -- ) 0 \ dummy argument LOCALS| Iperm row# mat{{ | frame| aa | Iperm{ row# } @ TO Iperm Nmax row# ?DO \ FDUP ( f: -- x x) \ mat{{ Iperm I }} DUP F@ ( -- adr[elt]) ( f: -- x x elt) \ F* F! f" mat{{ Iperm I }} = mat{{ Iperm I }} * aa " LOOP ( f" aa" ) |frame ; \ Usage: A{{ 2 }}row*x \ subtract a row times a constant from another row \ this is the innermost loop -- CODE for speed! : }}r1-r2*x ( M{{ r1 r2 -- ) ( f: x -- x) \ initialize assumed 0 0 LOCALS| I1 I2 r2 r1 mat{{ | \ local names frame| aa | \ local fvariable Iperm{ r1 } @ TO I1 Iperm{ r2 } @ TO I2 Nmax r2 ?DO \ begin loop f" mat{{ I1 I }} = mat{{ I1 I }} - mat{{ I2 I }} * aa" LOOP \ end loop ( f" aa") |frame ( f: -- x) ; \ v[i] = v[i] - v[j] * x : }v1-v2*x ( V{ r1 r2 -- f: x -- ) 0 0 LOCALS| Ir1 Ir2 r2 r1 v{ | frame| aa | Iperm{ r1 } @ TO Ir1 Iperm{ r2 } @ TO Ir2 f" v{ Ir1 } = v{ Ir1 } - v{ Ir2 } * aa" |frame ; \ permute row labels : mem_swap ( adr1 adr2 --) LOCALS| a2 a1 | a1 @ a2 @ a1 ! a2 ! ; : }swap ( I{ m n -- ) \ exchange 2 elts in an integer array LOCALS| N M I{ | I{ M } I{ N } mem_swap ; : initialize ( A{{ V{ -- A{{ V{ ) DUP 2@ DROP TO Nmax \ Nmax = # of equations Nmax 0 DO I Iperm{ I } ! LOOP \ init loop-label array f" Det = 1" \ init determinant ; : update_Det ( -- ) ( f: x -- x ) frame| aa | f" FABS( FDUP( - Det * aa ) )" ( f: -- D' = -x * D |D|) tiny F< ABORT" Ill-conditioned matrix" f" Det =" f" aa" |frame ; : triangularize ( A{{ V{ -- A{{ V{ ) \ assume INITIALIZEd 0 0 \ dummy arguments LOCALS| row Ipiv vec{ mat{{ | \ local names 0e0 frame| aa | Nmax 0 DO \ outer loop - by rows mat{{ I }}get_pivot \ find pivot in col I TO Ipiv \ pivot index Iperm{ I Ipiv }swap \ exchange rows Iperm{ I } @ TO row \ get current row# f" mat{{ row I }}" \ pivot elt -> fstack update_Det ( f: x -- x ) f" aa=" f" aa = fdup(1/aa)" ( f: -- 1/x) mat{{ I }}row*x \ row[i] = row[i] / pivot f" vec{ row } = vec{ row } * aa" \ vec{ row } DUP F@ F* F! \ V[i] = V[i] / pivot Nmax I 1+ ?DO \ middle loop - by rows Iperm{ I } @ TO row f" fdup(mat{{ row J }})" \ multiplier -> fstack mat{{ I J }}r1-r2*x \ row[i] = row[i]-row[j]*x vec{ I J }v1-v2*x \ same for V{ and drop x LOOP \ end middle loop LOOP \ end outer loop mat{{ vec{ ; \ push these addresses : back_solve ( A{{ V{ -- V{ ) \ assume INITIALIZEd 0 0e0 \ dummy arguments LOCALS| Jperm vec{ mat{{ | FRAME| aa | \ aa = sum 0 Nmax 2 - DO \ begin outer loop f" aa = 0" ( f: sum=0) Iperm{ I } @ TO Jperm \ permuted row index Nmax I 1+ DO \ begin inner loop f" aa = aa - mat{{ Jperm I }} * vec{ Iperm{ I } @ }" LOOP \ end inner loop f" vec{ Jperm } = vec{ Jperm } + aa" -1 +LOOP |FRAME vec{ ; \ solve linear equations Ax = V \ solution vector in V{ , matrix A{{ is overwritten : report ( V{ --) \ print solution LOCALS| vec{ | Nmax 0 DO CR ." x(" I . ." )= " f" vec{ Iperm{ I } @ } " F. LOOP ; : }}solve ( A{{ V{ --) initialize triangularize back_solve report ; FALSE [IF] Say: a{{ v{ }}solve ok Test case in file mat_ex4.f [THEN]