18 #ifndef __deal2__sparse_matrix_templates_h 19 #define __deal2__sparse_matrix_templates_h 22 #include <deal.II/base/config.h> 23 #include <deal.II/base/template_constraints.h> 24 #include <deal.II/base/parallel.h> 25 #include <deal.II/base/thread_management.h> 26 #include <deal.II/base/multithread_info.h> 27 #include <deal.II/base/utilities.h> 28 #include <deal.II/lac/sparse_matrix.h> 29 #include <deal.II/lac/trilinos_sparse_matrix.h> 30 #include <deal.II/lac/vector.h> 31 #include <deal.II/lac/full_matrix.h> 32 #include <deal.II/lac/compressed_simple_sparsity_pattern.h> 33 #include <deal.II/lac/vector_memory.h> 42 #include <deal.II/base/std_cxx1x/bind.h> 49 template <
typename number>
52 cols(0,
"SparseMatrix"),
59 template <
typename number>
63 cols(0,
"SparseMatrix"),
74 template <
typename number>
87 template <
typename number>
90 cols(0,
"SparseMatrix"),
99 template <
typename number>
103 cols(0,
"SparseMatrix"),
117 template <
typename number>
135 void zero_subrange (
const size_type
begin,
139 std::memset (dst+begin,0,(end-begin)*
sizeof(T));
146 template <
typename number>
164 internal::SparseMatrix::minimum_parallel_grain_size *
166 if (matrix_size>grain_size)
168 std_cxx1x::bind(&internal::SparseMatrix::template
169 zero_subrange<number>,
170 std_cxx1x::_1, std_cxx1x::_2,
173 else if (matrix_size > 0)
174 std::memset (&
val[0], 0, matrix_size*
sizeof(number));
181 template <
typename number>
199 template <
typename number>
228 template <
typename number>
240 template <
typename number>
252 template <
typename number>
262 template <
typename number>
272 template <
typename number>
281 if (std::fabs(
val[i]) > threshold)
288 template <
typename number>
306 while ((val_ptr != val_end_of_row) && (*colnum_ptr<row))
310 const number mean_value = (*val_ptr +
311 val[(*cols)(*colnum_ptr,row)]) / 2.0;
315 *val_ptr = mean_value;
316 set (*colnum_ptr, row, mean_value);
327 template <
typename number>
328 template <
typename somenumber>
344 template <
typename number>
345 template <
typename somenumber>
355 if (matrix(row,col) != 0)
356 set (row, col, matrix(row,col));
361 #ifdef DEAL_II_WITH_TRILINOS 363 template <
typename number>
373 std::vector < TrilinosScalar > value_cache;
374 std::vector<size_type> colnum_cache;
376 for (
size_type row = 0; row < matrix.
m(); ++row)
378 value_cache.resize(matrix.
n());
379 colnum_cache.resize(matrix.
n());
388 reinterpret_cast<TrilinosWrappers::types::int_type *>(&(colnum_cache[0])));
389 Assert (ierr==0, ExcTrilinosError(ierr));
392 value_cache.resize(ncols);
393 colnum_cache.resize(ncols);
407 template <
typename number>
408 template <
typename somenumber>
417 number *val_ptr = &
val[0];
418 const somenumber *matrix_ptr = &matrix.
val[0];
421 while (val_ptr != end_ptr)
422 *val_ptr++ += factor **matrix_ptr++;
439 template <
typename number,
442 void vmult_on_subrange (
const size_type begin_row,
444 const number *values,
445 const std::size_t *rowstart,
451 const number *val_ptr = &values[rowstart[begin_row]];
452 const size_type *colnum_ptr = &colnums[rowstart[begin_row]];
453 typename OutVector::iterator dst_ptr = dst.begin() + begin_row;
456 for (
size_type row=begin_row; row<end_row; ++row)
458 typename OutVector::value_type s = 0.;
459 const number *
const val_end_of_row = &values[rowstart[row+1]];
460 while (val_ptr != val_end_of_row)
461 s += *val_ptr++ * src(*colnum_ptr++);
465 for (
size_type row=begin_row; row<end_row; ++row)
467 typename OutVector::value_type s = *dst_ptr;
468 const number *
const val_end_of_row = &values[rowstart[row+1]];
469 while (val_ptr != val_end_of_row)
470 s += *val_ptr++ * src(*colnum_ptr++);
478 template <
typename number>
479 template <
typename number2>
484 const number2 *values,
485 const bool elide_zero_values,
486 const bool col_indices_are_sorted)
495 if (elide_zero_values ==
false && col_indices_are_sorted ==
true &&
502 Assert (col_indices[i] > col_indices[i-1],
503 ExcMessage(
"List of indices is unsorted or contains duplicates."));
520 const size_type diag = diag_pos - col_indices;
522 if (diag != n_cols && *diag_pos == row)
524 val_ptr[0] += *(values+(diag_pos-col_indices));
532 while (this_cols[counter]<col_indices[i] && counter<row_length_1)
535 Assert (this_cols[counter] == col_indices[i] || values[i] == 0,
536 ExcInvalidIndex(row,col_indices[i]));
542 for (
size_type i=post_diag; i<n_cols; ++i)
544 while (this_cols[counter]<col_indices[i] && counter<row_length_1)
547 Assert (this_cols[counter] == col_indices[i] || values[i] == 0,
548 ExcInvalidIndex(row,col_indices[i]));
553 Assert (counter < cols->row_length(row),
554 ExcMessage(
"Specified invalid column indices in add function."));
561 while (this_cols[counter]<col_indices[i] && counter<row_length_1)
564 Assert (this_cols[counter] == col_indices[i] || values[i] == 0,
565 ExcInvalidIndex(row,col_indices[i]));
569 Assert (counter < cols->row_length(row),
570 ExcMessage(
"Specified invalid column indices in add function."));
584 const number value = values[j];
588 if (elide_zero_values==
true && value == 0)
599 if (index < next_row_index && my_cols[index] == col_indices[j])
602 index =
cols->operator()(row, col_indices[j]);
610 Assert (value == 0., ExcInvalidIndex(row,col_indices[j]));
622 template <
typename number>
623 template <
typename number2>
628 const number2 *values,
629 const bool elide_zero_values)
632 Assert (row <
m(), ExcInvalidIndex1(row));
638 std::size_t index =
cols->
rowstart[row], next_index = index;
639 const std::size_t next_row_index =
cols->
rowstart[row+1];
641 if (elide_zero_values ==
true)
645 const number value = values[j];
655 if (index != next_row_index && my_cols[index] == col_indices[j])
658 next_index =
cols->operator()(row, col_indices[j]);
666 Assert (
false, ExcInvalidIndex(row,col_indices[j]));
681 const number value = values[j];
684 if (index != next_row_index && my_cols[index] == col_indices[j])
685 goto set_value_checked;
687 next_index =
cols->operator()(row, col_indices[j]);
691 Assert (value == 0., ExcInvalidIndex(row,col_indices[j]));
705 template <
typename number>
706 template <
class OutVector,
class InVector>
709 const InVector &src)
const 719 std_cxx1x::bind (&internal::SparseMatrix::vmult_on_subrange
720 <number,InVector,OutVector>,
721 std_cxx1x::_1, std_cxx1x::_2,
725 std_cxx1x::cref(src),
728 internal::SparseMatrix::minimum_parallel_grain_size);
733 template <
typename number>
734 template <
class OutVector,
class InVector>
737 const InVector &src)
const 753 dst(p) +=
val[j] * src(i);
760 template <
typename number>
761 template <
class OutVector,
class InVector>
764 const InVector &src)
const 774 std_cxx1x::bind (&internal::SparseMatrix::vmult_on_subrange
775 <number,InVector,OutVector>,
776 std_cxx1x::_1, std_cxx1x::_2,
780 std_cxx1x::cref(src),
783 internal::SparseMatrix::minimum_parallel_grain_size);
788 template <
typename number>
789 template <
class OutVector,
class InVector>
792 const InVector &src)
const 805 dst(p) +=
val[j] * src(i);
825 template <
typename number,
827 number matrix_norm_sqr_on_subrange (
const size_type begin_row,
829 const number *values,
830 const std::size_t *rowstart,
836 for (
size_type i=begin_row; i<end_row; ++i)
839 for (
size_type j=rowstart[i]; j<rowstart[i+1] ; j++)
840 s += values[j] * v(colnums[j]);
850 template <
typename number>
851 template <
typename somenumber>
861 parallel::accumulate_from_subranges<number>
862 (std_cxx1x::bind (&internal::SparseMatrix::matrix_norm_sqr_on_subrange
864 std_cxx1x::_1, std_cxx1x::_2,
868 internal::SparseMatrix::minimum_parallel_grain_size);
888 template <
typename number,
890 number matrix_scalar_product_on_subrange (
const size_type begin_row,
892 const number *values,
893 const std::size_t *rowstart,
900 for (
size_type i=begin_row; i<end_row; ++i)
903 for (
size_type j=rowstart[i]; j<rowstart[i+1] ; j++)
904 s += values[j] * v(colnums[j]);
914 template <
typename number>
915 template <
typename somenumber>
926 parallel::accumulate_from_subranges<number>
927 (std_cxx1x::bind (&internal::SparseMatrix::matrix_scalar_product_on_subrange
929 std_cxx1x::_1, std_cxx1x::_2,
934 internal::SparseMatrix::minimum_parallel_grain_size);
939 template <
typename number>
940 template <
typename numberB,
typename numberC>
945 const bool rebuild_sparsity_C)
const 947 const bool use_vector = V.
size() ==
n() ?
true :
false;
957 if (rebuild_sparsity_C ==
true)
962 ExcMessage (
"Can't use the same sparsity pattern for " 963 "different matrices if it is to be rebuilt."));
965 ExcMessage (
"Can't use the same sparsity pattern for " 966 "different matrices if it is to be rebuilt."));
981 for (
size_type i = 0; i < csp.n_rows(); ++i)
986 for (; rows != end_rows; ++rows)
1002 csp.add_entries (i, new_cols, end_new_cols,
true);
1005 sp_C.copy_from (csp);
1018 unsigned int max_n_cols_B = 0;
1020 max_n_cols_B = std::max (max_n_cols_B, sp_B.
row_length(i));
1021 std::vector<numberC> new_entries(max_n_cols_B);
1031 for (; rows != end_rows; ++rows)
1041 C.
add (i, *new_cols, A_val *
1043 (use_vector ? V(col) : 1));
1050 numberC *new_ptr = &new_entries[0];
1051 const numberB *B_val_ptr =
1053 const numberB *
const end_cols = &B.
val[sp_B.
rowstart[col+1]];
1054 for (; B_val_ptr != end_cols; ++B_val_ptr)
1055 *new_ptr++ = A_val **B_val_ptr * (use_vector ? V(col) : 1);
1057 C.
add (i, new_ptr-&new_entries[0], new_cols, &new_entries[0],
1066 template <
typename number>
1067 template <
typename numberB,
typename numberC>
1072 const bool rebuild_sparsity_C)
const 1074 const bool use_vector = V.
size() ==
m() ?
true :
false;
1084 if (rebuild_sparsity_C ==
true)
1089 ExcMessage (
"Can't use the same sparsity pattern for " 1090 "different matrices if it is to be rebuilt."));
1092 ExcMessage (
"Can't use the same sparsity pattern for " 1093 "different matrices if it is to be rebuilt."));
1099 sp_C.reinit (0,0,0);
1123 for (; rows != end_rows; ++rows)
1132 csp.add_entries (row, new_cols, end_new_cols,
true);
1135 sp_C.copy_from (csp);
1148 unsigned int max_n_cols_B = 0;
1150 max_n_cols_B = std::max (max_n_cols_B, sp_B.
row_length(i));
1151 std::vector<numberC> new_entries(max_n_cols_B);
1165 const numberB *
const end_cols = &B.
val[sp_B.
rowstart[i+1]];
1167 for (; rows != end_rows; ++rows)
1174 C.
add (row, i, A_val *
1176 (use_vector ? V(i) : 1));
1181 numberC *new_ptr = &new_entries[0];
1182 const numberB *B_val_ptr =
1184 for (; B_val_ptr != end_cols; ++B_val_ptr)
1185 *new_ptr++ = A_val **B_val_ptr * (use_vector ? V(i) : 1);
1187 C.
add (row, new_ptr-&new_entries[0], new_cols, &new_entries[0],
1195 template <
typename number>
1204 for (
size_type row=0; row<n_rows; ++row)
1213 template <
typename number>
1224 for (
size_type row=0; row<n_rows; ++row)
1228 while (val_ptr != val_end_of_row)
1238 template <
typename number>
1247 for (
const number *ptr = &
val[0];
1251 return std::sqrt (norm_sqr);
1271 template <
typename number,
1274 number residual_sqr_on_subrange (
const size_type begin_row,
1276 const number *values,
1277 const std::size_t *rowstart,
1285 for (
size_type i=begin_row; i<end_row; ++i)
1288 for (
size_type j=rowstart[i]; j<rowstart[i+1] ; j++)
1289 s -= values[j] * u(colnums[j]);
1299 template <
typename number>
1300 template <
typename somenumber>
1312 Assert (&u != &dst, ExcSourceEqualsDestination());
1315 std::sqrt (parallel::accumulate_from_subranges<number>
1316 (std_cxx1x::bind (&internal::SparseMatrix::residual_sqr_on_subrange
1318 std_cxx1x::_1, std_cxx1x::_2,
1322 std_cxx1x::ref(dst)),
1324 internal::SparseMatrix::minimum_parallel_grain_size));
1329 template <
typename number>
1330 template <
typename somenumber>
1334 const number om)
const 1343 somenumber *dst_ptr = dst.
begin();
1344 const somenumber *src_ptr = src.
begin();
1359 for (
size_type i=0; i<
n; ++i, ++dst_ptr, ++src_ptr, ++rowstart_ptr)
1360 *dst_ptr = om **src_ptr /
val[*rowstart_ptr];
1362 for (
size_type i=0; i<
n; ++i, ++dst_ptr, ++src_ptr, ++rowstart_ptr)
1363 *dst_ptr = *src_ptr /
val[*rowstart_ptr];
1368 template <
typename number>
1369 template <
typename somenumber>
1374 const std::vector<std::size_t> &pos_right_of_diagonal)
const 1388 somenumber *dst_ptr = &dst(0);
1393 if (pos_right_of_diagonal.size() != 0)
1395 Assert (pos_right_of_diagonal.size() == dst.
size(),
1399 for (
size_type row=0; row<
n; ++row, ++dst_ptr, ++rowstart_ptr)
1401 *dst_ptr = src(row);
1402 const std::size_t first_right_of_diagonal_index =
1403 pos_right_of_diagonal[row];
1404 Assert (first_right_of_diagonal_index <= *(rowstart_ptr+1),
1407 for (
size_type j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
1413 *dst_ptr /=
val[*rowstart_ptr];
1418 for ( ; rowstart_ptr!=&
cols->
rowstart[
n]; ++rowstart_ptr, ++dst_ptr)
1419 *dst_ptr *= om*(2.-om)*
val[*rowstart_ptr];
1423 dst_ptr = &dst(n-1);
1424 for (
int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
1426 const size_type end_row = *(rowstart_ptr+1);
1427 const size_type first_right_of_diagonal_index
1428 = pos_right_of_diagonal[row];
1430 for (
size_type j=first_right_of_diagonal_index; j<end_row; ++j)
1435 *dst_ptr /=
val[*rowstart_ptr];
1444 for (
size_type row=0; row<
n; ++row, ++dst_ptr, ++rowstart_ptr)
1446 *dst_ptr = src(row);
1454 const size_type first_right_of_diagonal_index
1462 for (
size_type j=(*rowstart_ptr)+1; j<first_right_of_diagonal_index; ++j)
1468 *dst_ptr /=
val[*rowstart_ptr];
1473 for (
size_type row=0; row<
n; ++row, ++rowstart_ptr, ++dst_ptr)
1474 *dst_ptr *= (2.-om)*
val[*rowstart_ptr];
1478 dst_ptr = &dst(n-1);
1479 for (
int row=n-1; row>=0; --row, --rowstart_ptr, --dst_ptr)
1481 const size_type end_row = *(rowstart_ptr+1);
1482 const size_type first_right_of_diagonal_index
1485 static_cast<size_type>(row)) -
1488 for (
size_type j=first_right_of_diagonal_index; j<end_row; ++j)
1492 *dst_ptr /=
val[*rowstart_ptr];
1497 template <
typename number>
1498 template <
typename somenumber>
1502 const number om)
const 1512 template <
typename number>
1513 template <
typename somenumber>
1517 const number om)
const 1527 template <
typename number>
1528 template <
typename somenumber>
1531 const number om)
const 1540 somenumber s = dst(row);
1545 s -=
val[j] * dst(col);
1554 template <
typename number>
1555 template <
typename somenumber>
1558 const number om)
const 1568 somenumber s = dst(row);
1584 template <
typename number>
1585 template <
typename somenumber>
1588 const std::vector<size_type> &permutation,
1589 const std::vector<size_type> &inverse_permutation,
1590 const number om)
const 1597 Assert (
m() == permutation.size(),
1599 Assert (
m() == inverse_permutation.size(),
1604 const size_type row = permutation[urow];
1605 somenumber s = dst(row);
1610 if (inverse_permutation[col] < urow)
1612 s -=
val[j] * dst(col);
1622 template <
typename number>
1623 template <
typename somenumber>
1626 const std::vector<size_type> &permutation,
1627 const std::vector<size_type> &inverse_permutation,
1628 const number om)
const 1635 Assert (
m() == permutation.size(),
1637 Assert (
m() == inverse_permutation.size(),
1643 const size_type row = permutation[urow];
1644 somenumber s = dst(row);
1648 if (inverse_permutation[col] > urow)
1649 s -=
val[j] * dst(col);
1659 template <
typename number>
1660 template <
typename somenumber>
1664 const number om)
const 1690 template <
typename number>
1691 template <
typename somenumber>
1695 const number om)
const 1705 somenumber s = b(row);
1717 template <
typename number>
1718 template <
typename somenumber>
1722 const number om)
const 1730 for (
int row=
m()-1; row>=0; --row)
1732 somenumber s = b(row);
1744 template <
typename number>
1745 template <
typename somenumber>
1749 const number om)
const 1757 template <
typename number>
1758 template <
typename somenumber>
1761 const number om)
const 1783 if (i>j) s +=
val[j] * dst(p);
1791 for (
int i=n-1; i>=0; i--)
1799 if (static_cast<size_type>(i)<j) s +=
val[j] * dst(p);
1809 template <
typename number>
1819 template <
typename number>
1821 const unsigned int precision,
1822 const bool scientific,
1823 const unsigned int width_,
1824 const char *zero_string,
1825 const double denominator)
const 1830 unsigned int width = width_;
1832 std::ios::fmtflags old_flags = out.flags();
1833 unsigned int old_precision = out.precision (precision);
1837 out.setf (std::ios::scientific, std::ios::floatfield);
1839 width = precision+7;
1843 out.setf (std::ios::fixed, std::ios::floatfield);
1845 width = precision+2;
1852 out << std::setw(width)
1853 <<
val[
cols->operator()(i,j)] * denominator <<
' ';
1855 out << std::setw(width) << zero_string <<
' ';
1861 out.precision(old_precision);
1862 out.flags (old_flags);
1867 template <
typename number>
1869 const double threshold)
const 1879 else if (std::fabs(
val[
cols->operator()(i,j)]) > threshold)
1890 template <
typename number>
1898 out <<
'[' <<
max_len <<
"][";
1900 out.write (reinterpret_cast<const char *>(&
val[0]),
1901 reinterpret_cast<const char *>(&
val[max_len])
1902 - reinterpret_cast<const char *>(&
val[0]));
1910 template <
typename number>
1933 in.read (reinterpret_cast<char *>(&
val[0]),
1934 reinterpret_cast<char *>(&
val[max_len])
1935 - reinterpret_cast<char *>(&
val[0]));
1941 template <
typename number>
1947 template <
typename number>
1951 return max_len*
static_cast<std::size_t
>(
sizeof(number)) +
sizeof(*this);
1955 DEAL_II_NAMESPACE_CLOSE
somenumber matrix_norm_square(const Vector< somenumber > &v) const
Iterator lower_bound(Iterator first, Iterator last, const T &val)
void precondition_Jacobi(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1.) const
static const number & conjugate(const number &x)
void Tvmult(OutVector &dst, const InVector &src) const
void vmult(OutVector &dst, const InVector &src) const
void PSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
#define AssertDimension(dim1, dim2)
void Tvmult_add(OutVector &dst, const InVector &src) const
numbers::NumberTraits< number >::real_type real_type
void mmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
const Epetra_CrsMatrix & trilinos_matrix() const
void SSOR(Vector< somenumber > &v, const number omega=1.) const
DEAL_VOLATILE unsigned int counter
size_type get_row_length(const size_type row) const
::ExceptionBase & ExcMessage(std::string arg1)
void set(const size_type i, const size_type j, const number value)
void block_read(std::istream &in)
size_type n_actually_nonzero_elements(const double threshold=0.) const
size_type n_nonzero_elements() const
#define AssertThrow(cond, exc)
static real_type abs(const number &x)
void Tmmult(SparseMatrix< numberC > &C, const SparseMatrix< numberB > &B, const Vector< number > &V=Vector< number >(), const bool rebuild_sparsity_pattern=true) const
void SOR(Vector< somenumber > &v, const number om=1.) const
void apply_to_subranges(const RangeType &begin, const typename identity< RangeType >::type &end, const Function &f, const unsigned int grainsize)
bool is_finite(const double x)
const_iterator begin() const
void print_pattern(std::ostream &out, const double threshold=0.) const
::ExceptionBase & ExcIO()
real_type l1_norm() const
void precondition_TSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
SparseMatrix< number > & operator=(const SparseMatrix< number > &)
unsigned int global_dof_index
#define Assert(cond, exc)
std::size_t memory_consumption() const
real_type linfty_norm() const
const SparsityPattern & get_sparsity_pattern() const
static bool equal(const T *p1, const T *p2)
void precondition_SSOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number omega=1., const std::vector< std::size_t > &pos_right_of_diagonal=std::vector< std::size_t >()) const
unsigned int row_length(const size_type row) const
void precondition_SOR(Vector< somenumber > &dst, const Vector< somenumber > &src, const number om=1.) const
::ExceptionBase & ExcInvalidConstructorCall()
size_type n_nonzero_elements() const
somenumber residual(Vector< somenumber > &dst, const Vector< somenumber > &x, const Vector< somenumber > &b) const
somenumber matrix_scalar_product(const Vector< somenumber > &u, const Vector< somenumber > &v) const
void add(const size_type i, const size_type j, const number value)
void TPSOR(Vector< somenumber > &v, const std::vector< size_type > &permutation, const std::vector< size_type > &inverse_permutation, const number om=1.) const
const_iterator end() const
SparseMatrix< number > & copy_from(const SparseMatrix< somenumber > &source)
virtual void reinit(const SparsityPattern &sparsity)
void Jacobi_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
void vmult_add(OutVector &dst, const InVector &src) const
types::global_dof_index size_type
::ExceptionBase & ExcNumberNotFinite()
void TSOR(Vector< somenumber > &v, const number om=1.) const
void SSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
SmartPointer< const SparsityPattern, SparseMatrix< number > > cols
void block_write(std::ostream &out) const
::ExceptionBase & ExcNotImplemented()
::ExceptionBase & ExcNotInitialized()
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
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
void SOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
real_type frobenius_norm() const
unsigned int row_length(const size_type row) const
void compress(::VectorOperation::values)
real_type linfty_norm() const
static const size_type invalid_entry
void TSOR_step(Vector< somenumber > &v, const Vector< somenumber > &b, const number om=1.) const