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

[Equil] Eliminate direct calls to LAPACK functions

Browse source
This commit is contained in:
Ray Speth 2016-04-10 23:57:03 -04:00
parent 961ca06cd5
commit a7d854f3b9
3 changed files with 24 additions and 48 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

28
src/equil/BasisOptimize.cpp
View file

@ -3,7 +3,6 @@
* stoichiometric coefficient matrix (see /ref equil functions)
*/
#include "cantera/equil/MultiPhase.h"
#include "cantera/numerics/ctlapack.h"
using namespace std;
@ -80,7 +79,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
mphase->getMoles(molNum.data());
// Other workspace
vector_fp sm(ne*ne, 0.0);
DenseMatrix sm(ne, ne);
vector_fp ss(ne, 0.0);
vector_fp sa(ne, 0.0);
if (formRxnMatrix.size() < nspecies*ne) {
@ -136,7 +135,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
size_t jl = jr;
for (j = 0; j < ne; ++j) {
size_t jj = orderVectorElements[j];
sm[j + jr*ne] = mphase->nAtoms(kk,jj);
sm(j, jr) = mphase->nAtoms(kk,jj);
}
if (jl > 0) {
// Compute the coefficients of JA column of the the upper
@ -145,7 +144,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
for (j = 0; j < jl; ++j) {
ss[j] = 0.0;
for (size_t i = 0; i < ne; ++i) {
ss[j] += sm[i + jr*ne] * sm[i + j*ne];
ss[j] += sm(i, jr) * sm(i, j);
}
ss[j] /= sa[j];
}
@ -154,7 +153,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// columns
for (j = 0; j < jl; ++j) {
for (size_t i = 0; i < ne; ++i) {
sm[i + jr*ne] -= ss[j] * sm[i + j*ne];
sm(i, jr) -= ss[j] * sm(i, j);
}
}
}
@ -163,8 +162,7 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// It will be used in the denominator in future row calcs.
sa[jr] = 0.0;
for (size_t ml = 0; ml < ne; ++ml) {
double tmp = sm[ml + jr*ne];
sa[jr] += tmp * tmp;
sa[jr] += pow(sm(ml, jr), 2);
}
// IF NORM OF NEW ROW .LT. 1E-3 REJECT
@ -224,11 +222,13 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
// Note the rearrangement of elements need only be done once in the problem.
// It's actually very similar to the top of this program with ne being the
// species and nc being the elements!!
sm.resize(nComponents, nComponents);
for (size_t k = 0; k < nComponents; ++k) {
size_t kk = orderVectorSpecies[k];
for (size_t j = 0; j < nComponents; ++j) {
size_t jj = orderVectorElements[j];
sm[j + k*ne] = mphase->nAtoms(kk, jj);
sm(j, k) = mphase->nAtoms(kk, jj);
}
}
@ -240,15 +240,9 @@ size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn, MultiPhase* mphase,
formRxnMatrix[j + i * ne] = - mphase->nAtoms(kk, jj);
}
}
// 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("BasisOptimize", "factorization returned an error condition");
}
ct_dgetrs(ctlapack::NoTranspose, nComponents, nNonComponents, &sm[0], ne,
&ipiv[0], &formRxnMatrix[0], ne, info);
// // Use LU factorization to calculate the reaction matrix
solve(sm, formRxnMatrix.data(), nNonComponents, ne);
if (BasisOptimize_print_lvl >= 1) {
writelog(" ---\n");

17
src/equil/vcs_elem.cpp
View file

@ -6,7 +6,7 @@
*/
#include "cantera/equil/vcs_solve.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/numerics/ctlapack.h"
#include "cantera/numerics/DenseMatrix.h"
namespace Cantera
{
@ -223,25 +223,18 @@ int VCS_SOLVE::vcs_elcorr(double aa[], double x[])
// Ok, do the general case. Linear algebra problem is of length nc, not ne,
// as there may be degenerate rows when nc .ne. ne.
DenseMatrix A(m_numComponents, m_numComponents);
for (size_t i = 0; i < m_numComponents; ++i) {
x[i] = m_elemAbundances[i] - m_elemAbundancesGoal[i];
if (fabs(x[i]) > 1.0E-13) {
retn = 1;
}
for (size_t j = 0; j < m_numComponents; ++j) {
aa[j + i*m_numElemConstraints] = - m_formulaMatrix(i,j);
A(j, i) = - m_formulaMatrix(i,j);
}
}
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);
solve(A, x, 1, m_numElemConstraints);
// Now apply the new direction without creating negative species.
double par = 0.5;

27
src/equil/vcs_solve_TP.cpp
View file

@ -13,7 +13,7 @@
#include "cantera/base/ctexceptions.h"
#include "cantera/base/clockWC.h"
#include "cantera/base/stringUtils.h"
#include "cantera/numerics/ctlapack.h"
#include "cantera/numerics/DenseMatrix.h"
#include <cstdio>
@ -2128,6 +2128,7 @@ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[],
size_t k;
size_t juse = npos;
size_t jlose = npos;
DenseMatrix C;
clockWC tickTock;
if (m_debug_print_lvl >= 2) {
plogf(" ");
@ -2418,9 +2419,10 @@ L_END_LOOP:
// the rearrangement of elements need only be done once in the problem. It's
// actually very similar to the top of this program with ne being the
// species and nc being the elements!!
C.resize(ncTrial, ncTrial);
for (size_t j = 0; j < ncTrial; ++j) {
for (size_t i = 0; i < ncTrial; ++i) {
sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,i);
C(i, j) = m_formulaMatrix(j,i);
}
}
for (size_t i = 0; i < m_numRxnTot; ++i) {
@ -2431,14 +2433,7 @@ L_END_LOOP:
}
// Solve the linear system to calculate the reaction matrix,
// m_stoichCoeffRxnMatrix.
int info;
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.ptrColumn(0), m_numElemConstraints, info);
solve(C, m_stoichCoeffRxnMatrix.ptrColumn(0), m_numRxnTot, m_numElemConstraints);
// NOW, if we have interfacial voltage unknowns, what we did was just wrong
// -> hopefully it didn't blow up. Redo the problem. Search for inactive E
@ -2459,9 +2454,9 @@ L_END_LOOP:
for (size_t j = 0; j < ncTrial; ++j) {
for (size_t i = 0; i < ncTrial; ++i) {
if (i == jlose) {
sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,juse);
C(i, j) = m_formulaMatrix(j,juse);
} else {
sm[i + j*m_numElemConstraints] = m_formulaMatrix(j,i);
C(i, j) = m_formulaMatrix(j,i);
}
}
}
@ -2476,13 +2471,7 @@ L_END_LOOP:
}
}
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);
solve(C, aw, 1, m_numElemConstraints);
size_t i = k - ncTrial;
for (size_t j = 0; j < ncTrial; j++) {
m_stoichCoeffRxnMatrix(j,i) = aw[j];