diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
index 4c833cd40..9f9c7a865 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVD.C
@@ -48,11 +48,11 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
     scalar scale = 0;
     scalar s = 0;
     scalar anorm = 0;
-    label l=0;
+    label l = 0;
 
     for (label i=0; i<Un; i++)
     {
-        l = i+2;
+        l = i + 2;
         rv1[i] = scale*g;
         g = s = scale = 0;
 
@@ -72,7 +72,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                 }
 
                 scalar f = U_(i, i);
-                g = -sign(Foam::sqrt(s), f);
+                g = -sign(sqrt(s), f);
                 scalar h = f*g - s;
                 U_(i, i) = f - g;
 
@@ -118,7 +118,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                 }
 
                 scalar f = U_[i][l-1];
-                g = -sign(Foam::sqrt(s),f);
+                g = -sign(sqrt(s), f);
                 scalar h = f*g - s;
                 U_[i][l-1] = f - g;
 
@@ -150,6 +150,8 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
         anorm = max(anorm, mag(S_[i]) + mag(rv1[i]));
     }
 
+    anorm *= SMALL;
+
     for (label i=Un-1; i >= 0; i--)
     {
         if (i < Un-1)
@@ -178,7 +180,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
             for (label j=l; j<Un;j++)
             {
-                V_(i, j) = V_(j, i) = 0.0;
+                V_(i, j) = V_(j, i) = 0;
             }
         }
 
@@ -194,7 +196,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
         for (label j=l; j<Un; j++)
         {
-            U_(i, j) = 0.0;
+            U_(i, j) = 0;
         }
 
         if (g != 0)
@@ -226,7 +228,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
         {
             for (label j=i; j<Um; j++)
             {
-                U_(j, i) = 0.0;
+                U_(j, i) = 0;
             }
         }
 
@@ -235,32 +237,40 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
     for (label k=Un-1; k >= 0; k--)
     {
-        for (label its = 0; its < 35; its++)
+        for (label its = 0; its < 30; its++)
         {
             bool flag = true;
 
             label mn;
             for (l = k; l >= 0; l--)
             {
-                mn = l-1;
-                if (mag(rv1[l]) + anorm == anorm)
+                mn = l - 1;
+
+                if (l == 0 || mag(rv1[l]) <= anorm)
                 {
                     flag = false;
                     break;
                 }
-                if (mag(S_[mn]) + anorm == anorm) break;
+
+                if (mag(S_[mn]) <= anorm)
+                {
+                    break;
+                }
             }
 
             if (flag)
             {
-                scalar c = 0.0;
-                s = 1.0;
+                scalar c = 0;
+                s = 1;
                 for (label i=l; i<k+1; i++)
                 {
                     scalar f = s*rv1[i];
                     rv1[i] = c*rv1[i];
 
-                    if (mag(f) + anorm == anorm) break;
+                    if (mag(f) <= anorm)
+                    {
+                        break;
+                    }
 
                     g = S_[i];
                     scalar h = sqrtSumSqr(f, g);
@@ -283,18 +293,20 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
             if (l == k)
             {
-                if (z < 0.0)
+                if (z < 0)
                 {
                     S_[k] = -z;
-
-                    for (label j=0; j<Un; j++) V_(j, k) = -V_(j, k);
+                    for (label j=0; j<Un; j++)
+                    {
+                        V_(j, k) = -V_(j, k);
+                    }
                 }
                 break;
             }
-            if (its == 34)
+            if (its == 29)
             {
                 WarningInFunction
-                    << "No convergence in 35 SVD iterations"
+                    << "No convergence in 30 SVD iterations"
                     << endl;
             }
 
@@ -303,15 +315,15 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
             scalar y = S_[mn];
             g = rv1[mn];
             scalar h = rv1[k];
-            scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2.0*h*y);
+            scalar f = ((y - z)*(y + z) + (g - h)*(g + h))/(2*h*y);
             g = sqrtSumSqr(f, scalar(1));
             f = ((x - z)*(x + z) + h*((y/(f + sign(g, f))) - h))/x;
-            scalar c = 1.0;
-            s = 1.0;
+            scalar c = 1;
+            s = 1;
 
             for (label j=l; j <= mn; j++)
             {
-                label i=j + 1;
+                label i = j + 1;
                 g = rv1[i];
                 y = S_[i];
                 h = s*g;
@@ -352,7 +364,7 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
                     U_[jj][i] = z*c - y*s;
                 }
             }
-            rv1[l] = 0.0;
+            rv1[l] = 0;
             rv1[k] = f;
             S_[k] = x;
         }
@@ -364,7 +376,6 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
     {
         if (S_[i] <= minS)
         {
-            //Info<< "Removing " << S_[i] << " < " << minS << endl;
             S_[i] = 0;
             nZeros_++;
         }
@@ -372,28 +383,6 @@ Foam::SVD::SVD(const scalarRectangularMatrix& A, const scalar minCondition)
 
     // Now multiply out to find the pseudo inverse of A, VSinvUt_
     multiply(VSinvUt_, V_, inv(S_), U_.T());
-
-    // test SVD
-    /*
-    scalarRectangularMatrix SVDA(A.m(), A.n());
-    multiply(SVDA, U_, S_, transpose(V_));
-    scalar maxDiff = 0;
-    scalar diff = 0;
-    for (label i=0; i<A.m(); i++)
-    {
-        for (label j=0; j<A.n(); j++)
-        {
-            diff = mag(A(i, j) - SVDA(i, j));
-            if (diff > maxDiff) maxDiff = diff;
-        }
-    }
-    Info<< "Maximum discrepancy between A and svd(A) = " << maxDiff << endl;
-
-    if (maxDiff > 4)
-    {
-        Info<< "singular values " << S_ << endl;
-    }
-    */
 }
 
 
diff --git a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
index c3fdb9bb1..9ff7ade03 100644
--- a/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
+++ b/src/OpenFOAM/matrices/scalarMatrices/SVD/SVDI.H
@@ -28,7 +28,8 @@ License
 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 >= 0 ? a : -a) : (a >= 0 ? -a : a);
