\ backsub Solves for linear systems via LU factorization \ Forth Scientific Library Algorithm #35 \ Solves Ax = b (A in LU form) \ Presumes that the matrix has been converted in LU form (using LUFACT) \ before being called. \ This code is an ANS Forth program requiring: \ 1. The Floating-Point word set \ 2. Uses words 'Private:', 'Public:' and 'Reset_Search_Order' \ to control the visibility of internal code. \ 3. Uses the words 'FLOAT' and Array to create floating point arrays. \ 4. The word '}' to dereference a one-dimensional array. \ 5. Uses the words 'DARRAY' and '&!' to set array pointers. \ 6. The compilation of the test code is controlled by the VALUE TEST-CODE? \ and the conditional compilation words in the Programming-Tools wordset \ see, \ Baker, L., 1989; C Tools for Scientists and Engineers, \ McGraw-Hill, New York, 324 pages, ISBN 0-07-003355-2 \ (c) Copyright 1994 Everett F. Carter. Permission is granted by the \ author to use this software for any application provided this \ copyright notice is preserved. \ ( defective test code commented out 22Jan05 cgm ) CR .( BACKSUB V1.2 18 April 1995 EFC ) \ Private: FLOAT DARRAY b{ FLOAT DMATRIX a{{ INTEGER DARRAY pivot{ FLOAT DARRAY sx{ FLOAT DARRAY sy{ : sdot ( k n first -- ) ( F: -- t ) \ based upon a standard BLAS routine 0.0E0 ?DO a{{ I 2 PICK }} F@ b{ I } F@ F* F+ LOOP DROP ; : backsub-init ( 'dlu 'b -- n ) & b{ &! 2DUP ->PIVOT{ & pivot{ &! 2DUP ->MATRIX{{ & a{{ &! ->N @ ; : solve-Ly=b ( n -- ) -1 SWAP 0 DO b{ pivot{ I } @ } DUP F@ b{ I } F@ F! DUP 0 < IF FDUP F0= 0= IF DROP I THEN ELSE I OVER ?DO a{{ J I }} F@ b{ I } F@ F* F- LOOP THEN b{ I } F! LOOP DROP ; : solve-Ux=y ( n -- ) 0 OVER 1- DO b{ I } F@ I OVER 1- < IF DUP I 1+ DO a{{ J I }} F@ b{ I } F@ F* F- LOOP THEN a{{ I I }} F@ F/ b{ I } F! -1 +LOOP DROP ; : solve-UTy=b ( n -- ) -1 SWAP 0 DO b{ I } F@ DUP 0 < IF FDUP F0= 0= IF DROP I THEN ELSE I OVER ?DO a{{ I J }} F@ b{ I } F@ F* F- LOOP THEN a{{ I I }} F@ F/ b{ I } F! LOOP DROP ; : solve-LTx=y ( n -- ) 0 OVER 1- DO b{ I } F@ I OVER 1- < IF DUP I 1+ DO a{{ I J }} F@ b{ I } F@ F* F- LOOP THEN b{ I } F! -1 +LOOP DROP ; Public: : backsub ( 'dlu 'b -- ) backsub-init DUP solve-Ly=b solve-Ux=y ; Reset_Search_Order TEST-CODE? [IF] \ test code ============================================== 3 3 FLOAT matrix mat{{ 3 3 FLOAT matrix lmat{{ 3 FLOAT array x{ 3 INTEGER array piv{ LUMATRIX lub : init-vals1 ( -- ) % 1.0 mat{{ 0 0 }} F! % 0.0 mat{{ 0 1 }} F! % 5.0 mat{{ 0 2 }} F! % 3.0 mat{{ 1 0 }} F! % 2.0 mat{{ 1 1 }} F! % 4.0 mat{{ 1 2 }} F! % 1.0 mat{{ 2 0 }} F! % 1.0 mat{{ 2 1 }} F! % 6.0 mat{{ 2 2 }} F! % 0.0 x{ 0 } F! % 4.0 x{ 1 } F! % 2.0 x{ 2 } F! ; : init-vals2 ( -- ) % 2.0 mat{{ 0 0 }} F! % 3.0 mat{{ 0 1 }} F! % -1.0 mat{{ 0 2 }} F! % 4.0 mat{{ 1 0 }} F! % 4.0 mat{{ 1 1 }} F! % -3.0 mat{{ 1 2 }} F! % -2.0 mat{{ 2 0 }} F! % 3.0 mat{{ 2 1 }} F! % -1.0 mat{{ 2 2 }} F! % 5.0 x{ 0 } F! % 3.0 x{ 1 } F! % 1.0 x{ 2 } F! ; : transpose ( &a n -- ) DUP 1- 0 DO DUP I DO OVER I J }} DUP F@ 2 PICK J I }} F@ F! OVER J I }} F! LOOP LOOP 2DROP ; : backsub-test ( -- ) CR ." --------------problem 1------------------ " CR init-vals1 lub lmat{{ piv{ 3 LUMATRIX-INIT ." A: " CR 3 3 mat{{ }}fprint CR ." B: " 3 x{ }fprint CR mat{{ lub lufact CR 3 3 lmat{{ }}fprint lub x{ backsub CR ." solution (should be 0 2 0 ): " 3 x{ }fprint CR ." --------------problem 2----------------- " CR init-vals2 lub lmat{{ piv{ 3 LUMATRIX-INIT ." A: " CR 3 3 mat{{ }}fprint CR ." B: " 3 x{ }fprint CR mat{{ lub lufact CR 3 3 lmat{{ }}fprint lub x{ backsub CR ." solution (should be 1 2 3 ): " 3 x{ }fprint CR ; \ : backsubt-test ( -- ) \ solves the same problem in transposed form \ CR ." --------------problem 1------------------ " CR \ init-vals1 \ mat{{ 3 transpose \ lub lmat{{ piv{ 3 LUMATRIX-INIT \ ." A: " CR \ 3 3 mat{{ }}fprint \ CR ." B: " 3 x{ }fprint \ CR \ mat{{ lub lufact \ CR \ 3 3 lmat{{ }}fprint \ lub x{ backsub-t \ CR ." solution (should be 0 2 0 ): " \ 3 x{ }fprint \ CR ." --------------problem 2----------------- " CR \ init-vals2 \ mat{{ 3 transpose \ lub lmat{{ piv{ 3 LUMATRIX-INIT \ ." A: " CR \ 3 3 mat{{ }}fprint \ CR ." B: " 3 x{ }fprint \ CR \ mat{{ lub lufact \ CR \ 3 3 lmat{{ }}fprint \ lub x{ backsub-t \ CR ." solution (should be 1 2 3 ): " \ 3 x{ }fprint \ CR \ ; [THEN]