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<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;
 }
 
 
diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
index 2f4a9dcb19..a384f22920 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.H
@@ -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;
 };
 
 
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::scalar> 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<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;
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<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]);