CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   LUDecompose, LUBacksubstitue in class scalarMatrix (http://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 11:04.