C++ Array2D::nRows方法代码示例

ChebyshevRate::ChebyshevRate(double Tmin, double Tmax, double Pmin, double Pmax,
                             const Array2D& coeffs)
    : Tmin_(Tmin)
    , Tmax_(Tmax)
    , Pmin_(Pmin)
    , Pmax_(Pmax)
    , nP_(coeffs.nColumns())
    , nT_(coeffs.nRows())
    , chebCoeffs_(coeffs.nColumns() * coeffs.nRows(), 0.0)
    , dotProd_(coeffs.nRows())
{
    double logPmin = std::log10(Pmin);
    double logPmax = std::log10(Pmax);
    double TminInv = 1.0 / Tmin;
    double TmaxInv = 1.0 / Tmax;

    TrNum_ = - TminInv - TmaxInv;
    TrDen_ = 1.0 / (TmaxInv - TminInv);
    PrNum_ = - logPmin - logPmax;
    PrDen_ = 1.0 / (logPmax - logPmin);

    for (size_t t = 0; t < nT_; t++) {
        for (size_t p = 0; p < nP_; p++) {
            chebCoeffs_[nP_*t + p] = coeffs(t,p);
        }
    }
}

NMatrix make_G(const NMatrix &P, const Array2D<SS_GTO> &orbitals)
{
    const int norbitals = orbitals.nRows();
    
    NMatrix G(norbitals, norbitals);
    
    for (int i=0; i<norbitals; ++i)
    {
        for (int j=0; j<norbitals; ++j)
        {
            double v = 0;
            
            for (int k=0; k<norbitals; ++k)
            {
                for (int l=0; l<norbitals; ++l)
                {
                    v += P(k,l)*
                            ( electron_integral(orbitals(i,j), orbitals(k,l))
                              - 0.5*electron_integral(orbitals(i,l),orbitals(k,j)) );
                }
            }
            
            G(i,j) = v;
        }
    }
    
    return G;
}

static NMatrix make_overlap_matrix(const Array2D<SS_GTO> &orbitals)
{
    const int norbitals = orbitals.nRows();
    
    NMatrix overlap_matrix(norbitals, norbitals);

    //calculate the overlap matrix - this is the overlap integral
    //of all pairs of orbitals
    for (int i=0; i<norbitals; ++i)
    {
        for (int j=0; j<norbitals; ++j)
        {
            overlap_matrix(i,j) = overlap_integral( orbitals(i,j) );
        }
    }

    return overlap_matrix;
}

  /*
   * @param s          output stream 
   * @param title      plot title 
   * @param names      vector of variable names
   * @param data       N x M data array.
   *                        data(n,m) is the m^th value of the n^th variable.
   */
  void outputExcel(std::ostream &s, const std::string &title, 
		   const std::vector<std::string>& names,  
		   const Array2D& data) {
    int i,j;
    int npts = static_cast<int>(data.nColumns());
    int nv = static_cast<int>(data.nRows());
    s << title + "," << endl;
    for (i = 0; i < nv; i++) {
      s << names[i] << ",";
    }
    s << endl;
    for (i = 0; i < npts; i++) {
      for (j = 0; j < nv; j++) {
	s << data(j,i) << ",";
      }
      s << endl;
    }
  }

  /*
   * @param s        output stream 
   * @param title    plot title 
   * @param names    vector of variable names 
   * @param data      N x M data array. 
   *                 data(n,m) is the m^th value of the n^th variable.
   */
  void outputTEC(std::ostream &s, const std::string &title,
		 const std::vector<std::string>& names,  
		 const Array2D& data) {
    int i,j;
    int npts = static_cast<int>(data.nColumns());
    int nv = static_cast<int>(data.nRows());
    s << "TITLE     = \"" + title + "\"" << endl;
    s << "VARIABLES = " << endl;
    for (i = 0; i < nv; i++) {
      s << "\"" << names[i] << "\"" << endl;
    }
    s << "ZONE T=\"zone1\"" << endl;
    s << " I=" << npts << ",J=1,K=1,F=POINT" << endl;
    s << "DT=( ";
    for (i = 0; i < nv; i++) s << " SINGLE";
    s << " )" << endl;
    for (i = 0; i < npts; i++) {
      for (j = 0; j < nv; j++) {
	s << data(j,i) << " ";
      }
      s << endl;
    }
  }

