From 3b7de01705491ed1cdd35f8045e8e819b8c241d9 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Mon, 18 Jul 2016 13:15:25 +0100 Subject: [PATCH] SVD: VSinvUt is not always needed and is now calculated and returned by the 'VSinvUt()' function rather than being calculated on construction and stored as member data. The convergence warning has be replaced with the 'convergence()' member function which returns 'true' if the SVD iteration converged, otherwise 'false'. --- .../matrices/scalarMatrices/SVD/SVD.C | 40 ++++++++++--------- .../matrices/scalarMatrices/SVD/SVD.H | 11 +++-- .../matrices/scalarMatrices/SVD/SVDI.H | 9 ++++- .../matrices/scalarMatrices/scalarMatrices.C | 2 +- .../CentredFitSnGrad/CentredFitSnGradData.C | 15 +++---- .../schemes/FitData/FitData.C | 3 +- 6 files changed, 47 insertions(+), 33 deletions(-) diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C index 9f9c7a865e..8b090ff275 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C @@ -35,7 +35,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) U_(A), V_(A.n(), A.n()), S_(A.n()), - VSinvUt_(A.n(), A.m()), + converged_(true), nZeros_(0) { // SVDcomp to find U_, V_ and S_ - the singular values @@ -117,10 +117,10 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition) s += U_(i, k)*U_(i, k); } - scalar f = U_[i][l-1]; + scalar f = U_(i, l-1); g = -sign(sqrt(s), f); scalar h = f*g - s; - U_[i][l-1] = f - g; + U_(i, l-1) = f - g; for (label k=l-1; k S_; - //- The matrix product V S^(-1) U^T - scalarRectangularMatrix VSinvUt_; + //- Convergence flag + bool converged_; //- The number of zero singular values label nZeros_; @@ -100,14 +100,17 @@ public: //- Return the singular values inline const scalarDiagonalMatrix& S() const; - //- Return VSinvUt (the pseudo inverse) - inline const scalarRectangularMatrix& VSinvUt() const; + //- Return the minimum non-zero singular value + inline bool converged() const; //- Return the number of zero singular values inline label nZeros() const; //- Return the minimum non-zero singular value inline scalar minNonZeroS() const; + + //- Return the matrix product V S^(-1) U^T (the pseudo inverse) + scalarRectangularMatrix VSinvUt() const; }; diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H index 9ff7ade03e..65bc0bd94a 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H +++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H @@ -40,26 +40,31 @@ 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 const Foam::scalarRectangularMatrix& Foam::SVD::VSinvUt() const + +inline bool Foam::SVD::converged() const { - return VSinvUt_; + return converged_; } + inline Foam::label Foam::SVD::nZeros() const { return nZeros_; } + inline Foam::scalar Foam::SVD::minNonZeroS() const { scalar minS = S_[0]; diff --git a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C index 0eb52a1be7..bbdb5118f1 100644 --- a/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C +++ b/src/OpenFOAM/matrices/scalarMatrices/scalarMatrices.C @@ -321,7 +321,7 @@ void Foam::multiply } -Foam::RectangularMatrix Foam::SVDinv +Foam::scalarRectangularMatrix Foam::SVDinv ( const scalarRectangularMatrix& A, scalar minCondition diff --git a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C index 519bf593d1..5c5980d9b2 100644 --- a/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C +++ b/src/finiteVolume/finiteVolume/snGradSchemes/CentredFitSnGrad/CentredFitSnGradData.C @@ -134,17 +134,18 @@ void Foam::CentredFitSnGradData::calcFit for (int iIt = 0; iIt < 8 && !goodFit; iIt++) { SVD svd(B, SMALL); + scalarRectangularMatrix invB(svd.VSinvUt()); for (label i=0; ilinearLimitFactor()*wLin) - && (mag(wts[0]*wts[1]*svd.VSinvUt()(0, 1) - (1 - wLin) + && (mag(wts[0]*wts[1]*invB(0, 1) - (1 - wLin) ) < this->linearLimitFactor()*(1 - wLin)) && coeffsi[0] < 0 && coeffsi[1] > 0 && mag(coeffsi[0] + deltaCoeff) < 0.5*deltaCoeff @@ -163,10 +164,10 @@ void Foam::CentredFitSnGradData::calcFit << " deltaCoeff " << deltaCoeff << nl << " sing vals " << svd.S() << nl << "Components of goodFit:\n" - << " wts[0]*wts[0]*svd.VSinvUt()(0, 0) = " - << wts[0]*wts[0]*svd.VSinvUt()(0, 0) << nl - << " wts[0]*wts[1]*svd.VSinvUt()(0, 1) = " - << wts[0]*wts[1]*svd.VSinvUt()(0, 1) + << " wts[0]*wts[0]*invB(0, 0) = " + << wts[0]*wts[0]*invB(0, 0) << nl + << " wts[0]*wts[1]*invB(0, 1) = " + << wts[0]*wts[1]*invB(0, 1) << " dim = " << this->dim() << endl; wts[0] *= 10; diff --git a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C index eb6da09f19..6fe277786b 100644 --- a/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C +++ b/src/finiteVolume/interpolation/surfaceInterpolation/schemes/FitData/FitData.C @@ -208,13 +208,14 @@ void Foam::FitData::calcFit for (int iIt = 0; iIt < 8 && !goodFit; iIt++) { SVD svd(B, SMALL); + scalarRectangularMatrix invB(svd.VSinvUt()); scalar maxCoeff = 0; label maxCoeffi = 0; for (label i=0; i maxCoeff) { maxCoeff = mag(coeffsi[i]);