Can you double check whether you have posted what you intended? I get b1 %. A1 4.15625 1.375 _0.635417 L1 forsub b1 5 1.375 _0.635417 U1 baksub L1 forsub b1 0.517014 0.941667 _0.0635417
which does not correspond to the result you displayed. Thanks, -- Raul On Mon, Aug 31, 2015 at 9:54 PM, Thomas McGuire <tmcguir...@gmail.com> wrote: > 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