static NMatrix make_core_fock_matrix(const Array2D<SS_GTO> &orbitals,
                                     const QVector<PointCharge> &charges)
{
    const int norbitals = orbitals.nRows();
    const int ncharges = charges.count();

    NMatrix fock_matrix(norbitals, norbitals);
    NMatrix kinetic_matrix(norbitals, norbitals);
    NMatrix nuclear_matrix(norbitals, norbitals, 0);
    
    //kinetic energy
    for (int i=0; i<norbitals; ++i)
    {
        for (int j=0; j<norbitals; ++j)
        {
            kinetic_matrix(i,j) = kinetic_integral( orbitals(i,j) ); 
        }
    }
    
    qDebug() << "KINETIC MATRIX\n" << kinetic_matrix.toString();
    
    //nuclear-electron energy
    for (int k=0; k<ncharges; ++k)
    {
        NMatrix my_nuclear_matrix(norbitals, norbitals);
    
        for (int i=0; i<norbitals; ++i)
        {
            for (int j=0; j<norbitals; ++j)
            {
                my_nuclear_matrix(i,j) = potential_integral(charges[k], orbitals(i,j));
            }
        }

        nuclear_matrix += my_nuclear_matrix;
        
        qDebug() << "NUCLEAR MATRIX" << k << "\n" << my_nuclear_matrix.toString();
    }

    qDebug() << "NUCLEAR MATRIX\n" << nuclear_matrix.toString();
    
    for (int i=0; i<norbitals; ++i)
    {
        for (int j=0; j<norbitals; ++j)
        {
            fock_matrix(i,j) = kinetic_matrix(i,j) + nuclear_matrix(i,j);
        }
    }
    
    for (int i=0; i<norbitals; ++i)
    {
        for (int j=0; j<norbitals; ++j)
        {
            for (int k=0; k<norbitals; ++k)
            {
                for (int l=0; l<norbitals; ++l)
                {
                    qDebug() << i+1 << j+1 << k+1 << l+1 
                             << electron_integral(orbitals(i,j), orbitals(k,l));
                }
            }
        }
    }
    
    return fock_matrix;
}

void getMatrixValues(const XML_Node& node,
                     const std::vector<std::string>& keyStringRow,
                     const std::vector<std::string>& keyStringCol,
                     Array2D& retnValues, const bool convert,
                     const bool matrixSymmetric)
{
    if (keyStringRow.size() > retnValues.nRows()) {
        throw CanteraError("getMatrixValues",
                           "size of key1 greater than numrows");
    } else if (keyStringCol.size() > retnValues.nColumns()) {
        throw CanteraError("getMatrixValues",
                           "size of key2 greater than num cols");
    } else if (matrixSymmetric && retnValues.nRows() != retnValues.nColumns()) {
        throw CanteraError("getMatrixValues",
                           "nrow != ncol for a symmetric matrix");
    }

    /*
     * Get the attributes field, units, from the XML node
     * and determine the conversion factor, funit.
     */
    doublereal funit = 1.0;
    if (convert && node["units"] != "") {
        funit = toSI(node["units"]);
    }

    vector<string> v;
    getStringArray(node, v);
    for (size_t i = 0; i < v.size(); i++) {
        size_t icolon = v[i].find(":");
        if (icolon == string::npos) {
            throw CanteraError("getMatrixValues","Missing two colons ("
                               +v[i]+")");
        }
        string key1 = v[i].substr(0,icolon);
        string rmm = v[i].substr(icolon+1, v[i].size());

        icolon = rmm.find(":");
        if (icolon == string::npos) {
            throw CanteraError("getMatrixValues","Missing one colon ("
                               +v[i]+")");
        }

        size_t irow = find(keyStringRow.begin(), keyStringRow.end(), key1)
                      - keyStringRow.begin();
        if (irow == keyStringRow.size()) {
            throw CanteraError("getMatrixValues","Row not matched by string: "
                               + key1);
        }

        string key2 = rmm.substr(0,icolon);
        size_t icol = find(keyStringCol.begin(), keyStringCol.end(), key2)
                      - keyStringCol.begin();
        if (icol == keyStringCol.size()) {
            throw CanteraError("getMatrixValues","Col not matched by string: "
                               + key2);
        }
        double dval = fpValueCheck(rmm.substr(icolon+1, rmm.size())) * funit;
        /*
         * Finally, insert the value;
         */
        retnValues(irow, icol) = dval;
        if (matrixSymmetric) {
            retnValues(icol, irow) = dval;
        }
    }
}