16 #ifndef dealii_sparse_matrix_ez_h 17 # define dealii_sparse_matrix_ez_h 33 template <
typename number>
35 template <
typename number>
103 template <
typename number>
170 static const unsigned short invalid_diagonal =
171 static_cast<unsigned short>(-1);
193 const unsigned short index);
245 const unsigned short index);
321 const unsigned int default_increment = 1);
356 unsigned int default_increment = 1,
415 template <
class StreamType>
432 std::vector<size_type> &used_by_line,
433 const bool compute_by_line)
const;
460 const bool elide_zero_values =
true);
489 template <
typename number2>
491 add(
const std::vector<size_type> &indices,
493 const bool elide_zero_values =
true);
500 template <
typename number2>
502 add(
const std::vector<size_type> &row_indices,
503 const std::vector<size_type> &col_indices,
505 const bool elide_zero_values =
true);
516 template <
typename number2>
519 const std::vector<size_type> &col_indices,
520 const std::vector<number2> &
values,
521 const bool elide_zero_values =
true);
532 template <
typename number2>
538 const bool elide_zero_values =
true,
539 const bool col_indices_are_sorted =
false);
562 template <
typename MatrixType>
564 copy_from(
const MatrixType &source,
const bool elide_zero_values =
true);
573 template <
typename MatrixType>
575 add(
const number factor,
const MatrixType &
matrix);
608 template <
typename somenumber>
610 vmult(Vector<somenumber> &dst,
const Vector<somenumber> &src)
const;
617 template <
typename somenumber>
619 Tvmult(Vector<somenumber> &dst,
const Vector<somenumber> &src)
const;
625 template <
typename somenumber>
627 vmult_add(Vector<somenumber> &dst,
const Vector<somenumber> &src)
const;
634 template <
typename somenumber>
636 Tvmult_add(Vector<somenumber> &dst,
const Vector<somenumber> &src)
const;
657 template <
typename somenumber>
660 const Vector<somenumber> &src,
661 const number omega = 1.)
const;
666 template <
typename somenumber>
669 const Vector<somenumber> & src,
670 const number om = 1.,
671 const std::vector<std::size_t> &pos_right_of_diagonal =
672 std::vector<std::size_t>())
const;
678 template <
typename somenumber>
681 const Vector<somenumber> &src,
682 const number om = 1.)
const;
688 template <
typename somenumber>
691 const Vector<somenumber> &src,
692 const number om = 1.)
const;
702 template <
typename MatrixTypeA,
typename MatrixTypeB>
705 const MatrixTypeB &B,
747 print(std::ostream &out)
const;
771 const unsigned int precision = 3,
772 const bool scientific =
true,
773 const unsigned int width = 0,
774 const char * zero_string =
" ",
775 const double denominator = 1.)
const;
815 <<
"The entry with index (" << arg1 <<
',' << arg2
816 <<
") does not exist.");
821 <<
"An entry with index (" << arg1 <<
',' << arg2
822 <<
") cannot be allocated.");
850 template <
typename somenumber>
853 const Vector<somenumber> &src,
862 template <
typename somenumber>
867 somenumber * partial_sum)
const;
874 template <
typename somenumber>
877 const Vector<somenumber> &v,
880 somenumber * partial_sum)
const;
913 template <
typename number>
915 const number & value)
922 template <
typename number>
929 template <
typename number>
938 template <
typename number>
942 const unsigned short i)
949 template <
typename number>
957 template <
typename number>
965 template <
typename number>
966 inline unsigned short 974 template <
typename number>
982 template <
typename number>
986 const unsigned short i)
1016 template <
typename number>
1043 template <
typename number>
1051 template <
typename number>
1059 template <
typename number>
1069 template <
typename number>
1074 return !(*
this == other);
1078 template <
typename number>
1090 template <
typename number>
1098 template <
typename number>
1106 template <
typename number>
1115 for (
size_type i = r.start; i < end; ++i)
1117 Entry *
const entry = &
data[i];
1118 if (entry->column == col)
1128 template <
typename number>
1133 return t->
locate(row, col);
1137 template <
typename number>
1154 while (i < end &&
data[i].column < col)
1158 if (i != end &&
data[i].column == col)
1171 if (end >=
row_info[row + 1].start)
1186 if (end >=
data.size())
1191 data.push_back(Entry());
1195 Entry *entry = &
data[i];
1197 Entry temp = *entry;
1202 entry->column = col;
1207 r.diagonal = i - r.start;
1235 template <
typename number>
1240 const bool elide_zero_values)
1247 if (elide_zero_values && value == 0.)
1249 Entry *entry =
locate(i, j);
1250 if (entry !=
nullptr)
1256 entry->value = value;
1262 template <
typename number>
1274 if (std::abs(value) == 0.)
1278 entry->value += value;
1282 template <
typename number>
1283 template <
typename number2>
1287 const bool elide_zero_values)
1290 for (
size_type i = 0; i < indices.size(); ++i)
1291 for (
size_type j = 0; j < indices.size(); ++j)
1292 if ((full_matrix(i, j) != 0) || (elide_zero_values ==
false))
1293 add(indices[i], indices[j], full_matrix(i, j));
1298 template <
typename number>
1299 template <
typename number2>
1302 const std::vector<size_type> &col_indices,
1304 const bool elide_zero_values)
1307 for (
size_type i = 0; i < row_indices.size(); ++i)
1308 for (
size_type j = 0; j < col_indices.size(); ++j)
1309 if ((full_matrix(i, j) != 0) || (elide_zero_values ==
false))
1310 add(row_indices[i], col_indices[j], full_matrix(i, j));
1315 template <
typename number>
1316 template <
typename number2>
1319 const std::vector<size_type> &col_indices,
1320 const std::vector<number2> &
values,
1321 const bool elide_zero_values)
1324 for (
size_type j = 0; j < col_indices.size(); ++j)
1325 if ((values[j] != 0) || (elide_zero_values ==
false))
1326 add(row, col_indices[j], values[j]);
1331 template <
typename number>
1332 template <
typename number2>
1338 const bool elide_zero_values,
1343 if ((std::abs(values[j]) != 0) || (elide_zero_values ==
false))
1344 add(row, col_indices[j], values[j]);
1349 template <
typename number>
1353 const Entry *entry =
locate(i, j);
1355 return entry->value;
1361 template <
typename number>
1365 const Entry *entry =
locate(i, j);
1367 return entry->value;
1373 template <
typename number>
1381 template <
typename number>
1388 template <
typename number>
1397 template <
typename number>
1406 template <
typename number>
1407 template <
typename MatrixType>
1410 const bool elide_zero_values)
1417 for (
size_type row = 0; row < M.m(); ++row)
1419 const typename MatrixType::const_iterator end_row = M.end(row);
1420 for (
typename MatrixType::const_iterator entry = M.begin(row);
1423 set(row, entry->column(), entry->value(), elide_zero_values);
1429 template <
typename number>
1430 template <
typename MatrixType>
1443 for (
size_type row = 0; row < M.m(); ++row)
1445 const typename MatrixType::const_iterator end_row = M.end(row);
1446 for (
typename MatrixType::const_iterator entry = M.begin(row);
1449 if (entry->value() != 0)
1450 add(row, entry->column(), factor * entry->value());
1456 template <
typename number>
1457 template <
typename MatrixTypeA,
typename MatrixTypeB>
1460 const MatrixTypeB &B,
1479 typename MatrixTypeB::const_iterator b1 = B.begin();
1480 const typename MatrixTypeB::const_iterator b_final = B.end();
1482 while (b1 != b_final)
1486 typename MatrixTypeB::const_iterator b2 = B.begin();
1487 while (b2 != b_final)
1492 const typename MatrixTypeA::value_type a = A.el(k, l);
1495 add(i, j, a * b1->value() * b2->value());
1506 std::vector<size_type> minrow(B.n(), B.m());
1507 std::vector<size_type> maxrow(B.n(), 0);
1508 while (b1 != b_final)
1511 if (r < minrow[b1->column()])
1512 minrow[b1->column()] = r;
1513 if (r > maxrow[b1->column()])
1514 maxrow[b1->column()] = r;
1518 typename MatrixTypeA::const_iterator ai = A.begin();
1519 const typename MatrixTypeA::const_iterator ae = A.end();
1523 const typename MatrixTypeA::value_type a = ai->value();
1533 b1 = B.begin(minrow[ai->row()]);
1534 const typename MatrixTypeB::const_iterator be1 =
1535 B.end(maxrow[ai->row()]);
1536 const typename MatrixTypeB::const_iterator be2 =
1537 B.end(maxrow[ai->column()]);
1541 const double b1v = b1->value();
1546 if (b1->column() == ai->row() && (b1v != 0.))
1550 typename MatrixTypeB::const_iterator b2 =
1551 B.begin(minrow[ai->column()]);
1554 if (b2->column() == ai->column())
1557 add(i, j, a * b1v * b2->value());
1570 template <
typename number>
1571 template <
class StreamType>
1578 std::vector<size_type> used_by_line;
1582 out <<
"SparseMatrixEZ:used entries:" << used << std::endl
1583 <<
"SparseMatrixEZ:allocated entries:" << allocated << std::endl
1584 <<
"SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1588 for (
size_type i = 0; i < used_by_line.size(); ++i)
1589 if (used_by_line[i] != 0)
1590 out <<
"SparseMatrixEZ:entries\t" << i <<
"\trows\t" 1591 << used_by_line[i] << std::endl;
size_type get_row_length(const size_type row) const
const types::global_dof_index invalid_size_type
void Tvmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
const_iterator & operator++()
static ::ExceptionBase & ExcEntryAllocationFailure(int arg1, int arg2)
#define DeclException2(Exception2, type1, type2, outsequence)
void add(const size_type i, const size_type j, const number value)
Contents is actually a matrix.
bool operator!=(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void compute_statistics(size_type &used, size_type &allocated, size_type &reserved, std::vector< size_type > &used_by_line, const bool compute_by_line) const
#define AssertIndexRange(index, range)
number operator()(const size_type i, const size_type j) const
SparseMatrixEZ< number > & copy_from(const MatrixType &source, const bool elide_zero_values=true)
void set(const size_type i, const size_type j, const number value, const bool elide_zero_values=true)
bool operator<(const SynchronousIterators< Iterators > &a, const SynchronousIterators< Iterators > &b)
void conjugate_add(const MatrixTypeA &A, const MatrixTypeB &B, const bool transpose=false)
void block_read(std::istream &in)
bool operator<(const const_iterator &) const
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
bool operator==(const AlignedVector< T > &lhs, const AlignedVector< T > &rhs)
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
static const size_type invalid
const Accessor & operator*() const
std::size_t memory_consumption() const
void threaded_matrix_norm_square(const Vector< somenumber > &v, const size_type begin_row, const size_type end_row, somenumber *partial_sum) const
void threaded_matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v, const size_type begin_row, const size_type end_row, somenumber *partial_sum) const
bool operator==(const const_iterator &) const
void print(std::ostream &out) const
static ::ExceptionBase & ExcNoDiagonal()
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
SparseMatrixEZ< number > & operator=(const SparseMatrixEZ< number > &)
Accessor(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
#define Assert(cond, exc)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
void reinit(const size_type n_rows, const size_type n_columns, size_type default_row_length=0, unsigned int default_increment=1, size_type reserve=0)
bool operator!=(const const_iterator &) const
const_iterator(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
#define DeclException0(Exception0)
#define DEAL_II_NAMESPACE_CLOSE
const SparseMatrixEZ< number > * matrix
std::vector< Entry > data
SynchronousIterators< Iterators > operator++(SynchronousIterators< Iterators > &a)
Entry * allocate(const size_type row, const size_type col)
RowInfo(const size_type start=Entry::invalid)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void print_statistics(StreamType &s, bool full=false)
static ::ExceptionBase & ExcIteratorPastEnd()
unsigned int global_dof_index
void block_write(std::ostream &out) const
void swap(MemorySpaceData< Number, MemorySpace > &, MemorySpaceData< Number, MemorySpace > &)
const Accessor * operator->() const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
unsigned short index() const
number el(const size_type i, const size_type j) const
#define DEAL_II_NAMESPACE_OPEN
const Entry * locate(const size_type row, const size_type col) const
const_iterator begin() const
static ::ExceptionBase & ExcInvalidEntry(int arg1, int arg2)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
void threaded_vmult(Vector< somenumber > &dst, const Vector< somenumber > &src, const size_type begin_row, const size_type end_row) const
~SparseMatrixEZ() override=default
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
std::enable_if< std::is_floating_point< T >::value &&std::is_floating_point< U >::value, typename ProductType< std::complex< T >, std::complex< U > >::type >::type operator*(const std::complex< T > &left, const std::complex< U > &right)
static const unsigned short invalid_diagonal
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
const_iterator end() const
#define AssertIsFinite(number)
unsigned int saved_default_row_length
std::vector< RowInfo > row_info
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
static ::ExceptionBase & ExcInternalError()
size_type n_nonzero_elements() const
void print_formatted(std::ostream &out, const unsigned int precision=3, const bool scientific=true, const unsigned int width=0, const char *zero_string=" ", const double denominator=1.) const