diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
index 0f78aef26e6e283500e164397f70a2ec3c80cf39..9383f233a2335316240f26bc00db6e825b26a9e7 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
@@ -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());
+}
+
+
// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
index 6bb5c5050aa5dd358173600d5817befd1ea08a51..005656b54c3d523f809a8625615e43b99d408f22 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
@@ -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,50 +82,65 @@ 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
+ );
+};
-} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-#include "SVDI.H"
+} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
deleted file mode 100644
index f23c90714f0bf4140dc816ddb5281b70632d874a..0000000000000000000000000000000000000000
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
+++ /dev/null
@@ -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;
-}
-
-
-// ************************************************************************* //
diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
index 1bd42a36080dad51a0be66547ee16e168659252a..5d82783a0c9569192bc7ca92ca4c81c5c8a0539f 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C
@@ -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);
}
diff --git a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
index fb41001b1f8009a25c3a12572172d697bd4b794c..d8c5af408e59df177a5d6202c658f6d160b9b3bb 100644
--- a/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
+++ b/src/OpenFOAM/primitives/SymmTensor/symmTensor/symmTensor.C
@@ -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()
);
}
diff --git a/src/OpenFOAM/primitives/Tensor/floats/tensor.C b/src/OpenFOAM/primitives/Tensor/floats/tensor.C
index d26aad40a4147af6c05f97316a986534f555d6cf..e1b5dcff68ebe807f61381ebc0c0390109588f9f 100644
--- a/src/OpenFOAM/primitives/Tensor/floats/tensor.C
+++ b/src/OpenFOAM/primitives/Tensor/floats/tensor.C
@@ -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);