+    return b >= 0 ? a : -a;
 }
 
 
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
index 92d41405e..7ff29e3ad 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.C
@@ -24,7 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "chemPointISAT.H"
-#include <limits>
+#include "SVD.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -49,14 +49,14 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
 
     for (label k=0; k<nCols-1; k++)
     {
-        scale = 0.0;
+        scale = 0;
         for (label i=k; i<nCols; i++)
         {
-            scale=max(scale, fabs(R(i, k)));
+            scale=max(scale, mag(R(i, k)));
         }
-        if (scale == 0.0)
+        if (scale == 0)
         {
-            c[k] = d[k] = 0.0;
+            c[k] = d[k] = 0;
         }
         else
         {
@@ -64,7 +64,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
             {
                 R(i, k) /= scale;
             }
-            sum = 0.0;
+            sum = 0;
             for (label i=k; i<nCols; i++)
             {
                 sum += sqr(R(i, k));
@@ -75,7 +75,7 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
             d[k] = -scale*sigma;
             for (label j=k+1; j<nCols; j++)
             {
-                sum=0.0;
+                sum=0;
                 for ( label i=k; i<nCols; i++)
                 {
                     sum += R(i, k)*R(i, j);
@@ -88,14 +88,16 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrDecompose
             }
         }
     }
+
     d[nCols-1] = R(nCols-1, nCols-1);
+
     // form R
     for (label i=0; i<nCols; i++)
     {
         R(i, i) = d[i];
         for ( label j=0; j<i; j++)
         {
-            R(i, j)=0.0;
+            R(i, j)=0;
         }
     }
 }
@@ -110,42 +112,47 @@ void Foam::chemPointISAT<CompType, ThermoType>::qrUpdate
     const Foam::scalarField &v
 )
 {
-    label k, i;
+    label k;
+
     scalarField w(u);
-    for (k=n-1;k>=0;k--)
+    for (k=n-1; k>=0; k--)
     {
-        if (w[k] != 0.0)
+        if (w[k] != 0)
         {
             break;
         }
     }
+
     if (k < 0)
     {
-        k=0;
+        k = 0;
     }
-    for (i=k-1;i>=0;i--)
+
+    for (label i=k-1; i>=0; i--)
     {
         rotate(R, i, w[i],-w[i+1], n);
-        if (w[i] == 0.0)
+        if (w[i] == 0)
         {
-            w[i] = fabs(w[i+1]);
+            w[i] = mag(w[i+1]);
         }
-        else if (fabs(w[i]) > fabs(w[i+1]))
+        else if (mag(w[i]) > mag(w[i+1]))
         {
-            w[i] = fabs(w[i])*sqrt(1.0+sqr(w[i+1]/w[i]));
+            w[i] = mag(w[i])*sqrt(1 + sqr(w[i+1]/w[i]));
         }
         else
         {
-            w[i] = fabs(w[i+1])*sqrt(1.0+sqr(w[i]/w[i+1]));
+            w[i] = mag(w[i+1])*sqrt(1 + sqr(w[i]/w[i+1]));
         }
     }
-    for (i=0;i<n;i++)
+
+    for (label i=0; i<n; i++)
     {
         R(0, i) += w[0]*v[i];
     }
-    for (i=0;i<k;i++)
+
+    for (label i=0; i<k; i++)
     {
-        rotate(R, i, R(i, i),-R(i+1, i), n);
+        rotate(R, i, R(i, i), -R(i+1, i), n);
     }
 }
 
