Commit 3577b3d4 authored by laurence's avatar laurence
Browse files

ENH: Add SymmetricSquareMatrix class and specialisations of LUDecompose and LUSolve.

parent 448bd7be
......@@ -23,7 +23,7 @@ License
\*---------------------------------------------------------------------------*/
#include "SquareMatrix.H"
#include "scalarMatrices.H"
#include "vector.H"
using namespace Foam;
......@@ -49,15 +49,15 @@ int main(int argc, char *argv[])
Info<< max(hmm) << endl;
Info<< min(hmm) << endl;
SquareMatrix<scalar> hmm2(3, 1.0);
SquareMatrix<scalar> hmm2(3, 3, 1.0);
hmm = hmm2;
Info<< hmm << endl;
SquareMatrix<scalar> hmm3(Sin);
//SquareMatrix<scalar> hmm3(Sin);
Info<< hmm3 << endl;
//Info<< hmm3 << endl;
SquareMatrix<scalar> hmm4;
......@@ -70,7 +70,65 @@ int main(int argc, char *argv[])
hmm4 = hmm5;
Info<< hmm5 << endl;
Info<< "End\n" << endl;
{
scalarSymmetricSquareMatrix symmMatrix(3, 3, 0);
symmMatrix(0, 0) = 4;
symmMatrix(1, 0) = 12;
symmMatrix(1, 1) = 37;
symmMatrix(2, 0) = -16;
symmMatrix(2, 1) = -43;
symmMatrix(2, 2) = 98;
Info<< "Symmetric Square Matrix = " << symmMatrix << endl;
Info<< "Inverse = " << inv(symmMatrix) << endl;
Info<< "Determinant = " << det(symmMatrix) << endl;
scalarSymmetricSquareMatrix symmMatrix2(symmMatrix);
LUDecompose(symmMatrix2);
Info<< "Inverse = " << invDecomposed(symmMatrix2) << endl;
Info<< "Determinant = " << detDecomposed(symmMatrix2) << endl;
scalarDiagonalMatrix rhs(3, 0);
rhs[0] = 1;
rhs[1] = 2;
rhs[2] = 3;
LUsolve(symmMatrix, rhs);
Info<< "Decomposition = " << symmMatrix << endl;
Info<< "Solution = " << rhs << endl;
}
{
scalarSquareMatrix squareMatrix(3, 3, 0);
squareMatrix[0][0] = 4;
squareMatrix[0][1] = 12;
squareMatrix[0][2] = -16;
squareMatrix[1][0] = 12;
squareMatrix[1][1] = 37;
squareMatrix[1][2] = -43;
squareMatrix[2][0] = -16;
squareMatrix[2][1] = -43;
squareMatrix[2][2] = 98;
Info<< nl << "Square Matrix = " << squareMatrix << endl;
scalarDiagonalMatrix rhs(3, 0);
rhs[0] = 1;
rhs[1] = 2;
rhs[2] = 3;
LUsolve(squareMatrix, rhs);
Info<< "Decomposition = " << squareMatrix << endl;
Info<< "Solution = " << rhs << endl;
}
Info<< "\nEnd\n" << endl;
return 0;
}
......
......@@ -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
......@@ -76,7 +76,7 @@ public:
// Member functions
//- Invert the diaganol matrix and return itself
//- Invert the diagonal matrix and return itself
DiagonalMatrix<Type>& invert();
};
......
......@@ -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
......@@ -221,7 +221,6 @@ template<class Form, class Type> Form operator*
const Matrix<Form, Type>&
);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
......
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "SymmetricSquareMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
template<class Type>
Foam::SymmetricSquareMatrix<Type> Foam::invDecomposed
(
const SymmetricSquareMatrix<Type>& matrix
)
{
SymmetricSquareMatrix<Type> inv(matrix.n(), matrix.n(), 0.0);
for (label i = 0; i < matrix.n(); ++i)
{
inv[i][i] = 1.0/matrix[i][i];
for (label j = 0; j < i; ++j)
{
scalar sum = 0.0;
for (label k = j; k < i; k++)
{
sum -= matrix[i][k]*inv[k][j];
}
inv[i][j] = sum/matrix[i][i];
}
}
return inv.T()*inv;
}
template<class Type>
Foam::SymmetricSquareMatrix<Type> Foam::inv
(
const SymmetricSquareMatrix<Type>& matrix
)
{
SymmetricSquareMatrix<Type> matrixTmp(matrix);
LUDecompose(matrixTmp);
return invDecomposed(matrixTmp);
}
template<class Type>
Foam::scalar Foam::detDecomposed(const SymmetricSquareMatrix<Type>& matrix)
{
scalar diagProduct = 1.0;
for (label i = 0; i < matrix.n(); ++i)
{
diagProduct *= matrix[i][i];
}
return sqr(diagProduct);
}
template<class Type>
Foam::scalar Foam::det(const SymmetricSquareMatrix<Type>& matrix)
{
SymmetricSquareMatrix<Type> matrixTmp = matrix;
LUDecompose(matrixTmp);
return detDecomposed(matrixTmp);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>.
Class
Foam::SymmetricSquareMatrix
Description
A templated 2D square symmetric matrix of objects of \<T\>, where the
n x n matrix dimension is known and used for subscript bounds checking, etc.
SourceFiles
SymmetricSquareMatrixI.H
SymmetricSquareMatrix.C
\*---------------------------------------------------------------------------*/
#ifndef SymmetricSquareMatrix_H
#define SymmetricSquareMatrix_H
#include "SquareMatrix.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class SymmetricSquareMatrix Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class SymmetricSquareMatrix
:
public Matrix<SymmetricSquareMatrix<Type>, Type>
{
public:
// Constructors
//- Null constructor.
inline SymmetricSquareMatrix();
//- Construct given number of rows/columns.
inline SymmetricSquareMatrix(const label n);
//- Construct with given number of rows/columns
inline SymmetricSquareMatrix(const label m, const label n);
//- Construct with given number of rows/columns
// and value for all elements.
inline SymmetricSquareMatrix(const label m, const label n, const Type&);
//- Construct from Istream.
inline SymmetricSquareMatrix(Istream&);
//- Clone
inline autoPtr<SymmetricSquareMatrix<Type> > clone() const;
//- Return subscript-checked row of Matrix.
inline Type& operator()(const label r, const label c);
//- Return subscript-checked row of constant Matrix.
inline const Type& operator()(const label r, const label c) const;
};
// Global functions
//- Return the LU decomposed SymmetricSquareMatrix inverse
template<class Type>
SymmetricSquareMatrix<Type> invDecomposed(const SymmetricSquareMatrix<Type>&);
//- Return the SymmetricSquareMatrix inverse
template<class Type>
SymmetricSquareMatrix<Type> inv(const SymmetricSquareMatrix<Type>&);
//- Return the LU decomposed SymmetricSquareMatrix det
template<class Type>
scalar detDecomposed(const SymmetricSquareMatrix<Type>&);
//- Return the SymmetricSquareMatrix det
template<class Type>
scalar det(const SymmetricSquareMatrix<Type>&);
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
# include "SymmetricSquareMatrixI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "SymmetricSquareMatrix.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix()
:
Matrix<SymmetricSquareMatrix<Type>, Type>()
{}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(const label n)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(n, n)
{}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label m,
const label n
)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(m, n)
{
if (m != n)
{
FatalErrorIn
(
"SymmetricSquareMatrix<Type>::SymmetricSquareMatrix"
"(const label m, const label n)"
) << "m != n for constructing a symmetric square matrix"
<< exit(FatalError);
}
}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix
(
const label m,
const label n,
const Type& t
)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(m, n, t)
{
if (m != n)
{
FatalErrorIn
(
"SymmetricSquareMatrix<Type>::SymmetricSquareMatrix"
"(const label m, const label n, const Type&)"
) << "m != n for constructing a symmetric square matrix"
<< exit(FatalError);
}
}
template<class Type>
inline Foam::SymmetricSquareMatrix<Type>::SymmetricSquareMatrix(Istream& is)
:
Matrix<SymmetricSquareMatrix<Type>, Type>(is)
{}
template<class Type>
inline Foam::autoPtr<Foam::SymmetricSquareMatrix<Type> >
Foam::SymmetricSquareMatrix<Type>::clone() const
{
return autoPtr<SymmetricSquareMatrix<Type> >
(
new SymmetricSquareMatrix<Type>(*this)
);
}
template<class Type>
inline Type& Foam::SymmetricSquareMatrix<Type>::operator()
(
const label r,
const label c
)
{
if (r > c)
{
return this->operator[](r)[c];
}
else
{
return this->operator[](c)[r];
}
}
template<class Type>
inline const Type& Foam::SymmetricSquareMatrix<Type>::operator()
(
const label r,
const label c
) const
{
if (r > c)
{
return this->operator[](r)[c];
}
else
{
return this->operator[](c)[r];
}
}
// ************************************************************************* //
......@@ -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
......@@ -134,6 +134,56 @@ void Foam::LUDecompose
}
void Foam::LUDecompose(scalarSymmetricSquareMatrix& matrix)
{
// Store result in upper triangular part of matrix
label size = matrix.n();
// Set upper triangular parts to zero.
for (label j = 0; j < size; j++)
{
for (label k = j + 1; k < size; k++)
{
matrix[j][k] = 0.0;
}
}
for (label j = 0; j < size; j++)
{
scalar d = 0.0;
for (label k = 0; k < j; k++)
{
scalar s = 0.0;
for (label i = 0; i < k; i++)
{
s += matrix[i][k]*matrix[i][j];
}
s = (matrix[j][k] - s)/matrix[k][k];
matrix[k][j] = s;
matrix[j][k] = s;
d += sqr(s);
}
d = matrix[j][j] - d;
if (d < 0.0)
{
FatalErrorIn("Foam::LUDecompose(scalarSymmetricSquareMatrix&)")
<< "Matrix is not symmetric positive-definite. Unable to "
<< "decompose."
<< abort(FatalError);
}
matrix[j][j] = sqrt(d);
}
}
// * * * * * * * * * * * * * * * Global Functions * * * * * * * * * * * * * //
void Foam::multiply
......
......@@ -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
......@@ -27,6 +27,10 @@ Class
Description
Scalar matrices
LUDecompose for scalarSymmetricSquareMatrix implements the Cholesky
decomposition method from JAMA, a public-domain library developed at NIST,
available at http://math.nist.gov/tnt/index.html
SourceFiles
scalarMatrices.C
scalarMatricesTemplates.C
......@@ -38,6 +42,7 @@ SourceFiles
#include "RectangularMatrix.H"
#include "SquareMatrix.H"
#include "SymmetricSquareMatrix.H"
#include "DiagonalMatrix.H"
#include "scalarField.H"
#include "labelList.H"
......@@ -49,21 +54,22 @@ namespace Foam
typedef RectangularMatrix<scalar> scalarRectangularMatrix;
typedef SquareMatrix<scalar> scalarSquareMatrix;
typedef SymmetricSquareMatrix<scalar> scalarSymmetricSquareMatrix;
typedef DiagonalMatrix<scalar> scalarDiagonalMatrix;
//- Solve the matrix using Gaussian elimination with pivoting,
// returning the solution in the source
template<class Type>
void solve(scalarSquareMatrix& matrix, Field<Type>& source);
void solve(scalarSquareMatrix& matrix, List<Type>& source);
//- Solve the matrix using Gaussian elimination with pivoting
// and return the solution
template<class Type>
void solve
(