# LUDecompose, LUBacksubstitue in class scalarMatrix

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

 July 1, 2009, 04:43 LUDecompose, LUBacksubstitue in class scalarMatrix #1 Senior Member   Join Date: Mar 2009 Posts: 248 Rep Power: 11 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 likes this.

 July 2, 2009, 02:52 #2 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 277 Rep Power: 11 Jaswinder, just a quick pointer: Have a look into the ODE solver classes and how scalarMatrix is used there. Henrik

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

 July 2, 2009, 06:52 #4 Super Moderator     Niklas Nordin Join Date: Mar 2009 Location: Stockholm, Sweden Posts: 693 Rep Power: 22 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, 03:37 #5 Senior Member   Henrik Rusche Join Date: Mar 2009 Location: Braunschweig, Niedersachsen, Germany Posts: 277 Rep Power: 11 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

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Claudene Siemens 1 July 17, 2008 05:49 unoder OpenFOAM Installation 11 January 30, 2008 21:30 ville OpenFOAM 4 December 11, 2006 14:20 leosding OpenFOAM Running, Solving & CFD 1 December 2, 2005 00:52 fabianpk OpenFOAM 0 October 3, 2005 04:26

All times are GMT -4. The time now is 06:43.