20#ifndef OPM_EXTRACT_MATRIX_HEADER_INCLUDED
21#define OPM_EXTRACT_MATRIX_HEADER_INCLUDED
34 template <
class Matrix>
35 void copySubMatrix(
const Matrix& A,
const std::vector<int>& indices, Matrix& B)
38 for (
auto row = B.begin(); row != B.end(); ++row) {
39 for (
auto col = row->begin(); col != row->end(); ++col) {
40 *col = A[indices[row.index()]][indices[col.index()]];
46 template <
class Matrix>
49 assert(std::is_sorted(indices.begin(), indices.end()));
52 const std::size_t n = indices.size();
53 std::size_t newrow = 0;
54 enum { NotIncluded = -1 };
55 std::vector<int> index_map(m.N(), NotIncluded);
56 for (
auto row = m.begin(); row != m.end(); ++row) {
57 const int row_index = row.index();
58 if (row_index == indices[newrow]) {
59 index_map[row_index] = newrow;
70 for (
auto row = m.begin(); row != m.end(); ++row) {
71 if (index_map[row.index()] != NotIncluded) {
72 for (
auto col = row->begin(); col != row->end(); ++col) {
73 if (index_map[col.index()] != NotIncluded) {
81 Matrix res(n, n, nnz, Matrix::row_wise);
82 auto from_row = m.begin();
83 for (
auto row = res.createbegin(); row != res.createend(); ++row) {
85 while (
static_cast<int>(from_row.index()) < indices[row.index()]) {
88 assert(
static_cast<int>(from_row.index()) == indices[row.index()]);
90 for (
auto from_col = from_row->begin(); from_col != from_row->end(); ++from_col) {
91 const int new_col = index_map[from_col.index()];
92 if (new_col != NotIncluded) {
103 template <
class Vector>
106 Vector res(indices.size());
107 for (std::size_t ii = 0; ii < indices.size(); ++ii) {
108 res[ii] = x[indices[ii]];
114 template <
class Vector>
115 void setGlobal(
const Vector& x,
const std::vector<int>& indices, Vector& global_x)
117 for (std::size_t ii = 0; ii < indices.size(); ++ii) {
118 global_x[indices[ii]] = x[ii];
124 template <
class Matrix>
128 if (m1.N() != m2.N())
return false;
129 if (m1.M() != m2.M())
return false;
130 if (m1.nonzeroes() != m2.nonzeroes())
return false;
132 auto row1 = m1.begin();
133 auto row2 = m2.begin();
134 for (; row1 != m1.end(); ++row1, ++row2) {
135 if (row2 == m2.end())
return false;
136 if (row1.index() != row2.index())
return false;
137 auto col1 = row1->begin();
138 auto col2 = row2->begin();
139 for (; col1 != row1->end(); ++col1, ++col2) {
140 if (col2 == row2->end())
return false;
141 if (col1.index() != col2.index())
return false;
142 if (*col1 != *col2)
return false;
Vector extractVector(const Vector &x, const std::vector< int > &indices)
Definition: extractMatrix.hpp:104
void copySubMatrix(const Matrix &A, const std::vector< int > &indices, Matrix &B)
Definition: extractMatrix.hpp:35
void setGlobal(const Vector &x, const std::vector< int > &indices, Vector &global_x)
Definition: extractMatrix.hpp:115
bool matrixEqual(const Matrix &m1, const Matrix &m2)
Definition: extractMatrix.hpp:125
Matrix extractMatrix(const Matrix &m, const std::vector< int > &indices)
Definition: extractMatrix.hpp:47
Definition: blackoilboundaryratevector.hh:37