Commit 857915da authored by Henry Weller's avatar Henry Weller
Browse files

LUscalarMatrix: Added processor-local matrix inverse function

parent f4a588d4
......@@ -152,6 +152,11 @@ int main(int argc, char *argv[])
LUscalarMatrix LU(squareMatrix);
scalarField x(LU.solve(source));
Info<< "LU solve residual " << (squareMatrix*x - source) << endl;
scalarSquareMatrix inv(3);
LU.inv(inv);
Info<< "LU inv " << inv << endl;
Info<< "LU inv*squareMatrix " << (inv*squareMatrix) << endl;
}
{
......@@ -169,6 +174,8 @@ int main(int argc, char *argv[])
Info<< "QR inverse solve residual "
<< (x - QR.inv()*source) << endl;
Info<< "QR inv *squareMatrix " << (QR.inv()*squareMatrix) << endl;
}
Info<< "\nEnd\n" << endl;
......
......@@ -412,4 +412,21 @@ void Foam::LUscalarMatrix::decompose(const scalarSquareMatrix& M)
}
void Foam::LUscalarMatrix::inv(scalarSquareMatrix& M) const
{
scalarField source(m());
for (label j=0; j<m(); j++)
{
source = Zero;
source[j] = 1;
LUBacksubstitute(*this, pivotIndices_, source);
for (label i=0; i<m(); i++)
{
M(i, j) = source[i];
}
}
}
// ************************************************************************* //
......@@ -126,6 +126,9 @@ public:
// returning the solution
template<class Type>
tmp<Field<Type>> solve(const Field<Type>& source) const;
//- Set M to the inverse of this square matrix
void inv(scalarSquareMatrix& M) const;
};
......
......@@ -269,7 +269,7 @@ void Foam::multiply
for (label i=0; i<A.m(); i++)
{
for (label g = 0; g < C.n(); g++)
for (label g=0; g<C.n(); g++)
{
for (label l=0; l<C.m(); l++)
{
......@@ -280,6 +280,47 @@ void Foam::multiply
}
void Foam::multiply
(
scalarSquareMatrix& ans, // value changed in return
const scalarSquareMatrix& A,
const DiagonalMatrix<scalar>& B,
const scalarSquareMatrix& C
)
{
if (A.m() != B.size())
{
FatalErrorInFunction
<< "A and B must have identical dimensions but A.m = "
<< A.m() << " and B.m = " << B.size()
<< abort(FatalError);
}
if (B.size() != C.m())
{
FatalErrorInFunction
<< "B and C must have identical dimensions but B.m = "
<< B.size() << " and C.m = " << C.m()
<< abort(FatalError);
}
const label size = A.m();
ans = scalarSquareMatrix(size, Zero);
for (label i=0; i<size; i++)
{
for (label g=0; g<size; g++)
{
for (label l=0; l<size; l++)
{
ans(i, g) += C(l, g)*A(i, l)*B[l];
}
}
}
}
Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
(
const scalarRectangularMatrix& A,
......
......@@ -145,6 +145,14 @@ void multiply
const scalarRectangularMatrix& C
);
void multiply
(
scalarSquareMatrix& answer, // value changed in return
const scalarSquareMatrix& A,
const DiagonalMatrix<scalar>& B,
const scalarSquareMatrix& C
);
//- Return the inverse of matrix A using SVD
scalarRectangularMatrix SVDinv
(
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment