From 626660102c0e92e511b9eb0ee00a10e03d015bbb Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Tue, 20 Sep 2011 23:41:17 +0000
Subject: [PATCH] Added a manual check for the satisfaction of Reaction
 Matrix's ability to conserve elements.

---
 Cantera/src/equil/vcs_solve_TP.cpp | 50 ++++++++++++++++++++++++++++++
 Cantera/src/equil/vcs_util.cpp     |  4 +--
 2 files changed, 52 insertions(+), 2 deletions(-)

diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp
index 1a169dd0f..59c4e9824 100644
--- a/Cantera/src/equil/vcs_solve_TP.cpp
+++ b/Cantera/src/equil/vcs_solve_TP.cpp
@@ -3623,9 +3623,59 @@ namespace VCSnonideal {
 	//plogf(" |  %6.2f", m_scSize[i]);
 	plogf("\n");
       }
+    
+
+      /*
+       *  Manual check on the satisfaction of the reaction matrix's ability
+       *  to conserve elements
+       */
+      double sum;
+      double sumMax = -1.0;
+      int iMax = -1;
+      int jMax = -1;
+      int n;
+      for (i = 0; i < m_numRxnTot; ++i) {
+	k = m_indexRxnToSpecies[i];
+	for (j = 0; j < ncTrial; ++j) {
+	  if (j == jlose) {
+	    sum = m_formulaMatrix[juse][k];
+	    for (n = 0; n < ncTrial; n++) {
+	      double numElements = m_formulaMatrix[juse][n];
+	      double coeff =  m_stoichCoeffRxnMatrix[i][n];
+	      sum += coeff * numElements;
+	    }
+	  } else {
+	    sum = m_formulaMatrix[j][k];
+	    for (n = 0; n < ncTrial; n++) {
+	      double numElements = m_formulaMatrix[j][n];
+	      double coeff =  m_stoichCoeffRxnMatrix[i][n];
+	      sum += coeff * numElements;
+	    }
+	  }
+	  if (fabs(sum) > sumMax) {
+	    sumMax = fabs(sum);
+	    iMax = i;
+	    jMax = j;
+	    if (j == jlose) {
+	      jMax = juse;
+	    }
+	  }
+	  if (fabs(sum) > 1.0E-6) {
+	    printf("we have a prob\n");
+	    exit(-1);
+	  }
+	}
+      }
+      plogf("   ---               largest error in Stoich coeff = %g at rxn = %d ", sumMax, iMax);
+      plogf("%-10.10s", m_speciesName[m_indexRxnToSpecies[iMax]].c_str());
+      plogf(" element = %d ", jMax);
+      plogf("%-5.5s", m_elementName[jMax].c_str());
+      plogf("\n");
       plogf("   "); for(i=0; i<77; i++) plogf("-"); plogf("\n");
     }
 #endif
+
+
     /* **************************************************** */
     /* **** EVALUATE DELTA N VALUES *********************** */
     /* **************************************************** */
diff --git a/Cantera/src/equil/vcs_util.cpp b/Cantera/src/equil/vcs_util.cpp
index 524ba39b1..61e79511e 100644
--- a/Cantera/src/equil/vcs_util.cpp
+++ b/Cantera/src/equil/vcs_util.cpp
@@ -550,7 +550,7 @@ namespace VCSnonideal {
     return 0;
   }
   //====================================================================================================================
-  //! Linear equation solution by Gauss-Jordan elimination for multiple rhs vectors
+  // Linear equation solution by Gauss-Jordan elimination for multiple rhs vectors
   /*
    * Solve a square matrix with multiple right hand sides
    *
@@ -558,7 +558,7 @@ namespace VCSnonideal {
    *     C X + B = 0;
    * \f]
    *
-   * This routine uses Gauss elimination and is optimized for the solution
+   * This routine uses Gauss-Jordan elimination with full pivoting and is optimized for the solution
    * of lots of rhs's. 
    *
    * @return Routine returns an integer representing success: