17 #ifndef __deal2__sparse_matrix_ez_h
18 #define __deal2__sparse_matrix_ez_h
21 #include <deal.II/base/config.h>
22 #include <deal.II/base/subscriptor.h>
23 #include <deal.II/base/smartpointer.h>
24 #include <deal.II/lac/exceptions.h>
30 template<
typename number>
class Vector;
101 template <
typename number>
129 const number &
value);
179 static const unsigned short
205 const unsigned short index);
212 size_type
row()
const;
219 unsigned short index()
const;
231 number
value()
const;
262 const unsigned short index);
379 const size_type default_row_length = 0,
380 const unsigned int default_increment = 1);
428 void reinit (
const size_type n_rows,
430 size_type default_row_length = 0,
431 unsigned int default_increment = 1,
432 size_type reserve = 0);
460 size_type
m ()
const;
468 size_type
n ()
const;
498 template <
class STREAM>
515 size_type &allocated,
517 std::vector<size_type> &used_by_line,
518 const bool compute_by_line)
const;
533 void set (
const size_type i,
const size_type j,
545 void add (
const size_type i,
576 template <
typename number2>
577 void add (
const std::vector<size_type> &indices,
579 const bool elide_zero_values =
true);
588 template <
typename number2>
589 void add (
const std::vector<size_type> &row_indices,
590 const std::vector<size_type> &col_indices,
592 const bool elide_zero_values =
true);
611 template <
typename number2>
612 void add (
const size_type row,
613 const std::vector<size_type> &col_indices,
614 const std::vector<number2> &values,
615 const bool elide_zero_values =
true);
634 template <
typename number2>
635 void add (
const size_type row,
636 const size_type n_cols,
637 const size_type *col_indices,
638 const number2 *values,
639 const bool elide_zero_values =
true,
640 const bool col_indices_are_sorted =
false);
668 template <
class MATRIX>
683 template <
class MATRIX>
684 void add (
const number factor,
710 const size_type j)
const;
717 number
el (
const size_type i,
718 const size_type j)
const;
729 template <
typename somenumber>
741 template <
typename somenumber>
751 template <
typename somenumber>
763 template <
typename somenumber>
790 template <
typename somenumber>
793 const number omega = 1.)
const;
799 template <
typename somenumber>
802 const number om = 1.,
803 const std::vector<std::size_t> &pos_right_of_diagonal = std::vector<std::size_t>())
const;
810 template <
typename somenumber>
813 const number om = 1.)
const;
820 template <
typename somenumber>
823 const number om = 1.)
const;
839 template <
class MATRIXA,
class MATRIXB>
842 const bool transpose =
false);
886 void print (std::ostream &out)
const;
928 const unsigned int precision = 3,
929 const bool scientific =
true,
930 const unsigned int width = 0,
931 const char *zero_string =
" ",
932 const double denominator = 1.)
const;
976 <<
"The entry with index (" << arg1 <<
',' << arg2
977 <<
") does not exist.");
981 <<
"An entry with index (" << arg1 <<
',' << arg2
982 <<
") cannot be allocated.");
992 const size_type col)
const;
1001 const size_type col);
1007 const size_type col);
1018 template <
typename somenumber>
1021 const size_type begin_row,
1022 const size_type end_row)
const;
1035 template <
typename somenumber>
1037 const size_type begin_row,
1038 const size_type end_row,
1039 somenumber *partial_sum)
const;
1052 template <
typename somenumber>
1055 const size_type begin_row,
1056 const size_type end_row,
1057 somenumber *partial_sum)
const;
1095 template <
typename number>
1098 const number &value)
1106 template <
typename number>
1115 template <
typename number>
1121 diagonal(invalid_diagonal)
1126 template <
typename number>
1131 const unsigned short i)
1139 template <
typename number>
1148 template <
typename number>
1153 return matrix->data[matrix->row_info[a_row].start+a_index].column;
1157 template <
typename number>
1167 template <
typename number>
1172 return matrix->data[matrix->row_info[a_row].start+a_index].value;
1176 template <
typename number>
1181 const unsigned short i)
1211 template <
typename number>
1219 ++(accessor.a_index);
1223 if (accessor.a_index >= accessor.matrix->row_info[accessor.a_row].length)
1225 accessor.a_index = 0;
1232 while (accessor.a_row < accessor.matrix->m()
1233 && accessor.matrix->row_info[accessor.a_row].length == 0);
1239 template <
typename number>
1248 template <
typename number>
1257 template <
typename number>
1268 template <
typename number>
1274 return ! (*
this == other);
1278 template <
typename number>
1291 template <
typename number>
1299 template <
typename number>
1307 template <
typename number>
1311 const size_type col)
1318 for (size_type i=r.
start; i<end; ++i)
1321 if (entry->
column == col)
1331 template <
typename number>
1335 const size_type col)
const
1338 return t->
locate(row,col);
1342 template <
typename number>
1346 const size_type col)
1354 size_type i = r.
start;
1361 while (i<end &&
data[i].column < col) ++i;
1364 if (i != end &&
data[i].column == col)
1386 for (size_type rn=row+1; rn<
row_info.size(); ++rn)
1392 if (end >=
data.size())
1403 Entry temp = *entry;
1422 for (size_type j = i+1; j <
end; ++j)
1429 std::swap (
data[j], temp);
1440 template <
typename number>
1461 entry->
value = value;
1467 template <
typename number>
1484 entry->
value += value;
1488 template <
typename number>
1489 template <
typename number2>
1492 const bool elide_zero_values)
1495 for (size_type i=0; i<indices.size(); ++i)
1496 for (size_type j=0; j<indices.size(); ++j)
1497 if ((full_matrix(i,j) != 0) || (elide_zero_values ==
false))
1498 add (indices[i], indices[j], full_matrix(i,j));
1503 template <
typename number>
1504 template <
typename number2>
1506 const std::vector<size_type> &col_indices,
1508 const bool elide_zero_values)
1511 for (size_type i=0; i<row_indices.size(); ++i)
1512 for (size_type j=0; j<col_indices.size(); ++j)
1513 if ((full_matrix(i,j) != 0) || (elide_zero_values ==
false))
1514 add (row_indices[i], col_indices[j], full_matrix(i,j));
1520 template <
typename number>
1521 template <
typename number2>
1523 const std::vector<size_type> &col_indices,
1524 const std::vector<number2> &values,
1525 const bool elide_zero_values)
1528 for (size_type j=0; j<col_indices.size(); ++j)
1529 if ((values[j] != 0) || (elide_zero_values ==
false))
1530 add (row, col_indices[j], values[j]);
1535 template <
typename number>
1536 template <
typename number2>
1538 const size_type n_cols,
1539 const size_type *col_indices,
1540 const number2 *values,
1541 const bool elide_zero_values,
1545 for (size_type j=0; j<n_cols; ++j)
1546 if ((values[j] != 0) || (elide_zero_values ==
false))
1547 add (row, col_indices[j], values[j]);
1553 template <
typename number>
1556 const size_type j)
const
1560 return entry->
value;
1566 template <
typename number>
1569 const size_type j)
const
1573 return entry->
value;
1574 Assert(
false, ExcInvalidEntry(i,j));
1579 template <
typename number>
1588 template <
typename number>
1596 template <
typename number>
1606 template <
typename number>
1616 template<
typename number>
1617 template <
class MATRIX>
1627 for (size_type row = 0; row < M.m(); ++row)
1629 const typename MATRIX::const_iterator end_row = M.end(row);
1630 for (
typename MATRIX::const_iterator entry = M.begin(row);
1631 entry != end_row; ++entry)
1632 if (entry->value() != 0)
1633 set(row, entry->column(), entry->value());
1639 template<
typename number>
1640 template <
class MATRIX>
1655 for (size_type row = 0; row < M.m(); ++row)
1657 const typename MATRIX::const_iterator end_row = M.end(row);
1658 for (
typename MATRIX::const_iterator entry = M.begin(row);
1659 entry != end_row; ++entry)
1660 if (entry->value() != 0)
1661 add(row, entry->column(), factor * entry->value());
1667 template<
typename number>
1668 template <
class MATRIXA,
class MATRIXB>
1672 const bool transpose)
1690 typename MATRIXB::const_iterator b1 = B.begin();
1691 const typename MATRIXB::const_iterator b_final = B.end();
1693 while (b1 != b_final)
1695 const size_type i = b1->column();
1696 const size_type k = b1->row();
1697 typename MATRIXB::const_iterator b2 = B.begin();
1698 while (b2 != b_final)
1700 const size_type j = b2->column();
1701 const size_type l = b2->row();
1703 const typename MATRIXA::value_type a = A.el(k,l);
1706 add (i, j, a * b1->value() * b2->value());
1717 std::vector<size_type> minrow(B.n(), B.m());
1718 std::vector<size_type> maxrow(B.n(), 0);
1719 while (b1 != b_final)
1721 const size_type r = b1->row();
1722 if (r < minrow[b1->column()])
1723 minrow[b1->column()] = r;
1724 if (r > maxrow[b1->column()])
1725 maxrow[b1->column()] = r;
1729 typename MATRIXA::const_iterator ai = A.begin();
1730 const typename MATRIXA::const_iterator ae = A.end();
1734 const typename MATRIXA::value_type a = ai->value();
1737 if (a == 0.)
continue;
1743 b1 = B.begin(minrow[ai->row()]);
1744 const typename MATRIXB::const_iterator
1745 be1 = B.end(maxrow[ai->row()]);
1746 const typename MATRIXB::const_iterator
1747 be2 = B.end(maxrow[ai->column()]);
1751 const double b1v = b1->value();
1756 if (b1->column() == ai->row() && (b1v != 0.))
1758 const size_type i = b1->row();
1760 typename MATRIXB::const_iterator
1761 b2 = B.begin(minrow[ai->column()]);
1764 if (b2->column() == ai->column())
1766 const size_type j = b2->row();
1767 add (i, j, a * b1v * b2->value());
1780 template <
typename number>
1781 template <
class STREAM>
1787 size_type allocated;
1789 std::vector<size_type> used_by_line;
1793 out <<
"SparseMatrixEZ:used entries:" << used << std::endl
1794 <<
"SparseMatrixEZ:allocated entries:" << allocated << std::endl
1795 <<
"SparseMatrixEZ:reserved entries:" << reserved << std::endl;
1799 for (size_type i=0; i< used_by_line.size(); ++i)
1800 if (used_by_line[i] != 0)
1801 out <<
"SparseMatrixEZ:entries\t" << i
1802 <<
"\trows\t" << used_by_line[i]
1809 DEAL_II_NAMESPACE_CLOSE
Entry * allocate(const size_type row, const size_type col)
const types::global_dof_index invalid_size_type
unsigned short index() const
const_iterator & operator++()
const_iterator end() const
const_iterator begin() const
void block_read(std::istream &in)
void threaded_matrix_norm_square(const Vector< somenumber > &v, const size_type begin_row, const size_type end_row, somenumber *partial_sum) const
const Accessor & operator*() const
void threaded_vmult(Vector< somenumber > &dst, const Vector< somenumber > &src, const size_type begin_row, const size_type end_row) const
types::global_dof_index size_type
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
void vmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
bool operator<(const const_iterator &) const
number el(const size_type i, const size_type j) const
std::vector< RowInfo > row_info
void print_statistics(STREAM &s, bool full=false)
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 is_finite(const double x)
void add(const size_type i, const size_type j, const number value)
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
static const size_type invalid
Iterator is invalid, probably due to an error.
DeclException0(ExcNoDiagonal)
const Entry * locate(const size_type row, const size_type col) 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
Accessor(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
unsigned int global_dof_index
#define Assert(cond, exc)
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
const Accessor * operator->() const
void set(const size_type i, const size_type j, const number value)
const_iterator(const SparseMatrixEZ< number > *matrix, const size_type row, const unsigned short index)
void Tvmult_add(Vector< somenumber > &dst, const Vector< somenumber > &src) const
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const SparseMatrixEZ< number > * matrix
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)
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
void conjugate_add(const MATRIXA &A, const MATRIXB &B, const bool transpose=false)
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
RowInfo(size_type start=Entry::invalid)
void vmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
::ExceptionBase & ExcIteratorPastEnd()
::ExceptionBase & ExcNumberNotFinite()
size_type get_row_length(const size_type row) const
void Tvmult(Vector< somenumber > &dst, const Vector< somenumber > &src) const
number operator()(const size_type i, const size_type j) const
void print(std::ostream &out) const
SparseMatrixEZ< number > & operator=(const SparseMatrixEZ< number > &)
size_type n_nonzero_elements() const
std::size_t memory_consumption() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static const unsigned short invalid_diagonal
bool operator!=(const const_iterator &) const
::ExceptionBase & ExcInternalError()
void block_write(std::ostream &out) const
bool operator==(const const_iterator &) const
std::vector< Entry > data
DeclException2(ExcInvalidEntry, int, int,<< "The entry with index ("<< arg1<< ','<< arg2<< ") does not exist.")
SparseMatrixEZ< number > & copy_from(const MATRIX &source)