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

Replaced mlequ with gaussj full pivoting in a few spots.

Browse source 
There were no changes above numerical roundoff in the test suite.
This commit is contained in:
Harry Moffat 2011-09-13 18:32:08 +00:00
parent b1b6dd945c
commit c58697dbd3
3 changed files with 66 additions and 3 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
Cantera/src/equil/vcs_internal.h
View file

@ -158,6 +158,8 @@ namespace VCSnonideal {
*/
int vcsUtil_mlequ(double *c, int idem, int n, double *b, int m);
int vcsUtil_gaussj(double *c, int idem, int n, double *b, int m);
//! Swap values in vector of doubles
/*!
* Switches the value of x[i1] with x[i2]

10
Cantera/src/equil/vcs_solve_TP.cpp
View file

@ -3517,7 +3517,11 @@ namespace VCSnonideal {
* Use Gauss-Jordon block elimination to calculate
* the reaction matrix, m_stoichCoeffRxnMatrix[][].
*/
j = vcsUtil_mlequ(sm, m_numElemConstraints, ncTrial, m_stoichCoeffRxnMatrix[0], m_numRxnTot);
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;
@ -3566,7 +3570,9 @@ namespace VCSnonideal {
}
}
}
j = vcsUtil_mlequ(sm, m_numElemConstraints, ncTrial, aw, 1);
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;

57
Cantera/src/equil/vcs_util.cpp
View file

@ -551,13 +551,57 @@ namespace VCSnonideal {
}
//====================================================================================================================
//! Linear equation solution by Gauss-Jordan elimination for multiple rhs vectors
/*
* 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.
*
* @return Routine returns an integer representing success:
* - 1 : Matrix is singluar
* - 0 : solution is OK
* The solution x[] is returned in the matrix b.
*
* @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, int idem, int n, double *b, int m) {
static int vcsUtil_gaussj(double *c, int idem, int n, double *b, int m) {
int i, j, k, l, ll;
int irow = -1;
int icol = -1;
bool needInverse = false;
double pivinv, dum;
#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<int> indxc(n);
std::vector<int> indxr(n);
std::vector<int> ipiv(n, 0);
@ -620,6 +664,17 @@ namespace VCSnonideal {
}
}
}
/*
* 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;
}
//====================================================================================================================