songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

ENH: refactor SVD to support tensors directly (#3313)

Browse Source
This commit is contained in:
Mark Olesen 2025-03-26 11:11:11 +01:00
parent a2df607998
commit e3c93a9c85
6 changed files with 195 additions and 150 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

158
src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -30,28 +31,87 @@ License
#include "scalarMatrices.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * //
#include "tensor.H"
#include "diagTensor.H"
Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
:
U_(A),
V_(A.n(), A.n()),
S_(A.n()),
converged_(true),
nZeros_(0)
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Same as implementation in scalarMatrices, but for tensor/diagTensor/tensor
static tensor do_multiply
(
const tensor& A,
const diagTensor& B,
const tensor& C
)
{
constexpr direction size = 3;
tensor ans(Foam::zero{});
for (direction i = 0; i < size; ++i)
{
for (direction g = 0; g < size; ++g)
{
for (direction l = 0; l < size; ++l)
{
ans(i, g) += C(l, g)*A(i, l)*B[l];
}
}
}
return ans;
}
template<class MatrixType, class DiagMatrixType, class WorkArrayType>
static std::pair<label, bool> SVDcomp
(
// input
const MatrixType& A,
// input
const scalar minCondition,
// output
MatrixType& U_,
// output
MatrixType& V_,
// output
DiagMatrixType& S_,
// scratch space
WorkArrayType& rv1
)
{
label nZeros_(0);
bool converged_(true);
// SVDcomp to find U_, V_ and S_ - the singular values
U_ = A;
const label Un = U_.n();
const label Um = U_.m();
scalarList rv1(Un);
scalar g = 0;
scalar scale = 0;
scalar s = 0;
scalar anorm = 0;
label l = 0;
const auto sign = [](scalar a, scalar b) -> scalar
{
//return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
return b >= 0 ? a : -a;
};
for (label i=0; i<Un; i++)
{
l = i + 2;
@ -87,7 +147,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
}
f = s/h;
for (label k=i; k<A.m(); k++)
for (label k=i; k<Um; k++)
{
U_(k, j) += f*U_(k, i);
}
@ -373,23 +433,91 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
// Zero singular values that are less than minCondition*maxS
const scalar minS = minCondition*S_[findMax(S_)];
forAll(S_, i)
for (auto& val : S_)
{
if (S_[i] <= minS)
if (val <= minS)
{
S_[i] = 0;
nZeros_++;
val = 0;
++nZeros_;
}
}
return { nZeros_, converged_ };
}
} // End namespace Foam
// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * //
Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
:
U_(),
V_(A.n(), A.n()),
S_(A.n()),
converged_(true),
nZeros_(0)
{
scalarList rv1(A.n());
// SVDcomp to find U_, V_ and S_ - the singular values
auto status = SVDcomp(A, minCondition, U_, V_, S_, rv1);
nZeros_ = status.first;
converged_ = status.second;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::scalar Foam::SVD::minNonZeroS() const
{
scalar minS = S_[0];
for (label i=1; i<S_.size(); i++)
{
scalar s = S_[i];
if (s > VSMALL && s < minS) minS = s;
}
return minS;
}
Foam::scalarRectangularMatrix Foam::SVD::VSinvUt() const
{
scalarRectangularMatrix VSinvUt;
multiply(VSinvUt, V_, inv(S_), U_.T());
multiply(VSinvUt, V_, Foam::inv(S_), U_.T());
return VSinvUt;
}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
Foam::scalarRectangularMatrix Foam::SVD::pinv
(
const scalarRectangularMatrix& A,
const scalar minCondition
)
{
SVD svd(A, minCondition);
return svd.VSinvUt();
}
Foam::Tensor<Foam::scalar>
Foam::SVD::pinv(const Tensor<scalar>& A, const scalar minCondition)
{
// SVDcomp to find U_, V_ and S_ - the singular values
tensor U_;
tensor V_;
diagTensor S_;
FixedList<scalar, 3> rv1;
SVDcomp(A, minCondition, U_, V_, S_, rv1);
return do_multiply(V_, Foam::inv(S_), U_.T());
}
// ************************************************************************* //

64
src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -30,13 +31,12 @@ Description
Singular value decomposition of a rectangular matrix.
SourceFiles
SVDI.H
SVD.C
\*---------------------------------------------------------------------------*/
#ifndef SVD_H
#define SVD_H
#ifndef Foam_SVD_H
#define Foam_SVD_H
#include "scalarMatrices.H"
@ -45,7 +45,8 @@ SourceFiles
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Forward Declarations
template<class Cmpt> class Tensor;
/*---------------------------------------------------------------------------*\
Class SVD Declaration
@ -53,7 +54,7 @@ namespace Foam
class SVD
{
// Private data
// Private Data
//- Rectangular matrix with the same dimensions as the input
scalarRectangularMatrix U_;
@ -71,7 +72,9 @@ class SVD
label nZeros_;
// Private Member Functions
public:
// Generated Methods
//- No copy construct
SVD(const SVD&) = delete;
@ -79,40 +82,59 @@ class SVD
//- No copy assignment
void operator=(const SVD&) = delete;
template<class T>
inline const T sign(const T& a, const T& b);
public:
// Constructors
//- Construct from a rectangular Matrix
SVD(const scalarRectangularMatrix& A, const scalar minCondition = 0);
explicit SVD
(
const scalarRectangularMatrix& A,
const scalar minCondition = 0
);
// Access functions
// Access Functions
//- Return U
inline const scalarRectangularMatrix& U() const;
const scalarRectangularMatrix& U() const noexcept { return U_; }
//- Return the square matrix V
inline const scalarRectangularMatrix& V() const;
const scalarRectangularMatrix& V() const noexcept { return V_; }
//- Return the singular values
inline const scalarDiagonalMatrix& S() const;
const scalarDiagonalMatrix& S() const noexcept { return S_; }
//- Return the minimum non-zero singular value
inline bool converged() const;
bool converged() const noexcept { return converged_; }
//- Return the number of zero singular values
inline label nZeros() const;
label nZeros() const noexcept { return nZeros_; }
//- Return the minimum non-zero singular value
inline scalar minNonZeroS() const;
scalar minNonZeroS() const;
// Member Functions
//- Return the matrix product V S^(-1) U^T (the pseudo inverse)
scalarRectangularMatrix VSinvUt() const;
// Static Member Functions
//- Return the pseudo inverse of the given matrix
static scalarRectangularMatrix pinv
(
const scalarRectangularMatrix& A,
const scalar minCondition = 0
);
//- Return the pseudo inverse of the given tensor
static Tensor<scalar> pinv
(
const Tensor<scalar>& A,
const scalar minCondition = 0
);
};
@ -122,10 +144,6 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "SVDI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

82
src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
View File

@ -1,82 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
-------------------------------------------------------------------------------
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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * * //
template<class T>
inline const T Foam::SVD::sign(const T& a, const T& b)
{
//return b >= 0 ? (a >= 0 ? a : -a) : (a >= 0 ? -a : a);
return b >= 0 ? a : -a;
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::scalarRectangularMatrix& Foam::SVD::U() const
{
return U_;
}
inline const Foam::scalarRectangularMatrix& Foam::SVD::V() const
{
return V_;
}
inline const Foam::scalarDiagonalMatrix& Foam::SVD::S() const
{
return S_;
}
inline bool Foam::SVD::converged() const
{
return converged_;
}
inline Foam::label Foam::SVD::nZeros() const
{
return nZeros_;
}
inline Foam::scalar Foam::SVD::minNonZeroS() const
{
scalar minS = S_[0];
for (label i=1; i<S_.size(); i++)
{
scalar s = S_[i];
if (s > VSMALL && s < minS) minS = s;
}
return minS;
}
// ************************************************************************* //

7
src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
View File

@ -271,7 +271,7 @@ void Foam::multiply
<< abort(FatalError);
}
ans = scalarRectangularMatrix(A.m(), C.n(), Zero);
ans = scalarRectangularMatrix(A.m(), C.n(), Foam::zero{});
for (label i = 0; i < A.m(); ++i)
{
@ -312,7 +312,7 @@ void Foam::multiply
const label size = A.m();
ans = scalarSquareMatrix(size, Zero);
ans = scalarSquareMatrix(size, Foam::zero{});
for (label i = 0; i < size; ++i)
{
@ -334,8 +334,7 @@ Foam::scalarRectangularMatrix Foam::SVDinv
scalar minCondition
)
{
SVD svd(A, minCondition);
return svd.VSinvUt();
return SVD::pinv(A, minCondition);
}

17
src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2020-2022 OpenCFD Ltd.
Copyright (C) 2020-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -345,19 +345,14 @@ Foam::symmTensor Foam::pinv(const symmTensor& st)
if (detval < ROOTVSMALL)
{
// Fall back to pseudo inverse
scalarRectangularMatrix mat(3, 3);
mat(0,0) = st.xx(); mat(0,1) = st.xy(); mat(0,2) = st.xz();
mat(1,0) = st.yx(); mat(1,1) = st.yy(); mat(1,2) = st.yz();
mat(2,0) = st.zx(); mat(2,1) = st.zy(); mat(2,2) = st.zz();
mat = SVDinv(mat);
tensor mat(st);
tensor t = SVD::pinv(mat);
return symmTensor
(
mat(0,0), mat(0,1), mat(0,2),
mat(1,1), mat(1,2),
mat(2,2)
t.xx(), t.xy(), t.xz(),
t.yy(), t.yz(),
t.zz()
);
}

17
src/OpenFOAM/primitives/Tensor/floats/tensor.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2020-2022 OpenCFD Ltd.
Copyright (C) 2020-2025 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -309,20 +309,7 @@ Foam::tensor Foam::pinv(const tensor& t)
if (detval < ROOTVSMALL)
{
// Fall back to pseudo inverse
scalarRectangularMatrix mat(3, 3);
mat(0,0) = t.xx(); mat(0,1) = t.xy(); mat(0,2) = t.xz();
mat(1,0) = t.yx(); mat(1,1) = t.yy(); mat(1,2) = t.yz();
mat(2,0) = t.zx(); mat(2,1) = t.zy(); mat(2,2) = t.zz();
mat = SVDinv(mat);
return tensor
(
mat(0,0), mat(0,1), mat(0,2),
mat(1,0), mat(1,1), mat(1,2),
mat(2,0), mat(2,1), mat(2,2)
);
return SVD::pinv(t);
}
return Foam::inv(t, detval);