Sorry I should have rerun the stuff rather than just pick through my debugging session.
The original Matrix ]A1 =: 3 3$1 2 3 2 _4 6 3 _9 _3 1 2 3 2 _4 6 3 _9 _3 LU decomposition triangular matrix L1 and U1 are L1 1 0 0 2 _8 0 3 _15 _12 U1 1 2 3 0 1 0 0 0 1 b1 =: 5 18 6 NB. this is the right hand side of a system of equations of the for Ax = b Matrix divide solution b1 %. A1 1 _1 2 Solution using the LU decomposition matrices ]y1 =: L1 forsub b1 5 _1 2 U1 baksub y1 1 _1 2 Tom McGuire On Aug 31, 2015, at 10:23 PM, Raul Miller <rauldmil...@gmail.com> wrote: > 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
