# LUDecompose, LUBacksubstitue in class scalarMatrix

 Register Blogs Members List Search Today's Posts Mark Forums Read

 July 1, 2009, 05:43 LUDecompose, LUBacksubstitue in class scalarMatrix #1 Senior Member   Join Date: Mar 2009 Posts: 248 Rep Power: 16 Dear Developers Good Day to All ! While using the scalarMatrix class (/src/OpenFOAM/matrices/scalarMatrix) I tried the LUDecompose and LUBacksubstitute functions. They work fine but when i tried to compare the results with the MATLAB answers , it seems the decomposed matrix has some special format to store the L and U parts. Here is what I did: Code: ```int main(int argc, char *argv[]) { scalarMatrix hmm(3); hmm[0][0] = -3.0; hmm[0][1] = 10.0; hmm[0][2] = -4.0; hmm[1][0] = 2.0; hmm[1][1] = 3.0; hmm[1][2] = 10.0; hmm[2][0] = 2.0; hmm[2][1] = 6.0; hmm[2][2] = 1.0; Info<< hmm << endl; labelList pivotIndices(hmm.n()); hmm.LUDecompose(hmm, pivotIndices); scalarField source(hmm.n(), 0); source[0] = 1; source[1] = 3; source[2] = 4; hmm.LUBacksubstitute(hmm, pivotIndices, source); Info << hmm << endl; Info << source << endl; Info << "End\n" << endl; return 0; }``` the matrix is a 3 x 3 square matrix hmm = ((-3 10 -4) (2 3 10) (2 6 1)) According to the class documentation LUDecompose will decompose the matrix with pivoting, which in my understanding implies MATLAB command [L,U] = lu(A). I called hmm.LUDecompose(hmm, pivotIndices), and assuming that the decomposition is saved back in the scalarMatrix hmm , the Info << hmm << endl; after decomposition gives: 3 3 ((2 6 1) (-1.5 19 -2.5) (1 -0.157895 8.60526)) and Info << source << endl; gives the solution 3(0.883792 0.370031 0.0122324) Now on the other hand using MATLAB one gets: Code: ```Lower Matrix L: 1 0 0 -0.6667 1 0 -0.6667 1.31034 1 Upper Matrix U: -3 10 -4 0 9.667 7.333 0 0 -11.275 Solution : 0.8837 0.370036 0.0122324``` The solution matches perfectly in both cases but the matrix returned by OpenFOAM has special storage format which is not clear to me. I hope Gurus might help me to understand this . I tried some additions and subtractions but could not find the storage pattern. Thanks for attention BR jaswi Tushar@cfd and wenxu like this.

 July 2, 2009, 03:52 #2 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Göttingen, Niedersachsen, Germany Posts: 280 Rep Power: 16 Jaswinder, just a quick pointer: Have a look into the ODE solver classes and how scalarMatrix is used there. Henrik

 July 2, 2009, 06:53 #3 Senior Member   Join Date: Mar 2009 Posts: 248 Rep Power: 16 Hi Henrik Thanks for the tip. I will follow that clue and post back . Best Regards jaswi

 July 2, 2009, 07:52 #4 Super Moderator     Niklas Nordin Join Date: Mar 2009 Location: Stockholm, Sweden Posts: 694 Rep Power: 27 a third variant i tried it in octave for fun... [l,u,p] = lu(hmm) and I got this l = 1.00000 0.00000 0.00000 -0.66667 1.00000 0.00000 -0.66667 0.76316 1.00000 u = -3.00000 10.00000 -4.00000 0.00000 12.66667 -1.66667 0.00000 0.00000 8.60526 p = 1 0 0 0 0 1 0 1 0 This is what 'help lu' says: Compute the LU decomposition of A, using subroutines from LAPACK.

 July 8, 2009, 04:37 #5 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Göttingen, Niedersachsen, Germany Posts: 280 Rep Power: 16 Dear Jaswi, you need to understand what pivotIndices is for, which takes the value (2 2 2) in my case. Have a look into Numerical Recipes (Press et al.). Henrik