C++ SparseMatrix::innerIndexPtr方法代码示例

ColCompressedMatrix convert_from_Eigen(const Eigen::SparseMatrix<double> &m)
{
	assert(m.rows() == m.cols());
	assert(m.isCompressed());
	return ColCompressedMatrix(m.rows, m.nonZeros(), 
		m.valuePtr(), m.outerIndexPtr(), m.innerIndexPtr());
}

    /// Solve the linear system Ax = b, with A being the
    /// combined derivative matrix of the residual and b
    /// being the residual itself.
    /// \param[in] residual   residual object containing A and b.
    /// \return               the solution x
    NewtonIterationBlackoilSimple::SolutionVector
    NewtonIterationBlackoilSimple::computeNewtonIncrement(const LinearisedBlackoilResidual& residual) const
    {
        typedef LinearisedBlackoilResidual::ADB ADB;
        const int np = residual.material_balance_eq.size();
        ADB mass_res = residual.material_balance_eq[0];
        for (int phase = 1; phase < np; ++phase) {
            mass_res = vertcat(mass_res, residual.material_balance_eq[phase]);
        }
        const ADB well_res = vertcat(residual.well_flux_eq, residual.well_eq);
        const ADB total_residual = collapseJacs(vertcat(mass_res, well_res));

        Eigen::SparseMatrix<double, Eigen::RowMajor> matr;
        total_residual.derivative()[0].toSparse(matr);

        SolutionVector dx(SolutionVector::Zero(total_residual.size()));
        Opm::LinearSolverInterface::LinearSolverReport rep
            = linsolver_->solve(matr.rows(), matr.nonZeros(),
                                matr.outerIndexPtr(), matr.innerIndexPtr(), matr.valuePtr(),
                                total_residual.value().data(), dx.data(), parallelInformation_);

        // store iterations
        iterations_ = rep.iterations;

        if (!rep.converged) {
            OPM_THROW(LinearSolverProblem,
                      "FullyImplicitBlackoilSolver::solveJacobianSystem(): "
                      "Linear solver convergence failure.");
        }
        return dx;
    }

ConvertToMklResult to_mkl(const Eigen::SparseMatrix<double> &Ain,
                          sparse_status_t &status) {
  ConvertToMklResult result;

  const int N = static_cast<int>(Ain.rows());
  // const-cast to work with C-api.
  int *row_starts = const_cast<int *>(Ain.outerIndexPtr());
  int *col_index = const_cast<int *>(Ain.innerIndexPtr());
  double *values = const_cast<double *>(Ain.valuePtr());

  result.descr.type = SPARSE_MATRIX_TYPE_GENERAL; /* Full matrix is stored */

  result.status =
      mkl_sparse_d_create_csr(&result.matrix, SPARSE_INDEX_BASE_ZERO, N, N,
                              row_starts, row_starts + 1, col_index, values);
  return result;
}

    void NewtonIterationBlackoilInterleaved::formInterleavedSystem(const std::vector<ADB>& eqs,
                                                                   const Eigen::SparseMatrix<double, Eigen::RowMajor>& A,
                                                                   Mat& istlA) const
    {
        const int np = eqs.size();

        // Find sparsity structure as union of basic block sparsity structures,
        // corresponding to the jacobians with respect to pressure.
        // Use addition to get to the union structure.
        Eigen::SparseMatrix<double> structure = eqs[0].derivative()[0];
        for (int phase = 0; phase < np; ++phase) {
            structure += eqs[phase].derivative()[0];
        }
        Eigen::SparseMatrix<double, Eigen::RowMajor> s = structure;

        // Create ISTL matrix with interleaved rows and columns (block structured).
        assert(np == 3);
        istlA.setSize(s.rows(), s.cols(), s.nonZeros());
        istlA.setBuildMode(Mat::row_wise);
        const int* ia = s.outerIndexPtr();
        const int* ja = s.innerIndexPtr();
        for (Mat::CreateIterator row = istlA.createbegin(); row != istlA.createend(); ++row) {
            int ri = row.index();
            for (int i = ia[ri]; i < ia[ri + 1]; ++i) {
                row.insert(ja[i]);
            }
        }
        const int size = s.rows();
        Span span[3] = { Span(size, 1, 0),
                         Span(size, 1, size),
                         Span(size, 1, 2*size) };
        for (int row = 0; row < size; ++row) {
            for (int col_ix = ia[row]; col_ix < ia[row + 1]; ++col_ix) {
                const int col = ja[col_ix];
                MatrixBlockType block;
                for (int p1 = 0; p1 < np; ++p1) {
                    for (int p2 = 0; p2 < np; ++p2) {
                        block[p1][p2] = A.coeff(span[p1][row], span[p2][col]);
                    }
                }
                istlA[row][col] = block;
            }
        }
    }