From 44e7dd91ad806b8c0329fb084f5ae442d3750229 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 21 May 2014 16:35:02 +0000 Subject: [PATCH] [Equil] Eliminate 'mlequ' and 'gaussj' in favor of using LAPACK --- include/cantera/equil/vcs_internal.h | 72 ------ src/equil/BasisOptimize.cpp | 108 +------- src/equil/vcs_elem.cpp | 17 +- src/equil/vcs_root1d.cpp | 17 +- src/equil/vcs_solve_TP.cpp | 40 +-- src/equil/vcs_util.cpp | 357 --------------------------- 6 files changed, 55 insertions(+), 556 deletions(-) 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;