\ -------------------------------------------------------------------- \ A matrix toolkit, containing: \ gaussj solve mat^-1 get-column get-row transpose mat* mat- mat+ mat+! \ fillmat }}absmat \ \ Forth Scientific Library Algorithm #48 \ gaussj ( 'A 'B r c -- bad? ) \ Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) \ of Numerical Recipes p 34. The input matrix A[] has r x r elements. \ B[] is an r x c input matrix containing the c right-hand side vectors. \ On output, A is replaced by its matrix inverse and B is replaced by the \ corresponding set of solution vectors. The flag is FALSE when gaussj was \ successful (it can fail because of too small pivots or memory problems). \ Note that the FSL's LU-method uses N^3/3 operation steps, which is much \ better than the Gauss-Jordan approach (N^3 steps). However, the Gauss-Jordan \ method computes the inverse matrix automatically. When the LU-method is used \ to do this it also needs N^3 steps. \ The Gauss-Jordan method should be most efficient when you have a (large) \ number of right-hand vectors (m). In this case you only need one call to get \ all solutions at once, versus m calls using lubksb. \ The Gauss-Jordan method is more convenient when iterative improvement of the \ solution is needed (see the SOLVE procedure). \ mat^-1 ( 'A n -- bad? ) \ Matrix inversion by Gauss-Jordan elimination. \ The input matrix A[] has n by n elements. On output, A is replaced by its \ matrix inverse. \ The flag is FALSE when mat^-1 was successful (it can fail because of too \ small pivots, a singular matrix, or memory problems). \ mat* ( 'A ra ca 'B rb cb xt-C -- ) \ Matrix multiplication, works for any set of (floating-point) matrices. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrices are not altered in any way. Bounds checks are done. \ get-column ( 'A ra ca xt-C c -- 'C ra 1 ) \ Cut column c out of a ra x ca matrix A and return the result as the ra x 1 \ matrix C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix A is not altered in any way. \ get-row ( 'A ra ca xt-C r -- 'C 1 ca ) \ Cut row r out of the ra x ca matrix A and return the result as the 1 x ca \ matrix C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix A is not altered in any way. \ transpose ( 'A ra ca xt-C -- ) \ Transpose a (floating-point) matrix A , that is, interchange its rows and \ columns. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix A is not altered in any way. \ mat+! ( 'A ra ca 'B rb cb -- ) \ Matrix addition of A to B. Bounds checks are done. \ mat- ( 'A ra ca 'B rb cb xt-C -- ) \ Matrix subtraction of B from A , with the result left in C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix is not altered in any way. Bounds checks are done. \ mat+ ( 'A ra ca 'B rb cb xt-C -- ) \ Matrix addition of B to A , with the result left in C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix is not altered in any way. Bounds checks are done. \ }}absmat ( 'A r c -- ) ( F: -- e ) \ Used on a 1-row or 1-column matrix A this gives the Euclidean 'length'. \ ( Gives square root of sum of squares of all elements. ) \ fillmat ( 'A r c -- ) ( F: e -- ) \ Initialize all matrix elements of A to the number e. \ solve ( 'A 'X 'Y n m MaxSteps -- steps bad? ) ( F: MaxError -- err cnv ) \ Returns the solution X to A*X=Y, where A is an n x m coefficient \ matrix, X is m x 1 , and Y is n x 1 , with m <= n. \ If m < n, solve returns an X that represents a least squares fit through \ all data points Y (n x 1). \ Again, A and Y are kept intact. \ Solve is able to solve sets of equations that are nearly singular, or \ "noisy", using a successive, automatic, refinement method. \ Refinement is done by passing in an X that is a good guess to the wanted \ solution vector. If you have no idea of the solution, use a zero-filled X. \ The maximum number of iterations is controlled with MaxSteps. \ Iteration stops when the error, measured by the norm of (A*X-Y), \ is less than MaxError; the final error is returned on the stack \ as err . The norm of the last correction to X is returned \ as cnv . The boolean bad? is false if a solution is reached, \ +1 if m > n , and +2 if a matrix inversion failed. \ This is an ANS Forth program requiring: \ 1. The Floating-Point word sets \ 2. Uses FSL words from fsl_util.xxx and (for the tests) from r250.xxx \ 3. Uses : F> FSWAP F< ; \ : F2DUP ( F: r1 r2 -- r1 r2 r1 r2 ) FOVER FOVER ; \ : 1/F ( F: r -- 1/r ) 1e FSWAP F/ ; \ : F+! ( addr -- ) ( F: r -- ) DUP F@ F+ F! ; \ : FSQR ( F: r1 -- r2 ) FDUP F* ; HERE 1+ VALUE seed ( note: '$' prefix means HEX ) : RANDOM seed $107465 * $234567 + DUP TO seed ; ( -- u ) : CHOOSE ( n -- +n ) RANDOM UM* NIP ; ( Paul Mennen, 1991 ) \ Note: the code uses 5 fp stack cells (iForth vsn 1.07) when executing \ the test words. \ See: 'Numerical recipes in Pascal, The Art of Scientific Computing', \ William H. Press, Brian P. Flannery, Saul A. Teukolsky and William \ T. Vetterling, Chapter 2 (2.1, 2.7): Solution of Linear Algebraic Equations. \ 1989; Cambridge University Press, Cambridge, ISBN 0-521-37516-9 \ (c) Copyright 1995 Marcel Hendrix. Permission is granted by the \ author to use this software for any application provided this \ copyright notice is preserved. \ CR .( GAUSSJ & MATRICES V1.0 6 May 1995 MH ) \ Vsn 1.1 Improved doc. (thanks to "C. Montgomery" ) \ Vsn 1.2 Improved doc. (thanks to Josef Gabriel ) CR .( GAUSSJ & MATRICES V1.1 19 Jan 1997 MHX ) CR .( GAUSSJ & MATRICES V1.2 11 Jul 2005 MHX ) Private: INTEGER DARRAY indxc{ \ used for bookkeeping on the pivoting INTEGER DARRAY indxr{ \ " " INTEGER DARRAY ipiv{ \ " " : search-pivot ( 'A n -- r c bad? ) -1 -1 LOCALS| irow icol n A{{ | 0e ( big ) n 0 DO \ outer loop of the search for a pivot element ipiv{ I } @ 0<> IF n 0 DO ipiv{ I } @ -1 = IF A{{ J I }} F@ FABS F2DUP F> IF FDROP ELSE FSWAP FDROP J TO irow I TO icol THEN ELSE ipiv{ I } ( singular matrix?) @ IF FDROP -1 -1 1 UNLOOP UNLOOP EXIT THEN THEN LOOP THEN LOOP ( big ) FDROP irow icol 0 ; 1e-100 FCONSTANT smallestpivot \ choose a number related to the float size Public: \ Linear equation solution by Gauss-Jordan elimination, equation (2.1.1) \ of Numerical Recipes p 34. The input matrix A[] has r x r elements. \ B[] is an r x c input matrix containing the c right-hand side vectors. \ On output, A is replaced by its matrix inverse and B is replaced by the \ corresponding set of solution vectors. The flag is FALSE when gaussj was \ successful (it can fail because of too small pivots or memory problems). : gaussj ( 'A 'B r c -- bad? ) 0 0 LOCALS| irow icol m n B{{ A{{ | & indxc{ n }malloc malloc-fail? & indxr{ n }malloc malloc-fail? OR & ipiv{ n }malloc malloc-fail? OR IF TRUE EXIT THEN n 0 DO -1 ipiv{ I } ! LOOP n 0 DO \ [i] main loop over columns to be reduced A{{ n search-pivot IF & ipiv{ }free & indxr{ }free & indxc{ }free 2DROP TRUE UNLOOP EXIT THEN TO icol TO irow 1 ipiv{ icol } +! \ We now have the pivot element, so we interchange rows, if needed, to \ put the pivot element on the diagonal. The columns are not physically \ interchanged, only relabeled: indexc^[i], the column of the ith pivot \ element, is the ith column that is reduced, while indexr^[i] is the \ row in which that pivot element was originally located. If indexr^[i] \ <> indexc^[i] there is an implied column interchange. With this form \ of bookkeeping, the solution b's will end up in the correct order, and \ the inverse matrix will be scrambled by columns. irow icol <> IF n 0 DO A{{ irow I }} DUP F@ A{{ icol I }} DUP F@ SWAP F! F! LOOP m 0 ?DO B{{ irow I }} DUP F@ B{{ icol I }} DUP F@ SWAP F! F! LOOP THEN irow indxr{ I } ! icol indxc{ I } ! A{{ icol DUP }} DUP F@ 1e F! FDUP FABS smallestpivot F< IF & ipiv{ }free & indxr{ }free & indxc{ }free FDROP TRUE UNLOOP EXIT THEN 1/F ( F: pivinv) n 0 DO FDUP A{{ icol I }} DUP F@ F* F! LOOP m 0 ?DO FDUP B{{ icol I }} DUP F@ F* F! LOOP FDROP n 0 DO I icol <> IF A{{ I icol }} DUP F@ ( F: dum) 0e F! n 0 DO A{{ J I }} DUP F@ FOVER A{{ icol I }} F@ F* F- F! LOOP m 0 ?DO B{{ J I }} DUP F@ FOVER B{{ icol I }} F@ F* F- F! LOOP FDROP THEN LOOP LOOP ( end main loop over the columns to be reduced ) \ Unscramble the solution in view of column interchanges 0 n 1- DO indxr{ I } @ indxc{ I } @ <> IF n 0 DO A{{ I indxr{ J } @ }} DUP F@ A{{ I indxc{ J } @ }} DUP F@ SWAP F! F! LOOP THEN -1 +LOOP & ipiv{ }free & indxr{ }free & indxc{ }free FALSE ; \ Matrix inversion by Gauss-Jordan elimination. \ The input matrix A[] has n x n elements. On output, a is replaced by its \ matrix inverse. \ The flag is FALSE when mat^-1 was successful (it can fail because of too \ small pivots, a singular matrix, or memory problems). \ This is very simple because gaussj supports c=0. : mat^-1 ( 'A n -- bad? ) HERE ( dummy 'B ) SWAP 0 gaussj ; \ Matrix multiplication, works for any set of (floating-point) matrices. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrices are not altered in any way. Bounds checks are done. : mat* ( 'A ra ca 'B rb cb xt-C -- ) 0 0 0 LOCALS| kk c3 r3 mult{{ c2 r2 b{{ c1 r1 a{{ | c1 r2 <> ABORT" mat* :: bounds mismatch" r1 TO r3 c2 TO c3 mult{{ r3 c3 }}malloc malloc-fail? ABORT" mat* :: out of memory" mult{{ EXECUTE TO mult{{ r3 0 DO I TO kk c3 0 DO 0e c1 0 ?DO a{{ kk I }} F@ b{{ I J }} F@ F* F+ ( accumulate on fstack) LOOP mult{{ J I }} F! LOOP LOOP ; \ Cut column c out of a ra x ca matrix and return the result as the ra x 1 \ matrix C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix A is not altered in any way. : get-column ( 'A ra ca xt-C c -- 'C ra 1 ) LOCALS| col b{{ c1 r1 a{{ | b{{ r1 1 }}malloc malloc-fail? ABORT" get-column :: out of memory" b{{ EXECUTE TO b{{ r1 0 ?DO a{{ I col }} F@ b{{ I 0 }} F! LOOP b{{ r1 1 ; \ Cut row r out of the ra x ca matrix A and return the result as a the 1 x ca \ matrix C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix A is not altered in any way. : get-row ( 'A ra ca xt-C r -- 'C 1 ca ) LOCALS| row b{{ c1 r1 a{{ | b{{ 1 c1 }}malloc malloc-fail? ABORT" get-row :: out of memory" b{{ EXECUTE TO b{{ c1 0 ?DO a{{ row I }} F@ b{{ 0 I }} F! LOOP b{{ 1 c1 ; \ Transpose a (floating-point) matrix A , that is, interchange its rows and \ columns. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix A is not altered in any way. : transpose ( 'A ra ca xt-C -- ) LOCALS| b{{ c1 r1 a{{ | b{{ c1 r1 }}malloc malloc-fail? ABORT" transpose :: out of memory" b{{ EXECUTE TO b{{ r1 0 ?DO c1 0 ?DO a{{ J I }} F@ b{{ I J }} F! LOOP LOOP ; \ Matrix addition of A to B. Bounds checks are done. : mat+! ( 'A ra ca 'B rb cb -- ) LOCALS| c2 r2 b{{ c1 r1 a{{ | r1 r2 <> c1 c2 <> OR ABORT" mat+! :: bounds mismatch" r2 0 ?DO c2 0 ?DO a{{ J I }} F@ b{{ J I }} F+! LOOP LOOP ; \ Matrix subtraction of B from A , with the result left in C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix is not altered in any way. Bounds checks are done. : mat- ( 'A ra ca 'B rb cb xt-C -- ) LOCALS| c{{ c2 r2 b{{ c1 r1 a{{ | r1 r2 <> c1 c2 <> OR ABORT" mat- :: bounds mismatch" c{{ r1 c1 }}malloc malloc-fail? ABORT" mat- :: out of memory" c{{ EXECUTE TO c{{ r1 0 ?DO c1 0 ?DO a{{ J I }} F@ b{{ J I }} F@ F- c{{ J I }} F! LOOP LOOP ; \ Matrix addition of B to A , with the result left in C. \ The result matrix C (execution token on the stack) must be freed first. \ The original matrix is not altered in any way. Bounds checks are done. : mat+ ( 'A r1 c1 'B r2 c2 'head -- ) LOCALS| c{{ c2 r2 b{{ c1 r1 a{{ | r1 r2 <> c1 c2 <> OR ABORT" mat+ :: bounds mismatch" c{{ r1 c1 }}malloc malloc-fail? ABORT" mat+ :: out of memory" c{{ EXECUTE TO c{{ r1 0 ?DO c1 0 ?DO a{{ J I }} F@ b{{ J I }} F@ F+ c{{ J I }} F! LOOP LOOP ; \ Used on a 1-row or 1-column matrix A this gives the Euclidean 'length'. \ ( Gives square root of sum of squares of all elements. ) : }}absmat ( 'A r c -- ) ( F: -- e ) LOCALS| c1 r1 a{{ | 0e r1 0 ?DO c1 0 ?DO a{{ J I }} F@ FSQR F+ LOOP LOOP FSQRT ; \ Initialize all matrix elements of A to the number e. : fillmat ( 'A r c -- ) ( F: e -- ) LOCALS| c1 r1 a{{ | r1 0 DO c1 0 DO FDUP a{{ J I }} F! LOOP LOOP FDROP ; \ solve ( 'A 'X 'Y n m MaxSteps -- steps bad? ) ( F: MaxError -- err cnv ) \ Returns the solution X to A*X=Y, where A is an n x m coefficient \ matrix, X is m x 1 , and Y is n x 1 , with m <= n. \ If m < n, solve returns an X that represents a least squares fit through \ all data points Y (n x 1). Again, A and Y are kept intact. \ Solve is able to solve sets of equations that are nearly singular, using a \ successive, automatic, refinement method. Refinement is done by passing in \ an X that is a good guess to the wanted solution vector. If you have no \ idea of the solution, use a zero-filled X. \ The maximum number of iterations is controlled with MaxSteps. \ Iteration stops when the error, measured by the norm of (A*X-Y), \ is less than MaxError; the final error is returned on the stack \ as err . The norm of the last correction to X is returned \ as cnv . The boolean bad? is false if a solution is reached, \ +1 if m > n , and +2 if a matrix inversion failed. \ After a suggestion of Dr. Jos Bergervoet, personal communication. \ Solve is more powerful than mprove (Press et al). Private: FLOAT DMATRIX At{{ FLOAT DMATRIX qd{{ FLOAT DMATRIX Q{{ FLOAT DMATRIX QAt{{ FLOAT DMATRIX Ax{{ FLOAT DMATRIX dif{{ FLOAT DMATRIX delta{{ FVARIABLE maxerror Public: \ Do not forget to zero-fill X when you have no idea of the solution at all! : solve ( 'A 'X 'Y n m MaxSteps -- steps bad? ) ( F: maxerror -- err cnv ) 1 LOCALS| iters MaxSteps m n y{{ x{{ A{{ | maxerror F! n m < IF -1 1 1e38 1e38 EXIT THEN A{{ n m & At{{ transpose At{{ m n A{{ n m & Q{{ mat* \ Q is (m x m) Q{{ m mat^-1 \ Q <- (At*A)^-1 IF & Q{{ }}free & At{{ }}free \ mat^-1 failed. -1 2 1e38 1e38 EXIT THEN Q{{ m m At{{ m n & QAt{{ mat* \ QAt <- (A*At)^-1 * At & Q{{ }}free & At{{ }}free BEGIN A{{ n m x{{ m 1 & Ax{{ mat* y{{ n 1 Ax{{ n 1 & dif{{ mat- QAt{{ m n dif{{ n 1 & delta{{ mat* delta{{ m 1 x{{ m 1 mat+! dif{{ n 1 }}absmat maxerror F@ F> iters MaxSteps < AND WHILE & delta{{ }}free & dif{{ }}free & Ax{{ }}free iters 1+ TO iters REPEAT iters ( steps taken) dif{{ n 1 }}absmat ( error) delta{{ m 1 }}absmat ( convergence) & QAt{{ }}free & Ax{{ }}free & dif{{ }}free & delta{{ }}free 0 ; Reset_Search_Order TEST-CODE? [IF] \ --------------------------------------------------------- \ Read ahead in a text file. This doesn't work with a terminal. \ A nice feature: the read text is interpreted, so { 1e 2e F+ } is valid. \ Data starts on the _next_ line. : READ-INFILE REFILL 0= ABORT" REFILL :: not possible" SOURCE EVALUATE ; \ Read a matrix (won't work from the terminal). The matrix head passed \ should be empty (free the contents first). The reading starts on the next \ line of the text file. \ Example: FLOAT MATRIX A{{ & A{{ 55 20 }}FREAD \ ... \ & A{{ }}free \ & A{{ 5 2 }}FREAD \ .... : }}FREAD ( 'head rows cols -- ) LOCALS| cols rows m{{ | m{{ rows cols }}malloc malloc-fail? ABORT" }}FREAD failed" m{{ EXECUTE TO m{{ rows 0 ?DO READ-INFILE ( coefficients) 0 cols 1- DO m{{ J I }} F! -1 +LOOP LOOP REFILL 0= ABORT" REFILL :: not possible" ; \ Let's use it FLOAT DMATRIX A{{ & A{{ 3 3 }}FREAD This field is not read... 1e 8e -7e 2e -3e 4e 3e 7e 1e FLOAT DMATRIX B{{ & B{{ 3 2 }}FREAD 0e 12e ( first row's solution -> 5,2,3) 16e 4e ( second row's solution -> 3,2,1) 32e 24e FLOAT DMATRIX a-row{{ FLOAT DMATRIX a-column{{ FLOAT DMATRIX differ{{ CR .( TEST-MAT shows get-row get-column transpose mat^-1 mat- }}absmat mat+!) : TEST-MAT CR ." (GET-ROW : Should be [1,8,-7])" A{{ 3 3 & a-row{{ 0 get-row ROT CR }}fprint & a-row{{ }}free CR ." (GET-COLUMN : Should be [8,-3,7]^T)" A{{ 3 3 & a-column{{ 1 get-column ROT CR }}fprint CR ." (TRANSPOSE : Should be [8,-3,7])" a-column{{ 3 1 & a-row{{ transpose 1 3 a-row{{ CR }}fprint & a-row{{ }}free & a-column{{ }}free & a-row{{ 3 3 }}malloc 3 0 DO 3 0 DO I J + CHOOSE 1+ S>F a-row{{ J I }} F! LOOP LOOP 3 3 a-row{{ CR ." A matrix " CR }}fprint a-row{{ 3 mat^-1 IF CR ." Inverse failed (*NOT* an error)!" ELSE 3 3 a-row{{ CR ." Its inverse " CR }}fprint THEN & a-row{{ }}free CR ." Matrix A - Matrix B --> Matrix C " A{{ 3 3 & a-row{{ 0 get-row ROT CR ." A: " }}fprint A{{ 3 3 & a-column{{ 1 get-row ROT ." B: " }}fprint a-row{{ 1 3 a-column{{ 1 3 & differ{{ mat- 1 3 differ{{ ." C: " }}fprint & a-row{{ }}free & a-column{{ }}free & differ{{ }}free CR ." Matrix A + Matrix B --> Matrix C " A{{ 3 3 & a-row{{ 1 get-row ROT CR ." A: " }}fprint A{{ 3 3 & a-column{{ 2 get-row ROT ." B: " }}fprint a-row{{ 1 3 a-column{{ 1 3 & differ{{ mat+ 1 3 differ{{ ." C: " }}fprint CR ." The length of vector C = " differ{{ 1 3 }}absmat F. CR ." C + [33,44,55] = " 33e a-row{{ 0 0 }} F! 44e a-row{{ 0 1 }} F! 55e a-row{{ 0 2 }} F! a-row{{ 1 3 differ{{ 1 3 mat+! 1 3 differ{{ }}fprint & differ{{ }}free & a-column{{ }}free & a-row{{ }}free ; CR .( 3EQS solves a 3 x 3 set of linear equations with 2 right-hand vectors.) FLOAT DMATRIX C{{ FLOAT DMATRIX oldA{{ FLOAT DMATRIX oldB{{ : 3EQS ( -- ) & oldA{{ 3 3 }}malloc A{{ oldA{{ 3 3 }}fcopy & oldB{{ 3 2 }}malloc B{{ oldB{{ 3 2 }}fcopy CR ." Original A{{ " CR 3 3 A{{ }}fprint CR ." Original B{{ " CR 3 2 B{{ }}fprint A{{ B{{ 3 2 gaussj ABORT" gaussj failed" CR ." A^-1 " CR 3 3 A{{ }}fprint CR ." Solution vectors " CR 3 2 B{{ }}fprint oldA{{ 3 3 A{{ 3 3 & C{{ mat* CR ." A^-1 x A " CR 3 3 C{{ }}fprint & C{{ }}free CR ." A^-1 x b[..,0] " oldB{{ 3 2 & a-column{{ 0 get-column ( 'a r c -- ) 2>R >R A{{ 3 3 R> 2R> & C{{ mat* CR 3 1 C{{ }}fprint & C{{ }}free & a-column{{ }}free CR ." A^-1 x b[..,1] " oldB{{ 3 2 & a-column{{ 1 get-column ( 'a r c -- ) 2>R >R A{{ 3 3 R> 2R> & C{{ mat* CR 3 1 C{{ }}fprint & C{{ }}free ( & a-column{{ }}free ) & differ{{ 3 1 }}malloc differ{{ 3 1 1e fillmat oldA{{ differ{{ a-column{{ 3 3 10 1e-17 solve ?DUP IF CR ." Solve failed with error " . DROP F2DROP ELSE CR ." Solved A for " CR 3 1 a-column{{ }}fprint ." -> " CR 3 1 differ{{ }}fprint CR ." Used " . ." iterations, convergence = " F. ." error = " F. THEN & differ{{ }}free & a-column{{ }}free oldB{{ B{{ 3 2 }}fcopy & oldB{{ }}free oldA{{ A{{ 3 3 }}fcopy & oldA{{ }}free ; CR .( SOLVE-IT finds a LSQ approximation through the data points in Y, given) CR .( a function described by the unknown coefficients in X. There are more) CR .( data points than unknowns, and noise is present.) CR .( Note that only 2 steps provide an adequate result already!) : SOLVE-IT ( -- ) & oldA{{ 4 3 }}malloc oldA{{ 4 3 0e fillmat 1e oldA{{ 0 0 }} F! \ 1 * x1 = 4 1e oldA{{ 1 1 }} F! \ 1 * x2 = 5 1e oldA{{ 2 2 }} F! \ 1 * x3 = 6 1e oldA{{ 3 0 }} F! \ an extra row (sum of above 3!) 1e oldA{{ 3 1 }} F! 1e oldA{{ 3 2 }} F! 1e-16 oldA{{ 0 2 }} F+! \ Add some noise to coefficients -1e-15 oldA{{ 1 0 }} F+! 1e-17 oldA{{ 2 1 }} F+! -1e-18 oldA{{ 3 1 }} F+! & differ{{ 3 1 }}malloc differ{{ 3 1 0e fillmat & oldB{{ 4 1 }}malloc 4e oldB{{ 0 0 }} F! 5e oldB{{ 1 0 }} F! 6e oldB{{ 2 0 }} F! 15e oldB{{ 3 0 }} F! 1e-16 oldB{{ 0 0 }} F+! \ Add some noise to Y vector -1e-15 oldB{{ 1 0 }} F+! 5e-17 oldB{{ 2 0 }} F+! -3e-16 oldB{{ 3 0 }} F+! CR ." The result we're looking for is [4,5,6]^T" CR oldA{{ differ{{ oldB{{ 4 3 10 1e-9 solve ?DUP IF CR ." Solve failed with error " . DROP F2DROP ELSE CR ." Solve (err < 1e-9) finished after " . ." steps, result = " 3 1 differ{{ CR }}fprint CR ." RMS difference with y vector = " FSWAP FS. 3 SPACES CR ." Last increment used for x = " FS. THEN CR oldA{{ differ{{ oldB{{ 4 3 10 1e-16 solve ?DUP IF CR ." Refinement step failed with error " . DROP F2DROP ELSE CR ." Refinement (err < 1e-16) finished after " . ." steps, result = " 3 1 differ{{ CR }}fprint CR ." RMS difference with y vector = " FSWAP FS. 3 SPACES CR ." Last increment used for x = " FS. THEN & oldB{{ }}free & differ{{ }}free & oldA{{ }}free ; [THEN] ( * End of File * )