- **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*)

LUDecompose, LUBacksubstitue in class scalarMatrixDear 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[])` 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:
Thanks for attention BR jaswi |

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

Hi Henrik
Thanks for the tip. I will follow that clue and post back . Best Regards jaswi |

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. |

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 23:16. |