Merge pull request #52 from gemshub/fix_LU

Fix lu

Author: gdmiron

CMakeLists.txt

@@ -2,7 +2,7 @@
 cmake_minimum_required(VERSION 3.16..4.2)
 
 # Set the name of the project
-project(GEMS3K VERSION 4.4.5 LANGUAGES CXX C)
+project(GEMS3K VERSION 4.5.0 LANGUAGES CXX C)
 
 set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CXX_STANDARD_REQUIRED ON)

GEMS3K/jama_cholesky.h

@@ -36,7 +36,7 @@ using namespace TNT;
 		L = chol.getL();
 
   	else
-        c out << "factorization was not complete.\n";
+		cout << "factorization was not complete.\n";
 
 	</pre>
 
@@ -115,8 +115,7 @@ Cholesky<Real>::Cholesky(const Array2D<Real> &A)
      for (int j = 0; j < n; j++)
 	 {
         double d = 0.0;
-        int k;
-        for ( k = 0; k < j; k++)
+        for (int k = 0; k < j; k++)
 		{
             Real s = 0.0;
             for (int i = 0; i < k; i++)
@@ -132,8 +131,8 @@ Cholesky<Real>::Cholesky(const Array2D<Real> &A)
 //         if( !isspd ) // Added SD 29/11/2006
 //             return;
 //        L_[j][j] = sqrt(d > 0.0 ? d : 0.0);
- L_[j][j] = sqrt(d > 0.0 ? d : 1e-60);  //  Test change!  DK 13.10.2006
-         for ( k = j+1; k < n; k++)
+            L_[j][j] = sqrt(d > 0.0 ? d : 1e-60);  //  Test change!  DK 13.10.2006
+         for (int k = j+1; k < n; k++)
 		 {
             L_[j][k] = 0.0;
          }
@@ -153,7 +152,7 @@ Cholesky<Real>::Cholesky(const Array2D<Real> &A)
 template <class Real>
 Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
 {
-	int k, n = L_.dim1();
+	int n = L_.dim1();
 	if (b.dim1() != n)
 		return Array1D<Real>();
 
@@ -162,7 +161,7 @@ Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
 
 
       // Solve L*y = b;
-      for ( k = 0; k < n; k++)
+      for (int k = 0; k < n; k++) 
 	  {
          for (int i = 0; i < k; i++)
                x[k] -= x[i]*L_[k][i];
@@ -171,7 +170,7 @@ Array1D<Real> Cholesky<Real>::solve(const Array1D<Real> &b)
       }
 
       // Solve L'*X = Y;
-      for ( k = n-1; k >= 0; k--)
+      for (int k = n-1; k >= 0; k--) 
 	  {
          for (int i = k+1; i < n; i++)
                x[k] -= x[i]*L_[i][k];

GEMS3K/jama_lu.h

@@ -3,7 +3,6 @@
 
 #include "tnt.h"
 #include <algorithm>
-#include "v_detail.h"
 using namespace TNT;
 using namespace std;
 
@@ -82,16 +81,16 @@ class LU
       // Outer loop.
 
       for (int j = 0; j < n; j++) {
-         int i;
+
          // Make a copy of the j-th column to localize references.
 
-         for ( i = 0; i < m; i++) {
+         for (int i = 0; i < m; i++) {
             LUcolj[i] = LU_[i][j];
          }
 
          // Apply previous transformations.
 
-         for ( i = 0; i < m; i++) {
+         for (int i = 0; i < m; i++) {
             LUrowi = LU_[i];
 
             // Most of the time is spent in the following dot product.
@@ -108,8 +107,8 @@ class LU
          // Find pivot and exchange if necessary.
 
          int p = j;
-         for ( i = j+1; i < m; i++) {
-            if (fabs(LUcolj[i]) > fabs(LUcolj[p])) {
+         for (int i = j+1; i < m; i++) {
+            if (abs(LUcolj[i]) > abs(LUcolj[p])) {
                p = i;
             }
          }
@@ -128,7 +127,7 @@ class LU
 
          // Compute multipliers.
 
-         if ( (j < m) && noZero(LU_[j][j]) ) {
+         if ( (j < m) && (LU_[j][j] != 0.0) ) {
             for (int ii = j+1; ii < m; ii++) {
                LU_[ii][j] /= LU_[j][j];
             }
@@ -144,7 +143,7 @@ class LU
 
    int isNonsingular () {
       for (int j = 0; j < n; j++) {
-         if ( approximatelyZero( LU_[j][j] ))
+         if (LU_[j][j] == 0)
             return 0;
       }
       return 1;
@@ -270,7 +269,7 @@ class LU
 
    Array1D<Real> solve (const Array1D<Real> &b)
    {
-     int k;
+
 	  /* Dimensions: A is mxn, X is nxk, B is mxk */
 
       if (b.dim1() != m) {
@@ -284,14 +283,14 @@ class LU
       Array1D<Real> x = permute_copy(b, piv1);
 
       // Solve L*Y = B(piv)
-      for ( k = 0; k < n; k++) {
+      for (int k = 0; k < n; k++) {
          for (int i = k+1; i < n; i++) {
                x[i] -= x[k]*LU_[i][k];
             }
          }
 
 	  // Solve U*X = Y;
-      for ( k = n-1; k >= 0; k--) {
+      for (int k = n-1; k >= 0; k--) {
             x[k] /= LU_[k][k];
       		for (int i = 0; i < k; i++)
             	x[i] -= x[k]*LU_[i][k];

GEMS3K/s_solmod3.cpp

@@ -28,7 +28,7 @@
 //-------------------------------------------------------------------
 
 #include "s_solmod.h"
-//#include "v_detail.h"
+#include "v_detail.h"
 
 //=============================================================================================
 // Implementation of ideal mixing model for multicomponent solid solutions
@@ -2433,7 +2433,7 @@ long int TBerman::ExcessPart()
             if( y0jsm <= 0. )
                continue; // skip - this moiety is not present on s site in this end member
 
-            if (y[s][m] == 0. && y0jsm > 0.) // DM 08.03.2024 // to prevent zerodivide in eq (5.2-3)
+            if (approximatelyZero(y[s][m]) && y0jsm > 0.) // DM 08.03.2024 // to prevent zerodivide in eq (5.2-3)
                 continue;
 
             // looking through the parameters list

GEMS3K/s_solmod5.cpp

@@ -31,6 +31,7 @@ using namespace std;
 //#include "verror.h"
 #include "s_solmod.h"
 #include "jsonconfig.h"
+#include "v_detail.h"
 
 //=============================================================================================
 // SIT model (NEA version) reimplementation for aqueous electrolyte solutions
@@ -1573,7 +1574,7 @@ void TPitzer::ETHETAS(double ZJ, double ZK, double I, double DH_term, double& et
    //*etheta = 0.0;
    //*ethetap = 0.0;
 
-   if (ZJ == ZK)
+    if (essentiallyEqual(ZJ, ZK))
       return /*(OK)*/;
 
    const double XCON = 6.0e0 * DH_term * sqrt(I);