diff --git a/applications/test/Matrix/Test-Matrix.C b/applications/test/Matrix/Test-Matrix.C index 9663ae8217..380576d695 100644 --- a/applications/test/Matrix/Test-Matrix.C +++ b/applications/test/Matrix/Test-Matrix.C @@ -116,17 +116,19 @@ int main(int argc, char *argv[]) squareMatrix[2][1] = -43; squareMatrix[2][2] = 98; + const scalarSquareMatrix squareMatrixCopy = squareMatrix; Info<< nl << "Square Matrix = " << squareMatrix << endl; - scalarDiagonalMatrix rhs(3, 0); - rhs[0] = 1; - rhs[1] = 2; - rhs[2] = 3; + Info<< "det = " << det(squareMatrixCopy) << endl; - LUsolve(squareMatrix, rhs); + labelList rhs(3, 0); + label sign; + LUDecompose(squareMatrix, rhs, sign); Info<< "Decomposition = " << squareMatrix << endl; - Info<< "Solution = " << rhs << endl; + Info<< "Pivots = " << rhs << endl; + Info<< "Sign = " << sign << endl; + Info<< "det = " << detDecomposed(squareMatrix, sign) << endl; } Info<< "\nEnd\n" << endl; diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C new file mode 100644 index 0000000000..439257fcb6 --- /dev/null +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.C @@ -0,0 +1,73 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "SquareMatrix.H" +#include "labelList.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +template +Foam::scalar Foam::detDecomposed +( + const SquareMatrix& matrix, + const label sign +) +{ + scalar diagProduct = 1.0; + + for (label i = 0; i < matrix.n(); ++i) + { + diagProduct *= matrix[i][i]; + } + + return sign*diagProduct; +} + + +template +Foam::scalar Foam::det(const SquareMatrix& matrix) +{ + SquareMatrix matrixTmp = matrix; + + labelList pivotIndices(matrix.n()); + label sign; + LUDecompose(matrixTmp, pivotIndices, sign); + + return detDecomposed(matrixTmp, sign); +} + + +template +Foam::scalar Foam::det(SquareMatrix& matrix) +{ + labelList pivotIndices(matrix.n()); + label sign; + LUDecompose(matrix, pivotIndices, sign); + + return detDecomposed(matrix, sign); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H index 660f39321e..531b649588 100644 --- a/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H +++ b/src/OpenFOAM/matrices/SquareMatrix/SquareMatrix.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -81,6 +81,21 @@ public: }; +// Global functions + +//- Return the LU decomposed SquareMatrix det +template +scalar detDecomposed(const SquareMatrix&, const label sign); + +//- Return the SquareMatrix det +template +scalar det(const SquareMatrix&); + +//- Return the SquareMatrix det and the LU decomposition in the original matrix +template +scalar det(SquareMatrix&); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam @@ -91,6 +106,12 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +#ifdef NoRepository +# include "SquareMatrix.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + #endif // ************************************************************************* // diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C index 68a6955d71..87a9d02347 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -33,9 +33,22 @@ void Foam::LUDecompose scalarSquareMatrix& matrix, labelList& pivotIndices ) +{ + label sign; + LUDecompose(matrix, pivotIndices, sign); +} + + +void Foam::LUDecompose +( + scalarSquareMatrix& matrix, + labelList& pivotIndices, + label& sign +) { label n = matrix.n(); scalar vv[n]; + sign = 1; for (register label i=0; i