Commit 3b7de017 authored by Henry Weller's avatar Henry Weller
Browse files

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'.
parent af6b67c7
......@@ -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<Un; k++)
{
......@@ -281,9 +281,9 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
for (label j=0; j<Um; j++)
{
scalar y = U_[j][mn];
scalar y = U_(j, mn);
scalar z = U_(j, i);
U_[j][mn] = y*c + z*s;
U_(j, mn) = y*c + z*s;
U_(j, i) = z*c - y*s;
}
}
......@@ -303,11 +303,10 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
}
break;
}
if (its == 29)
{
WarningInFunction
<< "No convergence in 30 SVD iterations"
<< endl;
converged_ = false;
}
scalar x = S_[l];
......@@ -339,10 +338,10 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
for (label jj = 0; jj < Un; jj++)
{
x = V_[jj][j];
z = V_[jj][i];
V_[jj][j] = x*c + z*s;
V_[jj][i] = z*c - x*s;
x = V_(jj, j);
z = V_(jj, i);
V_(jj, j) = x*c + z*s;
V_(jj, i) = z*c - x*s;
}
z = sqrtSumSqr(f, h);
......@@ -358,10 +357,10 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
for (label jj=0; jj < Um; jj++)
{
y = U_[jj][j];
z = U_[jj][i];
U_[jj][j] = y*c + z*s;
U_[jj][i] = z*c - y*s;
y = U_(jj, j);
z = U_(jj, i);
U_(jj, j) = y*c + z*s;
U_(jj, i) = z*c - y*s;
}
}
rv1[l] = 0;
......@@ -380,9 +379,14 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
nZeros_++;
}
}
}
// Now multiply out to find the pseudo inverse of A, VSinvUt_
multiply(VSinvUt_, V_, inv(S_), U_.T());
Foam::scalarRectangularMatrix Foam::SVD::VSinvUt() const
{
scalarRectangularMatrix VSinvUt;
multiply(VSinvUt, V_, inv(S_), U_.T());
return VSinvUt;
}
......
......@@ -62,8 +62,8 @@ class SVD
//- The singular values
DiagonalMatrix<scalar> 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;
};
......
......@@ -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];
......
......@@ -321,7 +321,7 @@ void Foam::multiply
}
Foam::RectangularMatrix<Foam::scalar> Foam::SVDinv
Foam::scalarRectangularMatrix Foam::SVDinv
(
const scalarRectangularMatrix& A,
scalar minCondition
......
......@@ -134,17 +134,18 @@ void Foam::CentredFitSnGradData<Polynomial>::calcFit
for (int iIt = 0; iIt < 8 && !goodFit; iIt++)
{
SVD svd(B, SMALL);
scalarRectangularMatrix invB(svd.VSinvUt());
for (label i=0; i<stencilSize; i++)
{
coeffsi[i] = wts[1]*wts[i]*svd.VSinvUt()(1, i)/scale;
coeffsi[i] = wts[1]*wts[i]*invB(1, i)/scale;
}
goodFit =
(
mag(wts[0]*wts[0]*svd.VSinvUt()(0, 0) - wLin)
mag(wts[0]*wts[0]*invB(0, 0) - wLin)
< this->linearLimitFactor()*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<Polynomial>::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;
......
......@@ -208,13 +208,14 @@ void Foam::FitData<FitDataType, ExtendedStencil, Polynomial>::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<stencilSize; i++)
{
coeffsi[i] = wts[0]*wts[i]*svd.VSinvUt()(0, i);
coeffsi[i] = wts[0]*wts[i]*invB(0, i);
if (mag(coeffsi[i]) > maxCoeff)
{
maxCoeff = mag(coeffsi[i]);
......
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