|
[Sponsors] |
LUDecompose, LUBacksubstitue in class scalarMatrix |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
July 1, 2009, 05:43 |
LUDecompose, LUBacksubstitue in class scalarMatrix
|
#1 |
Senior Member
Join Date: Mar 2009
Posts: 248
Rep Power: 18 |
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; } 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 Thanks for attention BR jaswi |
|
July 2, 2009, 03:52 |
|
#2 |
Senior Member
Henrik Rusche
Join Date: Mar 2009
Location: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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: 18 |
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: 693
Rep Power: 29 |
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: Wernigerode, Sachsen-Anhalt, Germany
Posts: 281
Rep Power: 18 |
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 | Search this Thread |
Display Modes | |
|
|
Similar Threads | ||||
Thread | Thread Starter | Forum | Replies | Last Post |
STAR-CCM+ Java class | Claudene | Siemens | 1 | July 17, 2008 06:49 |
Errors running allwmake in OpenFOAM141dev with WM_COMPILE_OPTION%3ddebug | unoder | OpenFOAM Installation | 11 | January 30, 2008 21:30 |
How to add a new class locally | ville | OpenFOAM | 4 | December 11, 2006 14:20 |
About UList and List class | leosding | OpenFOAM Running, Solving & CFD | 1 | December 2, 2005 00:52 |
Expanding a class | fabianpk | OpenFOAM | 0 | October 3, 2005 05:26 |