Method vcs_basopt

src/equil/vcs_solve_TP.cpp:2015–2497

}

int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[],
                          double ss[], double test, bool* const usedZeroedSpecies)
{
    size_t k;
    size_t juse = npos;
    size_t jlose = npos;
    DenseMatrix C;
    clockWC tickTock;
    if (m_debug_print_lvl >= 2) {
        plogf("   ");
        for (size_t i=0; i<77; i++) {
            plogf("-");
        }
        plogf("\n");
        plogf("   --- Subroutine BASOPT called to ");
        if (doJustComponents) {
            plogf("calculate the number of components\n");
        } else {
            plogf("reevaluate the components\n");
        }
        if (m_debug_print_lvl >= 2) {
            plogf("\n");
            plogf("   ---       Formula Matrix used in BASOPT calculation\n");
            plogf("   ---      Active |   ");
            for (size_t j = 0; j < m_nelem; j++) {
                plogf("  %1d  ", m_elementActive[j]);
            }
            plogf("\n");
            plogf("   ---     Species | ");
            for (size_t j = 0; j < m_nelem; j++) {
                writelog(" {:>8.8s}", m_elementName[j]);
            }
            plogf("\n");
            for (k = 0; k < m_nsp; k++) {
                writelog("   --- {:>11.11s} | ", m_speciesName[k]);
                for (size_t j = 0; j < m_nelem; j++) {
                    plogf(" %8.2g", m_formulaMatrix(k,j));
                }
                plogf("\n");
            }
            writelogendl();
        }
    }

    // Calculate the maximum value of the number of components possible. It's
    // equal to the minimum of the number of elements and the number of total
    // species.
    size_t ncTrial = std::min(m_nelem, m_nsp);
    m_numComponents = ncTrial;
    *usedZeroedSpecies = false;
    vector<int> ipiv(ncTrial);

    // Use a temporary work array for the mole numbers, aw[]
    std::copy(m_molNumSpecies_old.begin(),
              m_molNumSpecies_old.begin() + m_nsp, aw);

    // Take out the Voltage unknowns from consideration
    for (k = 0; k < m_nsp; k++) {