diff --git a/include/cantera/equil/vcs_internal.h b/include/cantera/equil/vcs_internal.h
index d8ea21c24..2f248d010 100644
--- a/include/cantera/equil/vcs_internal.h
+++ b/include/cantera/equil/vcs_internal.h
@@ -106,78 +106,6 @@ public:
*/
double vcsUtil_gasConstant(int mu_units);
-//! Invert an n x n matrix and solve m rhs's
-/*!
- * Solve a square matrix with multiple right hand sides
- *
- * \f[
- * C X + B = 0;
- * \f]
- *
- * This routine uses Gauss elimination and is optimized for the solution
- * of lots of rhs's. A crude form of row pivoting is used here.
- * The matrix C is destroyed during the solve.
- *
- * @return The solution x[] is returned in the matrix B.
- * Routine returns an integer representing success:
- * - 1 : Matrix is singular
- * - 0 : solution is OK
- *
- * @param c Matrix to be inverted. c is in fortran format, i.e., rows
- * are the inner loop. Row numbers equal to idem.
- * c[i+j*idem] = c_i_j = Matrix to be inverted:
- * - i = row number
- * - j = column number
- *
- * @param idem number of row dimensions in c
- * @param n Number of rows and columns in c
- * @param b Multiple RHS. Note, b is actually the negative of
- * most formulations. Row numbers equal to idem.
- * b[i+j*idem] = b_i_j = vectors of rhs's:
- * - i = row number
- * - j = column number
- * (each column is a new rhs)
- * @param m number of rhs's
- */
-int vcsUtil_mlequ(double* c, size_t idem, size_t n, double* b, size_t m);
-
-//! Invert an n x n matrix and solve m rhs's
-/*!
- * Solve a square matrix with multiple right hand sides
- *
- * \f[
- * C X + B = 0;
- * \f]
- *
- * This routine uses Gauss-Jordan elimination and is optimized for the solution
- * of lots of rhs's. Full row and column pivoting is used here. It's been
- * shown to be necessary in at least one case.
- * The matrix C is destroyed during the solve.
- *
- * @return The solution x[] is returned in the matrix B.
- * Routine returns an integer representing success:
- * - 1 : Matrix is singular
- * - 0 : solution is OK
- *
- * @param c Matrix to be inverted. c is in fortran format, i.e., rows
- * are the inner loop. Row numbers equal to idem.
- * c[i+j*idem] = c_i_j = Matrix to be inverted:
- * - i = row number
- * - j = column number
- *
- * @param idem number of row dimensions in c
- * @param n Number of rows and columns in c
- * @param b Multiple RHS. Note, b is actually the negative of
- * most formulations. Row numbers equal to idem.
- * b[i+j*idem] = b_i_j = vectors of rhs's:
- * - i = row number
- * - j = column number
- * (each column is a new rhs)
- * @param m number of rhs's
- */
-int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m);
-
-
//! Definition of the function pointer for the root finder
/*!
* see vcsUtil_root1d for a definition of how to use this.
diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp
index b6261facf..bd162d879 100644
--- a/src/equil/BasisOptimize.cpp
+++ b/src/equil/BasisOptimize.cpp
@@ -4,6 +4,7 @@
*/
#include "cantera/thermo/ThermoPhase.h"
#include "cantera/equil/MultiPhase.h"
+#include "cantera/numerics/ctlapack.h"
using namespace Cantera;
using namespace std;
@@ -37,31 +38,6 @@ static void print_stringTrunc(const char* str, int space, int alignment);
*/
static size_t amax(double* x, size_t j, size_t n);
-//! Invert an nxn matrix and solve m rhs's
-/*!
- * Solve C X + B = 0
- *
- * This routine uses Gauss elimination and is optimized for the solution of
- * lots of rhs's. A crude form of row pivoting is used here.
- *
- * @param c C is the matrix to be inverted
- * @param idem first dimension in the calling routine.
- * idem >= n must be true
- * @param n number of rows and columns in the matrix
- * @param b rhs of the matrix problem
- * @param m number of rhs to be solved for
- *
- * - c[i+j*idem] = c_i_j = Matrix to be inverted
- * - b[i+j*idem] = b_i_j = vectors of rhs's. Each column is a new rhs.
- *
- * Where j = column number and i = row number.
- *
- * @return Retuns 1 if the matrix is singular, or 0 if the solution is OK
- *
- * The solution is returned in the matrix b.
- */
-static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m);
-
size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
MultiPhase* mphase, std::vector& orderVectorSpecies,
std::vector& orderVectorElements,
@@ -357,18 +333,18 @@ size_t Cantera::BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
kk = orderVectorSpecies[k];
for (j = 0; j < nComponents; ++j) {
jj = orderVectorElements[j];
- formRxnMatrix[j + i * ne] = mphase->nAtoms(kk, jj);
+ formRxnMatrix[j + i * ne] = - mphase->nAtoms(kk, jj);
}
}
- /*
- * Use Gauss-Jordan block elimination to calculate
- * the reaction matrix
- */
- int ierr = mlequ(DATA_PTR(sm), ne, nComponents, DATA_PTR(formRxnMatrix), nNonComponents);
- if (ierr == 1) {
- writelog("ERROR: mlequ returned an error condition\n");
- throw CanteraError("basopt", "mlequ returned an error condition");
+ // Use LU factorization to calculate the reaction matrix
+ int info;
+ vector_int ipiv(nComponents);
+ ct_dgetrf(nComponents, nComponents, &sm[0], ne, &ipiv[0], info);
+ if (info) {
+ throw CanteraError("basopt", "factorization returned an error condition");
}
+ ct_dgetrs(ctlapack::NoTranspose, nComponents, nNonComponents, &sm[0], ne,
+ &ipiv[0], &formRxnMatrix[0], ne, info);
#ifdef DEBUG_MODE
if (Cantera::BasisOptimize_print_lvl >= 1) {
@@ -490,70 +466,6 @@ static size_t amax(double* x, size_t j, size_t n)
return largest;
}
-static int mlequ(double* c, size_t idem, size_t n, double* b, size_t m)
-{
- size_t i, j, k, l;
- double R;
-
- /*
- * Loop over the rows
- * -> At the end of each loop, the only nonzero entry in the column
- * will be on the diagonal. We can therfore just invert the
- * diagonal at the end of the program to solve the equation system.
- */
- for (i = 0; i < n; ++i) {
- if (c[i + i * idem] == 0.0) {
- /*
- * Do a simple form of row pivoting to find a non-zero pivot
- */
- bool foundPivot = false;
- for (k = i + 1; k < n; ++k) {
- if (c[k + i * idem] != 0.0) {
- foundPivot = true;
- break;
- }
- }
-
- if (!foundPivot) {
-#ifdef DEBUG_MODE
- writelogf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i);
-#endif
- return 1;
- }
-
- for (j = 0; j < n; ++j) {
- c[i + j * idem] += c[k + j * idem];
- }
- for (j = 0; j < m; ++j) {
- b[i + j * idem] += b[k + j * idem];
- }
- }
-
- for (l = 0; l < n; ++l) {
- if (l != i && c[l + i * idem] != 0.0) {
- R = c[l + i * idem] / c[i + i * idem];
- c[l + i * idem] = 0.0;
- for (j = i+1; j < n; ++j) {
- c[l + j * idem] -= c[i + j * idem] * R;
- }
- for (j = 0; j < m; ++j) {
- b[l + j * idem] -= b[i + j * idem] * R;
- }
- }
- }
- }
- /*
- * The negative in the last expression is due to the form of B upon
- * input
- */
- for (i = 0; i < n; ++i) {
- for (j = 0; j < m; ++j) {
- b[i + j * idem] = -b[i + j * idem] / c[i + i*idem];
- }
- }
- return 0;
-}
-
size_t Cantera::ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
MultiPhase* mphase,
std::vector& orderVectorSpecies,
diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp
index e2e2522be..961454b5a 100644
--- a/src/equil/vcs_elem.cpp
+++ b/src/equil/vcs_elem.cpp
@@ -4,11 +4,15 @@
* element abundances constraints and the algorithm for fixing violations
* of the element abundances constraints.
*/
+#include "cantera/base/ct_defs.h"
#include "cantera/equil/vcs_solve.h"
#include "cantera/equil/vcs_internal.h"
#include "cantera/base/ctexceptions.h"
+#include "cantera/numerics/ctlapack.h"
#include "math.h"
+using namespace Cantera;
+
namespace VCSnonideal
{
void VCS_SOLVE::vcs_elab()
@@ -266,14 +270,19 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
retn = 1;
}
for (size_t j = 0; j < m_numComponents; ++j) {
- aa[j + i*m_numElemConstraints] = m_formulaMatrix[j][i];
+ aa[j + i*m_numElemConstraints] = - m_formulaMatrix[j][i];
}
}
- int err = vcsUtil_mlequ(aa, m_numElemConstraints, m_numComponents, x, 1);
- if (err == 1) {
- plogf("vcs_elcorr ERROR: mlequ returned error condition\n");
+ int info;
+ vector_int ipiv(std::min(m_numComponents, m_numElemConstraints));
+ ct_dgetrf(m_numComponents, m_numComponents, aa, m_numElemConstraints,
+ &ipiv[0], info);
+ if (info) {
+ plogf("vcs_elcorr ERROR: matrix factorization\n");
return VCS_FAILED_CONVERGENCE;
}
+ ct_dgetrs(ctlapack::NoTranspose, m_numComponents, 1, aa,
+ m_numElemConstraints, &ipiv[0], x, m_numElemConstraints, info);
/*
* Now apply the new direction without creating negative species.
*/
diff --git a/src/equil/vcs_root1d.cpp b/src/equil/vcs_root1d.cpp
index 6917f9c5f..6fe2025f8 100644
--- a/src/equil/vcs_root1d.cpp
+++ b/src/equil/vcs_root1d.cpp
@@ -9,9 +9,12 @@
*/
#include "cantera/equil/vcs_internal.h"
+#include "cantera/numerics/ctlapack.h"
#include
+using namespace Cantera;
+
namespace VCSnonideal
{
@@ -136,6 +139,8 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax,
posStraddle = false;
}
}
+ int ipiv[3];
+ int info;
do {
/*
@@ -171,13 +176,15 @@ int vcsUtil_root1d(double xmin, double xmax, size_t itmax,
c[6] = SQUARE(x0);
c[7] = SQUARE(x1);
c[8] = SQUARE(x2);
- f[0] = - f0;
- f[1] = - f1;
- f[2] = - f2;
- retn = vcsUtil_mlequ(c, 3, 3, f, 1);
- if (retn == 1) {
+ f[0] = f0;
+ f[1] = f1;
+ f[2] = f2;
+
+ ct_dgetrf(3, 3, c, 3, ipiv, info);
+ if (info) {
goto QUAD_BAIL;
}
+ ct_dgetrs(ctlapack::NoTranspose, 3, 1, c, 3, ipiv, f, 3, info);
root = f[1]* f[1] - 4.0 * f[0] * f[2];
if (root >= 0.0) {
xn1 = (- f[1] + sqrt(root)) / (2.0 * f[2]);
diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp
index 11464be07..1eeefc776 100644
--- a/src/equil/vcs_solve_TP.cpp
+++ b/src/equil/vcs_solve_TP.cpp
@@ -16,6 +16,7 @@
#include "cantera/base/ctexceptions.h"
#include "cantera/base/clockWC.h"
#include "cantera/base/stringUtils.h"
+#include "cantera/numerics/ctlapack.h"
#include
@@ -2942,6 +2943,8 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
ncTrial = std::min(m_numElemConstraints, m_numSpeciesTot);
m_numComponents = ncTrial;
*usedZeroedSpecies = false;
+ vector_int ipiv(ncTrial);
+ int info;
/*
* Use a temporary work array for the mole numbers, aw[]
@@ -3242,22 +3245,18 @@ L_END_LOOP:
for (i = 0; i < m_numRxnTot; ++i) {
k = m_indexRxnToSpecies[i];
for (j = 0; j < ncTrial; ++j) {
- m_stoichCoeffRxnMatrix[i][j] = m_formulaMatrix[j][k];
+ m_stoichCoeffRxnMatrix[i][j] = - m_formulaMatrix[j][k];
}
}
- /*
- * Use Gauss-Jordan block elimination to calculate
- * the reaction matrix, m_stoichCoeffRxnMatrix[][].
- */
-
- j = vcsUtil_gaussj(sm, m_numElemConstraints, ncTrial, m_stoichCoeffRxnMatrix[0], m_numRxnTot);
- // j = vcsUtil_mlequ(sm, m_numElemConstraints, ncTrial, m_stoichCoeffRxnMatrix[0], m_numRxnTot);
-
-
- if (j == 1) {
- plogf("vcs_solve_TP ERROR: mlequ returned an error condition\n");
- return VCS_FAILED_CONVERGENCE;
+ // Solve the linear system to calculate the reaction matrix,
+ // m_stoichCoeffRxnMatrix[][].
+ ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info);
+ if (info) {
+ plogf("vcs_solve_TP ERROR: Error factorizing stoichiometric coefficient matrix\n");
+ return VCS_FAILED_CONVERGENCE;
}
+ ct_dgetrs(ctlapack::NoTranspose, ncTrial, m_numRxnTot, sm, m_numElemConstraints,
+ &ipiv[0], m_stoichCoeffRxnMatrix[0], m_numElemConstraints, info);
/*
* NOW, if we have interfacial voltage unknowns, what we did
@@ -3296,19 +3295,20 @@ L_END_LOOP:
k = m_indexRxnToSpecies[i];
for (j = 0; j < ncTrial; ++j) {
if (j == jlose) {
- aw[j] = m_formulaMatrix[juse][k];
+ aw[j] = - m_formulaMatrix[juse][k];
} else {
- aw[j] = m_formulaMatrix[j][k];
+ aw[j] = - m_formulaMatrix[j][k];
}
}
}
- j = vcsUtil_gaussj(sm, m_numElemConstraints, ncTrial, aw, 1);
- // j = vcsUtil_mlequ(sm, m_numElemConstraints, ncTrial, aw, 1);
- if (j == 1) {
- plogf("vcs_solve_TP ERROR: mlequ returned an error condition\n");
- return VCS_FAILED_CONVERGENCE;
+ ct_dgetrf(ncTrial, ncTrial, sm, m_numElemConstraints, &ipiv[0], info);
+ if (info) {
+ plogf("vcs_solve_TP ERROR: Error factorizing matrix\n");
+ return VCS_FAILED_CONVERGENCE;
}
+ ct_dgetrs(ctlapack::NoTranspose, ncTrial, 1, sm, m_numElemConstraints,
+ &ipiv[0], aw, m_numElemConstraints, info);
i = k - ncTrial;
for (j = 0; j < ncTrial; j++) {
m_stoichCoeffRxnMatrix[i][j] = aw[j];
diff --git a/src/equil/vcs_util.cpp b/src/equil/vcs_util.cpp
index bf5a0bf1c..5c58ac840 100644
--- a/src/equil/vcs_util.cpp
+++ b/src/equil/vcs_util.cpp
@@ -146,363 +146,6 @@ int vcs_max_int(const int* vector, int length)
return retn;
}
-#ifdef DEBUG_HKM
-static void mlequ_matrixDump(double* c, int idem, int n)
-{
- int i, j;
- printf("vcsUtil_mlequ() MATRIX DUMP --------------------------------------------------\n");
- printf(" ");
- for (j = 0; j < n; ++j) {
- printf(" % 3d ", j);
- }
- printf("\n");
- for (j = 0; j < n; ++j) {
- printf("-----------");
- }
- printf("\n");
- for (i = 0; i < n; ++i) {
- printf(" %3d | ", i);
- for (j = 0; j < n; ++j) {
- printf("% 10.3e ", c[i + j * idem]);
- }
- printf("\n");
- }
- for (j = 0; j < n; ++j) {
- printf("-----------");
- }
- printf("\n");
- printf("vcsUtil_mlequ() END MATRIX DUMP --------------------------------------------------\n");
-
-}
-#endif
-
-//! Swap rows in the c matrix and the b rhs matrix
-/*!
- * @param c Matrix of size nxn, row first
- * @param idem C storage dimension for the number of rows
- * @param n Size of the matrix
- * @param b RHS of the Ax=b problem to solve
- * @param m Number of rhs to solve
- * @param irowa first row to swap
- * @param irowb second row to swap
- */
-static void vcsUtil_swapRows(double* c, size_t idem, size_t n, double* b,
- size_t m, size_t irowa, size_t irowb)
-{
- if (irowa == irowb) {
- return;
- }
- for (size_t j = 0; j < n; j++) {
- std::swap(c[irowa + j * idem], c[irowb + j * idem]);
- }
- for (size_t j = 0; j < m; j++) {
- std::swap(b[irowa + j * idem], b[irowb + j * idem]);
- }
-}
-
-//! Swap rows in the c matrix and the b rhs matrix to lower the condition number of the matrix
-/*!
- * @param c Matrix of size nxn, row first
- * @param idem C storage dimension for the number of rows
- * @param n Size of the matrix
- * @param b RHS of the Ax=b problem to solve
- * @param m Number of rhs to solve
- */
-static void vcsUtil_mlequ_preprocess(double* c, size_t idem, size_t n,
- double* b, size_t m)
-{
- size_t j = 0;
- std::vector irowUsed(n, 0);
-
- for (j = 0; j < n; j++) {
- int numNonzero = 0;
- size_t inonzero = npos;
- for (size_t i = 0; i < n; i++) {
- if (c[i + j * idem] != 0.0) {
- numNonzero++;
- inonzero = i;
- }
- }
- if (numNonzero == 1) {
- if (inonzero != j) {
- if (irowUsed[inonzero] == 0) {
- vcsUtil_swapRows(c, idem, n, b, m, j, inonzero);
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
-#endif
- }
- }
- irowUsed[j] = 1;
- }
- }
-
- for (j = 0; j < n; j++) {
- if (c[j + j * idem] == 0.0) {
- int numNonzero = 0;
- size_t inonzero = npos;
- for (size_t i = 0; i < n; i++) {
- if (!irowUsed[i]) {
- if (c[i + j * idem] != 0.0) {
- if ((c[i + i * idem] == 0.0)
- || (c[j + i * idem] != 0.0)) {
- numNonzero++;
- inonzero = i;
- }
- }
- }
- }
- if (numNonzero == 1) {
- if (inonzero != j) {
- if (irowUsed[inonzero] == 0) {
- vcsUtil_swapRows(c, idem, n, b, m, j, inonzero);
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
-#endif
- }
- }
- irowUsed[j] = 1;
- }
- }
- }
-
- for (j = 0; j < n; j++) {
- if (c[j + j * idem] == 0.0) {
- int numNonzero = 0;
- size_t inonzero = npos;
- for (size_t i = 0; i < n; i++) {
- if (!irowUsed[i]) {
- if (c[i + j * idem] != 0.0) {
- if ((c[i + i * idem] == 0.0)
- || (c[j + i * idem] != 0.0)) {
- numNonzero++;
- inonzero = i;
- }
- }
- }
- }
- if (inonzero != npos) {
- if (inonzero != j) {
- if (irowUsed[inonzero] == 0) {
- vcsUtil_swapRows(c, idem, n, b, m, j, inonzero);
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
-#endif
- }
- }
- }
- }
- }
-}
-
-int vcsUtil_mlequ(double* c, size_t idem, size_t n, double* b, size_t m)
-{
- size_t k;
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
-#endif
- vcsUtil_mlequ_preprocess(c, idem, n, b, m);
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
- static int s_numCalls = 0;
- s_numCalls++;
-#endif
-
- double R;
- if (n > idem || n <= 0) {
- plogf("vcsUtil_mlequ ERROR: badly dimensioned matrix: %d %d\n", n, idem);
- return 1;
- }
-
-#ifdef DEBUG_HKM
- int dmatrix = 0;
- for (size_t i = 0; i < n; ++i) {
- bool notFound = true;
- for (size_t j = 0; j < n; ++j) {
- if (c[i + j * idem] != 0.0) {
- notFound = false;
- }
- }
- if (notFound) {
- printf(" vcsUtil_mlequ ERROR(): row %d is identically zero\n", i);
- }
- }
- for (size_t j = 0; j < n; ++j) {
- bool notFound = true;
- for (size_t i = 0; i < n; ++i) {
- if (c[i + j * idem] != 0.0) {
- notFound = false;
- }
- }
- if (notFound) {
- printf(" vcsUtil_mlequ ERROR(): column %d is identically zero\n", j);
- }
- }
- // if (s_numCalls >= 32) {
- // printf("vcsUtil_mlequ: we are here\n");
- // dmatrix = 1;
- // }
-
- if (dmatrix) {
- mlequ_matrixDump(c, idem, n);
- }
-#endif
- /*
- * Loop over the rows
- * -> At the end of each loop, the only nonzero entry in the column
- * will be on the diagonal. We can therfore just invert the
- * diagonal at the end of the program to solve the equation system.
- */
- for (size_t i = 0; i < n; ++i) {
- if (c[i + i * idem] == 0.0) {
- /*
- * Do a simple form of row pivoting to find a non-zero pivot
- */
- for (k = i + 1; k < n; ++k) {
- if (c[k + i * idem] != 0.0) {
- goto FOUND_PIVOT;
- }
- }
- plogf("vcsUtil_mlequ ERROR: Encountered a zero column: %d\n", i);
-#ifdef DEBUG_HKM
- plogf(" call # %d\n", s_numCalls);
-#endif
-#ifdef DEBUG_HKM
- mlequ_matrixDump(c, idem, n);
-#endif
- return 1;
-FOUND_PIVOT:
- ;
- for (size_t j = 0; j < n; ++j) {
- c[i + j * idem] += c[k + j * idem];
- }
- for (size_t j = 0; j < m; ++j) {
- b[i + j * idem] += b[k + j * idem];
- }
- }
-
- for (size_t l = 0; l < n; ++l) {
- if (l != i && c[l + i * idem] != 0.0) {
- R = c[l + i * idem] / c[i + i * idem];
- c[l + i * idem] = 0.0;
- for (size_t j = i + 1; j < n; ++j) {
- c[l + j * idem] -= c[i + j * idem] * R;
- }
- for (size_t j = 0; j < m; ++j) {
- b[l + j * idem] -= b[i + j * idem] * R;
- }
- }
- }
- }
- /*
- * The negative in the last expression is due to the form of B upon
- * input
- */
- for (size_t i = 0; i < n; ++i) {
- for (size_t j = 0; j < m; ++j) {
- b[i + j * idem] = -b[i + j * idem] / c[i + i * idem];
- }
- }
- return 0;
-}
-
-int vcsUtil_gaussj(double* c, size_t idem, size_t n, double* b, size_t m)
-{
- size_t i, j, k, l, ll;
- size_t irow = npos;
- size_t icol = npos;
- bool needInverse = false;
- double pivinv;
-#ifdef DEBUG_HKM
- static int s_numCalls = 0;
- s_numCalls++;
-#endif
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
-#endif
- /*
- * Preprocess the problem
- */
- vcsUtil_mlequ_preprocess(c, idem, n, b, m);
-
-#ifdef DEBUG_HKM
- // mlequ_matrixDump(c, idem, n);
-#endif
-
- std::vector indxc(n);
- std::vector indxr(n);
- std::vector ipiv(n, 0);
- doublereal big = 0.0;
- /*
- * This is the main loop over the columns to be reduced.
- */
- for (i = 0; i < n; i++) {
- big = 0.0;
- for (j = 0; j < n; j++) {
- if (ipiv[j] != 1) {
- for (k = 0; k < n; k++) {
- if (ipiv[k] == 0) {
- if (fabs(c[j + idem * k]) >= big) {
- big = fabs(c[j + idem * k]);
- irow = j;
- icol = k;
- }
- }
- }
- }
- }
- ++(ipiv[icol]);
- if (irow != icol) {
- vcsUtil_swapRows(c, idem, n, b, m, irow, icol);
- }
- indxr[i] = irow;
- indxc[i] = icol;
- if (c[icol + idem * icol] == 0.0) {
- plogf("vcsUtil_gaussj ERROR: Encountered a zero column: %d\n", i);
- return 1;
- }
- pivinv = 1.0 / c[icol + idem * icol];
- c[icol + idem * icol] = 1.0;
- for (l = 0; l < n; l++) {
- c[icol + idem * l] *= pivinv;
- }
- for (l = 0; l < m; l++) {
- b[icol + idem * l] *= pivinv;
- }
- for (ll = 0; ll < n; ll++) {
- if (ll != icol) {
- double dum = c[ll + idem * icol];
- c[ll + idem * icol] = 0;
- for (l = 0; l < n; l++) {
- c[ll + idem * l] -= c[icol + idem * l] * dum;
- }
- for (l = 0; l < m; l++) {
- b[ll + idem * l] -= b[icol + idem * l] * dum;
- }
- }
- }
- }
- if (needInverse) {
- for (l = n - 1; l != npos; l--) {
- if (indxr[l] != indxc[l]) {
- for (k = 0; k < n; k++) {
- std::swap(c[k + idem * indxr[l]], c[k + idem * indxr[l]]);
- }
- }
- }
- }
-
- /*
- * The negative in the last expression is due to the form of B upon
- * input
- */
- for (i = 0; i < n; ++i) {
- for (j = 0; j < m; ++j) {
- b[i + j * idem] = -b[i + j * idem];
- }
- }
- return 0;
-}
-
double vcsUtil_gasConstant(int mu_units)
{
double r;