@@ -160,360 +167,36 @@ void Foam::chemPointISAT<CompType, ThermoType>::rotate
     label n
 )
 {
-    label j;
     scalar c, fact, s, w, y;
-    if (a == 0.0)
+    if (a == 0)
     {
-        c=0.0;
-        s=(b >= 0.0 ? 1.0 : -1.0);
+        c = 0;
+        s = (b >= 0 ? 1.0 : -1.0);
     }
-    else if (fabs(a) > fabs(b))
+    else if (mag(a) > mag(b))
     {
         fact = b/a;
-        c=sign(a)/sqrt(1.0+(fact*fact));
-        s=fact*c;
+        c = sign(a)/sqrt(1 + sqr(fact));
+        s = fact*c;
     }
     else
     {
-        fact=a/b;
-        s=sign(b)/sqrt(1.0+(fact*fact));
-        c=fact*s;
+        fact = a/b;
+        s = sign(b)/sqrt(1 + sqr(fact));
+        c = fact*s;
     }
-    for (j=i;j<n;j++)
+    for (label j=i;j<n;j++)
     {
-        y=R(i, j);
-        w=R(i+1, j);
-        R(i, j)=c*y-s*w;
-        R(i+1, j)=s*y+c*w;
-    }
-}
-
-
-template<class CompType, class ThermoType>
-void Foam::chemPointISAT<CompType, ThermoType>::svd
-(
-    scalarSquareMatrix& A,
-    label m,
-    label n,
-    scalarDiagonalMatrix& d,
-    scalarSquareMatrix& V
-)
-{
-    // UPDATED VERSION NR3
-    bool flag;
-    label i, its, j, jj, k, l, nm;
-    scalar anorm, c, f, g, h, s, scale, x, y, z;
-    scalarField rv1(n);
-    scalar eps = std::numeric_limits<scalar>::epsilon();
-    g = scale = anorm = 0.0;
-
-    // Householder reduction to bidiagonal form
-    for ( i = 0; i<n; i++)
-    {
-        l=i+2; // change from i+1 to i+2
-        rv1[i] = scale*g;
-        g=s=scale=0.0;
-        if (i < m)
-        {
-            for (k=i;k<m;k++)
-            {
-                scale += fabs(A(k, i));
-            }
-            if (scale != 0.0)
-            {
-                for ( k=i;k<m;k++)
-                {
-                    A(k, i) /= scale;
-                    s += A(k, i)*A(k, i);
-                }
-                f = A(i, i);
-                g = -sign(f)*sqrt(s);
-                h = f*g-s;
-                A(i, i)=f-g;
-                for (j=l-1;j<n;j++)
-                {
-                    for (s=0.0,k=i;k<m;k++)
-                    {
-                        s += A(k, i)*A(k, j);
-                    }
-                    f = s/h;
-                    for (k=i; k<m;k++)
-                    {
-                        A(k, j) += f*A(k, i);
-                    }
-                }
-                for (k=i; k<m;k++)
-                {
-                    A(k, i) *= scale;
-                }
-            }
-        }
-        d[i] = scale * g;
-        g=s=scale=0.0;
-
-        if (i+1 <= m && i+1 != n)
-        {
-            for (k=l-1; k<n; k++)
-            {
-                scale += fabs(A(i, k));
-            }
-            if (scale != 0.0)
-            {
-                for (k=l-1; k<n; k++)
-                {
-                    A(i, k) /= scale;
-                    s += A(i, k)*A(i, k);
-                }
-                f = A(i, l-1);
-                g = -sign(f)*sqrt(s);
-                h = f*g-s;
-                A(i, l-1) = f-g;
-                for (k=l-1; k<n; k++)
-                {
-                    rv1[k] = A(i, k)/h;
-                }
-                for (j=l-1; j<m; j++)
-                {
-                    for (s=0.0,k=l-1; k<n; k++)
-                    {
-                        s += A(j, k)*A(i, k);
-                    }
-                    for (k=l-1; k<n; k++)
-                    {
-                        A(j, k) += s*rv1[k];
-                    }
-                }
-                for (k=l-1; k<n; k++)
-                {
-                    A(i, k) *= scale;
-                }
-            }
-        }
-        anorm = max(anorm, (fabs(d[i])+fabs(rv1[i])));
-    }
-
-    // Accumulation of right-hand transformations
-    for (i=n-1; i>=0; i--)
-    {
-        if (i < n-1)
-        {
-            if (g != 0.0)
-            {
-                for (j=l; j<n; j++)
-                {
-                    V(j, i) = (A(i, j)/A(i, l))/g;
-                }
-                for (j=l; j<n; j++)
-                {
-                    for (s=0.0,k=l; k<n; k++)
-                    {
-                        s += A(i, k)*V(k, j);
-                    }
-                    for (k=l; k<n; k++)
-                    {
-                        V(k, j) += s*V(k, i);
-                    }
-                }
-            }
-            for (j=l; j<n; j++)
-            {
-                V(i, j)=V(j, i)=0.0;
-            }
-        }
-        V(i, i) = 1.0;
-        g = rv1[i];
-        l = i;
-    }
-    // Accumulation of left-hand transformations
-    for (i = min(m, n)-1; i>=0; i--)
-    {
-        l=i+1;
-        g=d[i];
-        for (j=l; j<n; j++)
-        {
-            A(i, j) = 0.0;
-        }
-        if (g != 0.0)
-        {
-            g = 1.0/g;
-            for (j=l; j<n; j++)
-            {
-                for (s=0.0, k=l; k<m; k++)
-                {
-                    s+= A(k, i)*A(k, j);
-                }
-                f = (s/A(i, i))*g;
-                for (k=i; k<m; k++)
-                {
-                    A(k, j) += f*A(k, i);
-                }
-            }
-            for (j=i; j<m; j++)
-            {
-                A(j, i) *= g;
-            }
-        }
-        else
-        {
-            for (j=i; j<m; j++)
-            {
-                A(j, i)=0.0;
-            }
-        }
-        ++A(i, i);
-    }
-
-    // Diagonalization of the bidiagonal form :
-    // Loop over singular values, and over allowed iteration
-    for (k=n-1; k>=0; k--)
-    {
-        for (its=0; its<30; its++)
-        {
-            flag=true;
-            // Test for splitting (rv1[1] always zero)
-            for (l=k; l>=0; l--)
-            {
-                nm = l-1;
-                if (l == 0 || fabs(rv1[l]) <= eps*anorm)
-                {
-                    flag = false;
-                    break;
-                }
-                if (fabs(d[nm]) <= eps*anorm)
-                {
-                    break;
-                }
-            }
-            // Cancellation of rv1[l], if l>1
-            if (flag)
-            {
-                c = 0.0;
-                s = 1.0;
-                for (i=l; i<k+1; i++)
-                {
-                    f = s*rv1[i];
-                    rv1[i] = c*rv1[i];
-                    if (fabs(f) <= eps*anorm)
-                    {
-                        break;
-                    }
-                    g = d[i];
-                    h = pythag(f, g);
-                    d[i] = h;
-                    h = 1.0/h;
-                    c = g*h;
-                    s = -f*h;
-                    for (j=0; j<m; j++)
-                    {
-                        y = A(j, nm);
-                        z = A(j, i);
-                        A(j, nm) = y*c + z*s;
-                        A(j, i) = z*c - y*s;
-                    }
-                }
-            }
-
-            z = d[k];
-            if (l == k) // Convergence
-            {
-                if (z < 0.0) // Singular value is made nonnegative
-                {
-                    d[k] = -z;
-                    for (j=0; j<n; j++)
-                    {
-                        V(j, k) = -V(j, k);
-                    }
-                }
-                break;
-            }
-            if (its == 34)
-            {
-                WarningInFunction
-                    << "no convergence in 35 SVD iterations"
-                    << endl;
-            }
-
-            x = d[l];
-            nm = k-1;
-            y = d[nm];
-            g = rv1[nm];
-            h = rv1[k];
-            f = ((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
-            g = pythag(f,1.0);
-            f = ((x-z)*(x+z)+h*((y/(f+sign(f)*g))-h))/x;
-            c=s=1.0;
-            // Next QR transformation
-            for (j=l; j<=nm; j++)
-            {
-                i = j+1;
-                g = rv1[i];
-                y = d[i];
-                h = s*g;
-                g = c*g;
-                z = pythag(f, h);
-                rv1[j] = z;
-                c = f/z;
-                s = h/z;
-                f = x*c + g*s;
-                g = g*c - x*s;
-                h = y*s;
-                y *= c;
-                for (jj=0; jj<n; jj++)
-                {
-                    x = V(jj, j);
-                    z = V(jj, i);
-                    V(jj, j) = x*c + z*s;
-                    V(jj, i) = z*c - x*s;
-                }
-                z = pythag(f, h);
-                d[j] = z;
-                if (z)
-                {
-                    z = 1.0/z;
-                    c = f*z;
-                    s = h*z;
-                }
-                f = c*g + s*y;
-                x = c*y - s*g;
-                for (jj=0; jj<m; jj++)
-                {
-                    y = A(jj, j);
-                    z = A(jj, i);
-                    A(jj, j) = y*c + z*s;
-                    A(jj, i) = z*c - y*s;
-                }
-            }
-            rv1[l] = 0.0;
-            rv1[k] = f;
-            d[k] = x;
-        }
-    }
-}
-
-
-// pythag function used in svd
-// compute (a^2+b^2)^1/2 without descrutive underflow or overflow
-template<class CompType, class ThermoType>
-Foam::scalar
-Foam::chemPointISAT<CompType, ThermoType>::pythag(scalar a, scalar b)
-{
-    scalar absa, absb;
-    absa = fabs(a);
-    absb = fabs(b);
-    if (absa > absb)
-    {
-        return absa*sqrt(1.0+sqr(absb/absa));
-    }
-    else
-    {
-        return (absb == 0.0 ? 0.0 : absb*sqrt(1.0+sqr(absa/absb)));
+        y = R(i, j);
+        w = R(i+1, j);
+        R(i, j) = c*y-s*w;
+        R(i+1, j) = s*y+c*w;
     }
 }
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-
 template<class CompType, class ThermoType>
 Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
 (
@@ -565,37 +248,38 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
     label reduOrCompDim = completeSpaceSize;
     if (isMechRedActive)
     {
-        reduOrCompDim = nActiveSpecies_+2;
+        reduOrCompDim = nActiveSpecies_ + 2;
     }
 
-    // SVD decomposition A= U*D*V^T
-    scalarSquareMatrix Atmp(A);// A computed in ISAT.C
-    scalarSquareMatrix B(reduOrCompDim, Zero);
-    DiagonalMatrix<scalar> diag(reduOrCompDim, Zero);
-    svd(Atmp, reduOrCompDim, reduOrCompDim, diag, B);
+    // SVD decomposition A = U*D*V^T
+    SVD svdA(A);
 
-    // replace the value of vector diag by max(diag, 1/2), first ISAT paper,
-    // Pope
+    scalarDiagonalMatrix D(reduOrCompDim);
+    const scalarDiagonalMatrix& S = svdA.S();
+
+    // Replace the value of vector D by max(D, 1/2), first ISAT paper
     for (label i=0; i<reduOrCompDim; i++)
     {
-        diag[i] = max(diag[i], 0.5);
+        D[i] = max(S[i], 0.5);
     }
 
-    // rebuild A with max length, tol and scale factor before QR decomposition
-    scalarSquareMatrix Atilde(reduOrCompDim, reduOrCompDim);
+    // Rebuild A with max length, tol and scale factor before QR decomposition
+    scalarRectangularMatrix Atilde(reduOrCompDim);
 
-    // result stored in Atilde
-    multiply(Atilde, Atmp, diag, B.T());
+    // Result stored in Atilde
+    multiply(Atilde, svdA.U(), D, svdA.V().T());
 
-    for (label i=0; i<reduOrCompDim; i++)// on species loop
+    for (label i=0; i<reduOrCompDim; i++)
     {
-        for (label j=0; j<reduOrCompDim; j++)// species, T and p loop
+        for (label j=0; j<reduOrCompDim; j++)
         {
-            label compi=i;
+            label compi = i;
+
             if (isMechRedActive)
             {
                 compi = simplifiedToCompleteIndex(i);
             }
+
             // SF*A/tolerance
             // (where SF is diagonal with inverse of scale factors)
             // SF*A is the same as dividing each line by the scale factor
@@ -639,6 +323,7 @@ Foam::chemPointISAT<CompType, ThermoType>::chemPointISAT
    tolerance_ = p.tolerance();
 }
 
+
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class CompType, class ThermoType>
@@ -647,27 +332,30 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
     scalarField dphi(phiq-phi());
     bool isMechRedActive = chemistry_.mechRed()->active();
     label dim = (isMechRedActive) ? nActiveSpecies_ : completeSpaceSize()-2;
-    scalar epsTemp=0.0;
-    List<scalar> propEps(completeSpaceSize(),0.0);
+    scalar epsTemp=0;
+    List<scalar> propEps(completeSpaceSize(),0);
 
     for (label i=0; i<completeSpaceSize()-2; i++)
     {
-        scalar temp(0.0);
+        scalar temp = 0;
+
         // When mechanism reduction is inactive OR on active species multiply L
         // by dphi to get the distance in the active species direction else (for
         // inactive species), just multiply the diagonal element and dphi
         if
         (
             !(isMechRedActive)
-          ||(isMechRedActive && completeToSimplifiedIndex_[i]!=-1)
+          ||(isMechRedActive && completeToSimplifiedIndex_[i] != -1)
         )
         {
             label si=(isMechRedActive) ? completeToSimplifiedIndex_[i] : i;
+
             for (label j=si; j<dim; j++)// LT is upper triangular
             {
                 label sj=(isMechRedActive) ? simplifiedToCompleteIndex_[j] : j;
                 temp += LT_(si, j)*dphi[sj];
             }
+
             temp += LT_(si, nActiveSpecies_)*dphi[completeSpaceSize()-2];
             temp += LT_(si, nActiveSpecies_+1)*dphi[completeSpaceSize()-1];
         }
@@ -675,12 +363,15 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
         {
             temp = dphi[i]/(tolerance_*scaleFactor_[i]);
         }
+
         epsTemp += sqr(temp);
+
         if (printProportion_)
         {
             propEps[i] = temp;
         }
     }
+
     // Temperature
     epsTemp +=
         sqr
@@ -688,6 +379,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
             LT_(dim, dim)*dphi[completeSpaceSize()-2]
            +LT_(dim, dim+1)*dphi[completeSpaceSize()-1]
         );
+
     // Pressure
     epsTemp += sqr(LT_(dim+1, dim+1)*dphi[completeSpaceSize()-1]);
 
@@ -703,12 +395,12 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
         propEps[completeSpaceSize()-1] =
             sqr(LT_(dim+1, dim+1)*dphi[completeSpaceSize()-1]);
     }
-    if (sqrt(epsTemp) > 1.0+tolerance_)
+    if (sqrt(epsTemp) > 1 + tolerance_)
     {
         if (printProportion_)
         {
-            scalar max=-1.0;
-            label maxIndex=-1;
+            scalar max = -1;
+            label maxIndex = -1;
             for (label i=0; i<completeSpaceSize(); i++)
             {
                 if(max < propEps[i])
@@ -736,7 +428,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::inEOA(const scalarField& phiq)
             Info<< "Direction maximum impact to error in ellipsoid: "
                 << propName << endl;
             Info<< "Proportion to the total error on the retrieve: "
-                << max / (epsTemp+SMALL) << endl;
+                << max/(epsTemp+SMALL) << endl;
         }
         return false;
     }
@@ -754,13 +446,13 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
     const scalarField& Rphiq
 )
 {
-    scalar eps2 = 0.0;
+    scalar eps2 = 0;
     scalarField dR(Rphiq - Rphi());
     scalarField dphi(phiq - phi());
     const scalarField& scaleFactorV(scaleFactor());
     const scalarSquareMatrix& Avar(A());
     bool isMechRedActive = chemistry_.mechRed()->active();
-    scalar dRl = 0.0;
+    scalar dRl = 0;
     label dim = completeSpaceSize()-2;
     if (isMechRedActive)
     {
@@ -771,12 +463,13 @@ bool Foam::chemPointISAT<CompType, ThermoType>::checkSolution
     // included
     for (label i=0; i<completeSpaceSize()-2; i++)
     {
-        dRl = 0.0;
+        dRl = 0;
         if (isMechRedActive)
         {
             label si = completeToSimplifiedIndex_[i];
+
             // If this species is active
-            if (si!=-1)
+            if (si != -1)
             {
                 for (label j=0; j<dim; j++)
                 {
@@ -846,7 +539,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
             // corresponds to an inactive on the query side
             if
             (
-                completeToSimplifiedIndex_[i]!=-1
+                completeToSimplifiedIndex_[i] != -1
              && chemistry_.completeToSimplifiedIndex()[i] == -1
             )
             {
@@ -860,7 +553,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
             (
                 completeToSimplifiedIndex_[i] == -1
              && chemistry_.completeToSimplifiedIndex()[i] == -1
-             && dphi[i] != 0.0
+             && dphi[i] != 0
             )
             {
                 activeAdded++;
@@ -933,17 +626,17 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
             {
                 LT_(i, i)=
                     1.0
-                  / (tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
-                A_(i, i)=1.0;
+                   /(tolerance_*scaleFactor_[simplifiedToCompleteIndex_[i]]);
+                A_(i, i) = 1;
             }
         }
 
-        dim = nActiveSpecies_+2;
+        dim = nActiveSpecies_ + 2;
     }
 
     // beginning of grow algorithm
-    scalarField phiTilde(dim, 0.0);
-    scalar normPhiTilde = 0.0;
+    scalarField phiTilde(dim, 0);
+    scalar normPhiTilde = 0;
     // p' = L^T.(p-phi)
 
     for (label i=0; i<dim; i++)
@@ -953,7 +646,7 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
             label sj = j;
             if (isMechRedActive)
             {
-                sj=simplifiedToCompleteIndex_[j];
+                sj = simplifiedToCompleteIndex_[j];
             }
             phiTilde[i] += LT_(i, j)*dphi[sj];
         }
@@ -961,12 +654,15 @@ bool Foam::chemPointISAT<CompType, ThermoType>::grow(const scalarField& phiq)
         phiTilde[i] += LT_(i, dim-1)*dphi[completeSpaceSize()-1];
         normPhiTilde += sqr(phiTilde[i]);
     }
+
     scalar invSqrNormPhiTilde = 1.0/normPhiTilde;
     normPhiTilde = sqrt(normPhiTilde);
+
     // gamma = (1/|p'| - 1)/|p'|^2
     scalar gamma = (1/normPhiTilde - 1)*invSqrNormPhiTilde;
     scalarField u(gamma*phiTilde);
-    scalarField v(dim,0.0);
+    scalarField v(dim, 0);
+
     for ( label i=0; i<dim; i++)
     {
         for (register label j=0; j<=i;j++)
diff --git a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
index 2345b537e..386174327 100644
--- a/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
+++ b/src/thermophysicalModels/chemistryModel/chemistryModel/TDACChemistryModel/tabulation/ISAT/chemPointISAT/chemPointISAT.H
@@ -208,8 +208,8 @@ class chemPointISAT
         (
             scalarSquareMatrix& R,
             const label n,
-            const scalarField &u,
-            const scalarField &v
+            const scalarField& u,
+            const scalarField& v
         );
 
         void rotate
@@ -221,29 +221,6 @@ class chemPointISAT
             label n
         );
 
-        //- Singular Value Decomposition (SVD) for a square matrix
-        //  needed to compute the the length of the hyperellipsoid semi-axes
-        //  SVD decompose a matrix A into:
-        //  A = U * D * V^T ,
-        //  with the singular value in the diagonal matrix D and U and V
-        //  orthogonal A (scalarMatrix) the square matrix to apply the
-        //  decomposition on m (label) the number of line of the matrix A
-        //  n (label) the size of the matrix A
-        //  Output:  U (scalarMatrix) replace A
-        //  V (scalarMatrix) not the transpose V^T
-        //  d (scalarField) the diagonal element of matrix D
-        void svd
-        (
-            scalarSquareMatrix& A,
-            label m,
-            label n,
-            scalarDiagonalMatrix& d,
-            scalarSquareMatrix& V
-        );
-
-        //- Function used in svd function
-        scalar pythag(scalar a, scalar b);
-
 
 public: