CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   LUDecompose, LUBacksubstitue in class scalarMatrix (https://www.cfd-online.com/Forums/openfoam-programming-development/65925-ludecompose-lubacksubstitue-class-scalarmatrix.html)

 jaswi July 1, 2009 04:43

LUDecompose, LUBacksubstitue in class scalarMatrix

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

 henrik July 2, 2009 02:52

Jaswinder,

just a quick pointer: Have a look into the ODE solver classes and how scalarMatrix is used there.

Henrik

 jaswi July 2, 2009 05:53

Hi Henrik

Thanks for the tip. I will follow that clue and post back .

Best Regards
jaswi

 niklas July 2, 2009 06:52

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.

 henrik July 8, 2009 03:37

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

 All times are GMT -4. The time now is 12:32.