research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Equil] Eliminate 'mlequ' and 'gaussj' in favor of using LAPACK

Browse source
This commit is contained in:
Ray Speth 2014-05-21 16:35:02 +00:00
parent ef781b44e6
commit 44e7dd91ad
6 changed files with 55 additions and 556 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

72
include/cantera/equil/vcs_internal.h
View file

@ -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 <I>B</I>.
* 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 <I>B</I>.
* 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.

108
src/equil/BasisOptimize.cpp
View file

@ -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<size_t>& orderVectorSpecies,
std::vector<size_t>& 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<size_t>& orderVectorSpecies,

17
src/equil/vcs_elem.cpp
View file

@ -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.
*/

17
src/equil/vcs_root1d.cpp
View file

@ -9,9 +9,12 @@
*/
#include "cantera/equil/vcs_internal.h"
#include "cantera/numerics/ctlapack.h"
#include <cstdio>
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]);

40
src/equil/vcs_solve_TP.cpp
View file

@ -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 <cstdio>
@ -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];

357
src/equil/vcs_util.cpp
View file

@ -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<int> 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<size_t> indxc(n);
std::vector<size_t> indxr(n);
std::vector<int> 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;