diff --git a/src/equil/BasisOptimize.cpp b/src/equil/BasisOptimize.cpp index a91028926..0e90e0f78 100644 --- a/src/equil/BasisOptimize.cpp +++ b/src/equil/BasisOptimize.cpp @@ -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"); diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 9a3bfaa63..e090b84a5 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -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; diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index ff13af2de..a5c55511e 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -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 @@ -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];