I have been playing around with creating forward substitution and backward substitution J scripts to solve lower and upper triangular matrices of the type you would get by performing LU decomposition on a single square matrix. I could certainly just use matrix divide but this starts to slow down as the size of the matrix NxN gets bigger. I have come up with a partial J like solution but the script requires a for loop. I was wondering if someone has a mor J-ish solution? Here is what I have done so far.
BACKGROUND and SCRIPTS: >From the LU decomp paper on the J website: >http://www.jsoftware.com/jwiki/Essays/LU%20Decomposition There is a routine to perform the LU decomposition of a matrix. Rather than use this though I cheated and borrowed a worked example from the web: The original matrix before LU ]A1 =: 3 3$ 1 2 3 2 _4 6 3 _9 _3 1 2 3 2 _4 6 3 _9 _3 b1 =: 5 _1 2 ]L1 =: 3 3$1 0 0 2 _8 0 3 _15 _12 1 0 0 2 _8 0 3 _15 _12 ]U1 =: 3 3$3 4 5 0 2 8 0 0 10 3 4 5 0 2 8 0 0 10 The forward substitution and backward substitution routines are as follows: NB. Forward substitution - multiply upper triangular matrix by vector NB. Usage: matrix forsub vector forsub =: 4 : 0 a =. x b =. y z =. (#b)$z =. 0 0 0 z =. ((0 {b)% (<0 0) {a) 0} z for_i. 1 + i. (#b)-1 do. z1 =.((<i,i ){a)%~ ( (i{b) - (+/ (i{. (i {a))*(i{. z) ) ) z =. z1 i}z end. z ) NB. Back substitution - multiply lower triangular matrix by vector NB. Usage: matrix baksub vector baksub =: 4 : 0 B =. x Y =. y N =. _1 + #Y z =. (N+1)$0 z =. ((N {Y)% (<N,N) {B) N} z for_i. |. i. N do. z1 =. ((<i,i){B)%~ ((i{Y) - (+/ ((i+1)}. (i {B))*((i+1)}. z))) z =. z1 i} z end. z ) Solving linear equations with just matrix divide. b1 %. A1 1 _1 2 To solve this set of linear equations with the L and U decomps you do the following: ]y1 =: L1 forsub b1 5 _1 2 ]x1 =: U1 baksub y1 1 _1 2 which agrees with the matrix divide version These matrices were taken from the internet paper found here: https://www.utdallas.edu/dept/abp/PDF_Files/LinearAlgebra_Folder/LU_Decomposition.pdf ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
