17 #ifndef __deal2__trilinos_sparse_matrix_h 18 #define __deal2__trilinos_sparse_matrix_h 21 #include <deal.II/base/config.h> 23 #ifdef DEAL_II_WITH_TRILINOS 25 # include <deal.II/base/std_cxx1x/shared_ptr.h> 26 # include <deal.II/base/subscriptor.h> 27 # include <deal.II/base/index_set.h> 28 # include <deal.II/lac/full_matrix.h> 29 # include <deal.II/lac/exceptions.h> 30 # include <deal.II/lac/trilinos_vector.h> 31 # include <deal.II/lac/vector_view.h> 37 # define TrilinosScalar double 38 # include <Epetra_FECrsMatrix.h> 39 # include <Epetra_Map.h> 40 # include <Epetra_CrsGraph.h> 41 # include <Epetra_MultiVector.h> 42 # ifdef DEAL_II_WITH_MPI 43 # include <Epetra_MpiComm.h> 46 # include "Epetra_SerialComm.h" 80 std::size_t, std::size_t, std::size_t,
81 <<
"You tried to access row " << arg1
82 <<
" of a distributed sparsity pattern, " 83 <<
" but only rows " << arg2 <<
" through " << arg3
84 <<
" are stored locally and can be accessed.");
115 const size_type
index);
121 size_type
row()
const;
127 size_type
index()
const;
212 template <
bool Constess>
218 TrilinosScalar value()
const;
223 TrilinosScalar &value();
255 template <
bool Other>
261 TrilinosScalar value()
const;
289 operator TrilinosScalar ()
const;
295 const Reference &operator = (
const TrilinosScalar n)
const;
301 const Reference &operator += (
const TrilinosScalar n)
const;
308 const Reference &operator -= (
const TrilinosScalar n)
const;
315 const Reference &operator *= (
const TrilinosScalar n)
const;
322 const Reference &operator /= (
const TrilinosScalar n)
const;
354 Reference value()
const;
366 friend class Reference;
383 template <
bool Constness>
407 const size_type
index);
413 template <
bool Other>
456 bool operator < (const Iterator<Constness> &)
const;
467 size_type, size_type,
468 <<
"Attempt to access element " << arg2
469 <<
" of row " << arg1
470 <<
" which doesn't have that many elements.");
479 template <
bool Other>
friend class Iterator;
542 static const bool zero_addition_can_be_elided =
true;
585 const unsigned int n_max_entries_per_row);
599 const std::vector<unsigned int> &n_entries_per_row);
651 template<
typename SparsityType>
652 void reinit (
const SparsityType &sparsity_pattern);
708 template <
typename number>
709 void reinit (const ::SparseMatrix<number> &dealii_sparse_matrix,
710 const double drop_tolerance=1e-13,
711 const bool copy_values=
true,
712 const ::SparsityPattern *use_this_sparsity=0);
721 void reinit (
const Epetra_CrsMatrix &input_matrix,
722 const bool copy_values =
true);
753 const size_type n_max_entries_per_row = 0);
769 const std::vector<unsigned int> &n_entries_per_row);
800 SparseMatrix (
const Epetra_Map &row_parallel_partitioning,
801 const Epetra_Map &col_parallel_partitioning,
802 const size_type n_max_entries_per_row = 0);
833 SparseMatrix (
const Epetra_Map &row_parallel_partitioning,
834 const Epetra_Map &col_parallel_partitioning,
835 const std::vector<unsigned int> &n_entries_per_row);
874 template<
typename SparsityType>
875 void reinit (
const Epetra_Map ¶llel_partitioning,
876 const SparsityType &sparsity_pattern,
877 const bool exchange_data =
false);
900 template<
typename SparsityType>
901 void reinit (
const Epetra_Map &row_parallel_partitioning,
902 const Epetra_Map &col_parallel_partitioning,
903 const SparsityType &sparsity_pattern,
904 const bool exchange_data =
false);
934 template <
typename number>
935 void reinit (
const Epetra_Map ¶llel_partitioning,
936 const ::SparseMatrix<number> &dealii_sparse_matrix,
937 const double drop_tolerance=1e-13,
938 const bool copy_values=
true,
939 const ::SparsityPattern *use_this_sparsity=0);
962 template <
typename number>
963 void reinit (
const Epetra_Map &row_parallel_partitioning,
964 const Epetra_Map &col_parallel_partitioning,
965 const ::SparseMatrix<number> &dealii_sparse_matrix,
966 const double drop_tolerance=1e-13,
967 const bool copy_values=
true,
968 const ::SparsityPattern *use_this_sparsity=0);
999 const MPI_Comm &communicator = MPI_COMM_WORLD,
1000 const unsigned int n_max_entries_per_row = 0);
1017 const MPI_Comm &communicator,
1018 const std::vector<unsigned int> &n_entries_per_row);
1048 const IndexSet &col_parallel_partitioning,
1049 const MPI_Comm &communicator = MPI_COMM_WORLD,
1050 const size_type n_max_entries_per_row = 0);
1082 const IndexSet &col_parallel_partitioning,
1083 const MPI_Comm &communicator,
1084 const std::vector<unsigned int> &n_entries_per_row);
1125 template<
typename SparsityType>
1126 void reinit (
const IndexSet ¶llel_partitioning,
1127 const SparsityType &sparsity_pattern,
1128 const MPI_Comm &communicator = MPI_COMM_WORLD,
1129 const bool exchange_data =
false);
1152 template<
typename SparsityType>
1153 void reinit (
const IndexSet &row_parallel_partitioning,
1154 const IndexSet &col_parallel_partitioning,
1155 const SparsityType &sparsity_pattern,
1156 const MPI_Comm &communicator = MPI_COMM_WORLD,
1157 const bool exchange_data =
false);
1187 template <
typename number>
1188 void reinit (
const IndexSet ¶llel_partitioning,
1189 const ::SparseMatrix<number> &dealii_sparse_matrix,
1190 const MPI_Comm &communicator = MPI_COMM_WORLD,
1191 const double drop_tolerance=1e-13,
1192 const bool copy_values=
true,
1193 const ::SparsityPattern *use_this_sparsity=0);
1216 template <
typename number>
1217 void reinit (
const IndexSet &row_parallel_partitioning,
1218 const IndexSet &col_parallel_partitioning,
1219 const ::SparseMatrix<number> &dealii_sparse_matrix,
1220 const MPI_Comm &communicator = MPI_COMM_WORLD,
1221 const double drop_tolerance=1e-13,
1222 const bool copy_values=
true,
1223 const ::SparsityPattern *use_this_sparsity=0);
1234 size_type m ()
const;
1240 size_type n ()
const;
1256 unsigned int local_size ()
const;
1274 std::pair<size_type, size_type>
1275 local_range ()
const;
1282 bool in_local_range (
const size_type
index)
const;
1288 size_type n_nonzero_elements ()
const;
1294 unsigned int row_length (
const size_type
row)
const;
1305 bool is_compressed ()
const;
1315 size_type memory_consumption ()
const;
1340 operator = (
const double d);
1397 void compress (::VectorOperation::values operation);
1430 void set (const size_type i,
1432 const TrilinosScalar value);
1470 void set (const
std::vector<size_type> &indices,
1471 const
FullMatrix<TrilinosScalar> &full_matrix,
1472 const
bool elide_zero_values = false);
1481 void set (const
std::vector<size_type> &row_indices,
1482 const
std::vector<size_type> &col_indices,
1483 const
FullMatrix<TrilinosScalar> &full_matrix,
1484 const
bool elide_zero_values = false);
1513 void set (const size_type
row,
1514 const
std::vector<size_type> &col_indices,
1515 const
std::vector<TrilinosScalar> &values,
1516 const
bool elide_zero_values = false);
1545 void set (const size_type row,
1546 const size_type n_cols,
1547 const size_type *col_indices,
1548 const TrilinosScalar *values,
1549 const
bool elide_zero_values = false);
1566 void add (const size_type i,
1568 const TrilinosScalar value);
1606 void add (const
std::vector<size_type> &indices,
1607 const
FullMatrix<TrilinosScalar> &full_matrix,
1608 const
bool elide_zero_values = true);
1617 void add (const
std::vector<size_type> &row_indices,
1618 const
std::vector<size_type> &col_indices,
1619 const
FullMatrix<TrilinosScalar> &full_matrix,
1620 const
bool elide_zero_values = true);
1648 void add (const size_type row,
1649 const
std::vector<size_type> &col_indices,
1650 const
std::vector<TrilinosScalar> &values,
1651 const
bool elide_zero_values = true);
1678 void add (const size_type row,
1679 const size_type n_cols,
1680 const size_type *col_indices,
1681 const TrilinosScalar *values,
1682 const
bool elide_zero_values = true,
1683 const
bool col_indices_are_sorted = false);
1689 SparseMatrix &operator *= (const TrilinosScalar factor);
1695 SparseMatrix &operator /= (const TrilinosScalar factor);
1715 void add (const TrilinosScalar factor,
1757 void clear_row (const size_type row,
1758 const TrilinosScalar new_diag_value = 0);
1776 void clear_rows (const
std::vector<size_type> &rows,
1777 const TrilinosScalar new_diag_value = 0);
1810 TrilinosScalar operator () (const size_type i,
1811 const size_type j) const;
1850 TrilinosScalar el (const size_type i,
1851 const size_type j) const;
1864 TrilinosScalar diag_element (const size_type i) const;
1893 template<typename VectorType>
1894 void vmult (VectorType &dst,
1895 const VectorType &src) const;
1919 template <typename VectorType>
1920 void Tvmult (VectorType &dst,
1921 const VectorType &src) const;
1946 template<typename VectorType>
1947 void vmult_add (VectorType &dst,
1948 const VectorType &src) const;
1973 template <typename VectorType>
1974 void Tvmult_add (VectorType &dst,
1975 const VectorType &src) const;
2026 TrilinosScalar matrix_norm_square (const
VectorBase &v) const;
2063 TrilinosScalar matrix_scalar_product (const
VectorBase &u,
2185 TrilinosScalar l1_norm () const;
2202 TrilinosScalar linfty_norm () const;
2211 TrilinosScalar frobenius_norm () const;
2224 const Epetra_CrsMatrix &trilinos_matrix () const;
2233 const Epetra_CrsGraph &trilinos_sparsity_pattern () const;
2244 const Epetra_Map &domain_partitioner () const;
2255 const Epetra_Map &range_partitioner () const;
2264 const Epetra_Map &row_partitioner () const;
2275 const Epetra_Map &col_partitioner () const;
2286 const_iterator begin () const;
2291 const_iterator end () const;
2306 const_iterator begin (const size_type r) const;
2323 const_iterator end (const size_type r) const;
2349 iterator begin (const size_type r);
2366 iterator end (const size_type r);
2383 void write_ascii ();
2397 void print (
std::ostream &out,
2398 const
bool write_extended_trilinos_info = false) const;
2410 << "An error with error number " << arg1
2411 << " occurred while calling a Trilinos function");
2417 size_type, size_type,
2418 << "The entry with
index <" << arg1 << ',' << arg2
2419 << "> does not exist.");
2435 size_type, size_type, size_type, size_type,
2436 << "You tried to access element (" << arg1
2437 << "/" << arg2 << ")"
2438 << " of a distributed matrix, but only rows "
2439 << arg3 << " through " << arg4
2440 << " are stored locally and can be accessed.");
2446 size_type, size_type,
2447 << "You tried to access element (" << arg1
2448 << "/" << arg2 << ")"
2449 << " of a sparse matrix, but it appears to not"
2450 << " exist in the Trilinos sparsity pattern.");
2536 Epetra_CombineMode last_action;
2577 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
2586 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
2595 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
2602 const size_type row,
2603 const size_type index)
2605 AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2609 template <
bool Other>
2621 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
2628 const Accessor<false> &acc)
2630 accessor(
const_cast<Accessor<false>&
>(acc))
2635 Accessor<false>::Reference::operator TrilinosScalar ()
const 2637 return (*accessor.value_cache)[accessor.a_index];
2644 (*accessor.value_cache)[accessor.a_index] = n;
2645 accessor.matrix->set(accessor.row(), accessor.column(),
2646 static_cast<TrilinosScalar
>(*this));
2655 (*accessor.value_cache)[accessor.a_index] += n;
2656 accessor.matrix->set(accessor.row(), accessor.column(),
2657 static_cast<TrilinosScalar
>(*this));
2666 (*accessor.value_cache)[accessor.a_index] -= n;
2667 accessor.matrix->set(accessor.row(), accessor.column(),
2668 static_cast<TrilinosScalar
>(*this));
2677 (*accessor.value_cache)[accessor.a_index] *= n;
2678 accessor.matrix->set(accessor.row(), accessor.column(),
2679 static_cast<TrilinosScalar
>(*this));
2688 (*accessor.value_cache)[accessor.a_index] /= n;
2689 accessor.matrix->set(accessor.row(), accessor.column(),
2690 static_cast<TrilinosScalar
>(*this));
2697 const size_type row,
2698 const size_type index)
2708 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
2709 return Reference(*
this);
2714 template <
bool Constness>
2716 Iterator<Constness>::Iterator(MatrixType *matrix,
2717 const size_type row,
2718 const size_type index)
2720 accessor(matrix, row, index)
2724 template <
bool Constness>
2725 template <
bool Other>
2727 Iterator<Constness>::Iterator(
const Iterator<Other> &other)
2729 accessor(other.accessor)
2733 template <
bool Constness>
2735 Iterator<Constness> &
2746 if (accessor.a_index >= accessor.colnum_cache->size())
2748 accessor.a_index = 0;
2751 while ((accessor.a_row < accessor.matrix->m())
2753 (accessor.matrix->row_length(accessor.a_row) == 0))
2756 accessor.visit_present_row();
2762 template <
bool Constness>
2767 const Iterator<Constness> old_state = *
this;
2774 template <
bool Constness>
2776 const Accessor<Constness> &
2784 template <
bool Constness>
2786 const Accessor<Constness> *
2794 template <
bool Constness>
2799 return (accessor.a_row == other.accessor.a_row &&
2800 accessor.a_index == other.accessor.a_index);
2805 template <
bool Constness>
2810 return ! (*
this == other);
2815 template <
bool Constness>
2820 return (accessor.row() < other.accessor.row() ||
2821 (accessor.row() == other.accessor.row() &&
2822 accessor.index() < other.accessor.index()));
2826 template <
bool Constness>
2831 return (other < *
this);
2842 return const_iterator(
this, 0, 0);
2851 return const_iterator(
this, m(), 0);
2861 if (row_length(r) > 0)
2862 return const_iterator(
this, r, 0);
2878 for (size_type i=r+1; i<m(); ++i)
2879 if (row_length(i) > 0)
2880 return const_iterator(
this, i, 0);
2893 return iterator(
this, 0, 0);
2902 return iterator(
this, m(), 0);
2912 if (row_length(r) > 0)
2913 return iterator(
this, r, 0);
2929 for (size_type i=r+1; i<m(); ++i)
2930 if (row_length(i) > 0)
2931 return iterator(
this, i, 0);
2945 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 2946 begin = matrix->RowMap().MinMyGID();
2947 end = matrix->RowMap().MaxMyGID()+1;
2949 begin = matrix->RowMap().MinMyGID();
2950 end = matrix->RowMap().MaxMyGID()+1;
2953 return ((index >= static_cast<size_type>(begin)) &&
2954 (index < static_cast<size_type>(end)));
2973 Epetra_CombineMode mode = last_action;
2974 if (last_action == Zero)
2976 if ((operation==::VectorOperation::add) ||
2977 (operation==::VectorOperation::unknown))
2979 else if (operation==::VectorOperation::insert)
2985 ((last_action == Add) && (operation!=::VectorOperation::insert))
2987 ((last_action == Insert) && (operation!=::VectorOperation::add)),
2988 ExcMessage(
"operation and argument to compress() do not match"));
2993 ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(),
2998 ierr = matrix->OptimizeStorage ();
3012 compress(::VectorOperation::unknown);
3022 compress (::VectorOperation::unknown);
3024 const int ierr = matrix->PutScalar(d);
3043 const TrilinosScalar value)
3048 set (i, 1, &j, &value,
false);
3057 const bool elide_zero_values)
3059 Assert (indices.size() == values.
m(),
3061 Assert (values.
m() == values.
n(), ExcNotQuadratic());
3063 for (size_type i=0; i<indices.size(); ++i)
3064 set (indices[i], indices.size(), &indices[0], &values(i,0),
3073 const std::vector<size_type> &col_indices,
3075 const bool elide_zero_values)
3077 Assert (row_indices.size() == values.
m(),
3079 Assert (col_indices.size() == values.
n(),
3082 for (size_type i=0; i<row_indices.size(); ++i)
3083 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
3092 const std::vector<size_type> &col_indices,
3093 const std::vector<TrilinosScalar> &values,
3094 const bool elide_zero_values)
3096 Assert (col_indices.size() == values.size(),
3099 set (
row, col_indices.size(), &col_indices[0], &values[0],
3108 const size_type n_cols,
3109 const size_type *col_indices,
3110 const TrilinosScalar *values,
3111 const bool elide_zero_values)
3114 if (last_action == Add)
3116 ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(),
3119 Assert (ierr == 0, ExcTrilinosError(ierr));
3122 last_action = Insert;
3125 TrilinosScalar *col_value_ptr;
3128 TrilinosScalar short_val_array[100];
3130 std::vector<TrilinosScalar> long_val_array;
3131 std::vector<TrilinosWrappers::types::int_type> long_index_array;
3137 if (elide_zero_values ==
false)
3140 col_value_ptr =
const_cast<TrilinosScalar *
>(values);
3149 long_val_array.resize(n_cols);
3150 long_index_array.resize(n_cols);
3151 col_index_ptr = &long_index_array[0];
3152 col_value_ptr = &long_val_array[0];
3156 col_index_ptr = &short_index_array[0];
3157 col_value_ptr = &short_val_array[0];
3161 for (size_type j=0; j<n_cols; ++j)
3163 const double value = values[j];
3167 col_index_ptr[n_columns] = col_indices[j];
3168 col_value_ptr[n_columns] = value;
3185 if (row_partitioner().MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) ==
true)
3187 if (matrix->Filled() ==
false)
3189 ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(
3190 static_cast<TrilinosWrappers::types::int_type>(row),
3191 static_cast<int>(n_columns),const_cast<double *>(col_value_ptr),
3201 ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_columns,
3215 if (matrix->Filled() ==
false)
3217 ierr = matrix->InsertGlobalValues (1,
3219 n_columns, col_index_ptr,
3221 Epetra_FECrsMatrix::ROW_MAJOR);
3226 ierr = matrix->ReplaceGlobalValues (1,
3228 n_columns, col_index_ptr,
3230 Epetra_FECrsMatrix::ROW_MAJOR);
3233 Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
3243 const TrilinosScalar value)
3256 if (last_action == Insert)
3259 ierr = matrix->GlobalAssemble(*column_space_map,
3260 row_partitioner(),
false);
3262 Assert (ierr == 0, ExcTrilinosError(ierr));
3271 add (i, 1, &j, &value,
false);
3280 const bool elide_zero_values)
3282 Assert (indices.size() == values.
m(),
3284 Assert (values.
m() == values.
n(), ExcNotQuadratic());
3286 for (size_type i=0; i<indices.size(); ++i)
3287 add (indices[i], indices.size(), &indices[0], &values(i,0),
3296 const std::vector<size_type> &col_indices,
3298 const bool elide_zero_values)
3300 Assert (row_indices.size() == values.
m(),
3302 Assert (col_indices.size() == values.
n(),
3305 for (size_type i=0; i<row_indices.size(); ++i)
3306 add (row_indices[i], col_indices.size(), &col_indices[0],
3307 &values(i,0), elide_zero_values);
3315 const std::vector<size_type> &col_indices,
3316 const std::vector<TrilinosScalar> &values,
3317 const bool elide_zero_values)
3319 Assert (col_indices.size() == values.size(),
3322 add (row, col_indices.size(), &col_indices[0], &values[0],
3331 const size_type n_cols,
3332 const size_type *col_indices,
3333 const TrilinosScalar *values,
3334 const bool elide_zero_values,
3338 if (last_action == Insert)
3342 ierr = matrix->GlobalAssemble(*column_space_map,
3343 row_partitioner(),
false);
3351 TrilinosScalar *col_value_ptr;
3354 double short_val_array[100];
3356 std::vector<TrilinosScalar> long_val_array;
3357 std::vector<TrilinosWrappers::types::int_type> long_index_array;
3362 if (elide_zero_values ==
false)
3365 col_value_ptr =
const_cast<TrilinosScalar *
>(values);
3368 for (size_type j=0; j<n_cols; ++j)
3378 long_val_array.resize(n_cols);
3379 long_index_array.resize(n_cols);
3380 col_index_ptr = &long_index_array[0];
3381 col_value_ptr = &long_val_array[0];
3385 col_index_ptr = &short_index_array[0];
3386 col_value_ptr = &short_val_array[0];
3390 for (size_type j=0; j<n_cols; ++j)
3392 const double value = values[j];
3397 col_index_ptr[n_columns] = col_indices[j];
3398 col_value_ptr[n_columns] = value;
3410 if (row_partitioner().MyGID(static_cast<TrilinosWrappers::types::int_type>(row)) ==
true)
3412 ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_columns,
3426 ierr = matrix->SumIntoGlobalValues (1,
3430 Epetra_FECrsMatrix::ROW_MAJOR);
3436 std::cout <<
"------------------------------------------" 3438 std::cout <<
"Got error " << ierr <<
" in row " << row
3439 <<
" of proc " << row_partitioner().Comm().MyPID()
3440 <<
" when trying to add the columns:" << std::endl;
3442 std::cout << col_index_ptr[i] <<
" ";
3443 std::cout << std::endl << std::endl;
3444 std::cout <<
"Matrix row has the following indices:" << std::endl;
3445 int n_indices, *indices;
3446 trilinos_sparsity_pattern().ExtractMyRowView(row_partitioner().LID(static_cast<TrilinosWrappers::types::int_type>(row)),
3450 std::cout << indices[i] <<
" ";
3451 std::cout << std::endl << std::endl;
3453 ExcAccessToNonPresentElement(row, col_index_ptr[0]));
3456 Assert (ierr >= 0, ExcTrilinosError(ierr));
3468 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 3469 return matrix->NumGlobalRows();
3471 return matrix->NumGlobalRows64();
3481 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 3482 return matrix->NumGlobalCols();
3484 return matrix->NumGlobalCols64();
3494 return matrix -> NumMyRows();
3500 std::pair<SparseMatrix::size_type, SparseMatrix::size_type>
3503 size_type begin, end;
3504 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 3505 begin = matrix->RowMap().MinMyGID();
3506 end = matrix->RowMap().MaxMyGID()+1;
3508 begin = matrix->RowMap().MinMyGID64();
3509 end = matrix->RowMap().MaxMyGID64()+1;
3512 return std::make_pair (begin, end);
3521 #ifndef DEAL_II_USE_LARGE_INDEX_TYPE 3522 return matrix->NumGlobalNonzeros();
3524 return matrix->NumGlobalNonzeros64();
3530 template <
typename SparsityType>
3533 const SparsityType &sparsity_pattern,
3534 const MPI_Comm &communicator,
3535 const bool exchange_data)
3538 reinit (map, map, sparsity_pattern, exchange_data);
3543 template <
typename SparsityType>
3546 const IndexSet &col_parallel_partitioning,
3547 const SparsityType &sparsity_pattern,
3548 const MPI_Comm &communicator,
3549 const bool exchange_data)
3551 Epetra_Map row_map =
3553 Epetra_Map col_map =
3555 reinit (row_map, col_map, sparsity_pattern, exchange_data);
3560 template <
typename number>
3563 const ::SparseMatrix<number> &sparse_matrix,
3564 const MPI_Comm &communicator,
3565 const double drop_tolerance,
3566 const bool copy_values,
3567 const ::SparsityPattern *use_this_sparsity)
3570 reinit (map, map, sparse_matrix, drop_tolerance, copy_values,
3576 template <
typename number>
3579 const IndexSet &col_parallel_partitioning,
3580 const ::SparseMatrix<number> &sparse_matrix,
3581 const MPI_Comm &communicator,
3582 const double drop_tolerance,
3583 const bool copy_values,
3584 const ::SparsityPattern *use_this_sparsity)
3586 Epetra_Map row_map =
3588 Epetra_Map col_map =
3590 reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values,
3600 Assert (matrix->Filled(), ExcMatrixNotCompressed());
3601 return matrix->NormOne();
3610 Assert (matrix->Filled(), ExcMatrixNotCompressed());
3611 return matrix->NormInf();
3620 Assert (matrix->Filled(), ExcMatrixNotCompressed());
3621 return matrix->NormFrobenius();
3630 const int ierr = matrix->Scale (a);
3631 Assert (ierr == 0, ExcTrilinosError(ierr));
3645 const TrilinosScalar factor = 1./a;
3647 const int ierr = matrix->Scale (factor);
3648 Assert (ierr == 0, ExcTrilinosError(ierr));
3658 namespace SparseMatrix
3660 template <
typename VectorType>
3662 void check_vector_map_equality(
const Epetra_CrsMatrix &,
3669 void check_vector_map_equality(
const Epetra_CrsMatrix &m,
3674 ExcMessage (
"Column map of matrix does not fit with vector map!"));
3676 ExcMessage (
"Row map of matrix does not fit with vector map!"));
3685 template <
typename VectorType>
3689 const VectorType &src)
const 3691 Assert (&src != &dst, ExcSourceEqualsDestination());
3692 Assert (matrix->Filled(), ExcMatrixNotCompressed());
3696 internal::SparseMatrix::check_vector_map_equality(*matrix, src, dst);
3697 const size_type dst_local_size = dst.end() - dst.begin();
3698 AssertDimension (dst_local_size, static_cast<size_type>(matrix->RangeMap().NumMyElements()));
3699 (void)dst_local_size;
3700 const size_type src_local_size = src.end() - src.begin();
3701 AssertDimension (src_local_size, static_cast<size_type>(matrix->DomainMap().NumMyElements()));
3702 (void)src_local_size;
3704 Epetra_MultiVector tril_dst (View, matrix->RangeMap(), dst.begin(),
3705 matrix->DomainMap().NumMyPoints(), 1);
3706 Epetra_MultiVector tril_src (View, matrix->DomainMap(),
3707 const_cast<TrilinosScalar *
>(src.begin()),
3708 matrix->DomainMap().NumMyPoints(), 1);
3710 const int ierr = matrix->Multiply (
false, tril_src, tril_dst);
3711 Assert (ierr == 0, ExcTrilinosError(ierr));
3717 template <
typename VectorType>
3721 const VectorType &src)
const 3723 Assert (&src != &dst, ExcSourceEqualsDestination());
3724 Assert (matrix->Filled(), ExcMatrixNotCompressed());
3726 internal::SparseMatrix::check_vector_map_equality(*matrix, dst, src);
3727 const size_type dst_local_size = dst.end() - dst.begin();
3728 AssertDimension (dst_local_size, static_cast<size_type>(matrix->DomainMap().NumMyElements()));
3729 const size_type src_local_size = src.end() - src.begin();
3730 AssertDimension (src_local_size, static_cast<size_type>(matrix->RangeMap().NumMyElements()));
3732 Epetra_MultiVector tril_dst (View, matrix->DomainMap(), dst.begin(),
3733 matrix->DomainMap().NumMyPoints(), 1);
3734 Epetra_MultiVector tril_src (View, matrix->RangeMap(),
3735 const_cast<double *
>(src.begin()),
3736 matrix->DomainMap().NumMyPoints(), 1);
3738 const int ierr = matrix->Multiply (
true, tril_src, tril_dst);
3739 Assert (ierr == 0, ExcTrilinosError(ierr));
3745 template <
typename VectorType>
3749 const VectorType &src)
const 3751 Assert (&src != &dst, ExcSourceEqualsDestination());
3759 temp_vector.
reinit(dst.end()-dst.begin(),
true);
3762 vmult (temp_vector,
static_cast<const ::Vector<TrilinosScalar>&
>(src_view));
3763 if (dst_view.size() > 0)
3764 dst_view += temp_vector;
3769 template <
typename VectorType>
3773 const VectorType &src)
const 3775 Assert (&src != &dst, ExcSourceEqualsDestination());
3783 temp_vector.
reinit(dst.end()-dst.begin(),
true);
3786 Tvmult (temp_vector,
static_cast<const ::Vector<TrilinosScalar>&
>(src_view));
3787 if (dst_view.size() > 0)
3788 dst_view += temp_vector;
3797 Assert (row_partitioner().SameAs(domain_partitioner()),
3801 temp_vector.
reinit(v,
true);
3803 vmult (temp_vector, v);
3804 return temp_vector*v;
3814 Assert (row_partitioner().SameAs(domain_partitioner()),
3818 temp_vector.
reinit(v,
true);
3820 vmult (temp_vector, v);
3821 return u*temp_vector;
3841 const Epetra_CrsMatrix &
3844 return static_cast<const Epetra_CrsMatrix &
>(*matrix);
3850 const Epetra_CrsGraph &
3853 return matrix->Graph();
3862 return matrix->DomainMap();
3871 return matrix->RangeMap();
3880 return matrix->RowMap();
3889 return matrix->ColMap();
3917 DEAL_II_NAMESPACE_CLOSE
3920 #endif // DEAL_II_WITH_TRILINOS
const Reference & operator-=(const TrilinosScalar n) const
void reinit(const VectorBase &v, const bool fast=false)
#define DeclException2(Exception2, type1, type2, outsequence)
bool operator!=(const Iterator< Constness > &) const
#define AssertDimension(dim1, dim2)
real_type l2_norm() const
SparseMatrixIterators::Iterator< false > iterator
::types::global_dof_index size_type
void reinit(const size_type n, const bool fast=false)
const Epetra_CrsMatrix & trilinos_matrix() const
void Tvmult(VectorType &dst, const VectorType &src) const
const Reference & operator/=(const TrilinosScalar n) const
Accessor(MatrixType *matrix, const size_type row, const size_type index)
TrilinosScalar matrix_norm_square(const VectorBase &v) const
::ExceptionBase & ExcMessage(std::string arg1)
const Epetra_Map & row_partitioner() const
::types::global_dof_index size_type
const Reference & operator+=(const TrilinosScalar n) const
std_cxx1x::shared_ptr< std::vector< size_type > > colnum_cache
#define AssertThrow(cond, exc)
bool operator>(const Iterator< Constness > &) const
const_iterator begin() const
TrilinosScalar matrix_scalar_product(const VectorBase &u, const VectorBase &v) const
DeclException0(ExcBeyondEndOfMatrix)
void Tvmult_add(VectorType &dst, const VectorType &src) const
SparseMatrix & operator=(const double d)
const Accessor< Constness > * operator->() const
bool is_finite(const double x)
bool is_compressed() const
void reinit(const SparsityType &sparsity_pattern)
const Epetra_Map & col_partitioner() const
TrilinosScalar value() const
void vmult(VectorType &dst, const VectorType &src) const
#define DeclException1(Exception1, type1, outsequence)
void add(const size_type i, const size_type j, const TrilinosScalar value)
#define Assert(cond, exc)
TrilinosScalar linfty_norm() const
std_cxx1x::shared_ptr< std::vector< TrilinosScalar > > value_cache
SparseMatrix & operator/=(const TrilinosScalar factor)
void compress() DEAL_II_DEPRECATED
::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
unsigned int local_size() const
Reference(const Accessor< false > &accessor)
BlockCompressedSparsityPattern CompressedBlockSparsityPattern DEAL_II_DEPRECATED
const Epetra_Map & vector_partitioner() const
const Epetra_Map & domain_partitioner() const
AccessorBase(SparseMatrix *matrix, const size_type row, const size_type index)
Epetra_Map make_trilinos_map(const MPI_Comm &communicator=MPI_COMM_WORLD, const bool overlapping=false) const
const SparseMatrix MatrixType
const Epetra_Map & range_partitioner() const
TrilinosScalar residual(VectorBase &dst, const VectorBase &x, const VectorBase &b) const
SparseMatrix & operator*=(const TrilinosScalar factor)
Accessor< Constness > accessor
::ExceptionBase & ExcIteratorPastEnd()
SparseMatrixIterators::Iterator< true > const_iterator
const Reference & operator=(const TrilinosScalar n) const
::ExceptionBase & ExcNumberNotFinite()
size_type n_nonzero_elements() const
const Reference & operator*=(const TrilinosScalar n) const
DeclException3(ExcAccessToNonlocalRow, std::size_t, std::size_t, std::size_t,<< "You tried to access row "<< arg1<< " of a distributed sparsity pattern, "<< " but only rows "<< arg2<< " through "<< arg3<< " are stored locally and can be accessed.")
Iterator< Constness > & operator++()
const Epetra_CrsGraph & trilinos_sparsity_pattern() const
bool operator<(const Iterator< Constness > &) const
bool operator==(const Iterator< Constness > &) const
::types::global_dof_index size_type
#define DeclException4(Exception4, type1, type2, type3, type4, outsequence)
Accessor< Constness >::MatrixType MatrixType
const_iterator end() const
Accessor(MatrixType *matrix, const size_type row, const size_type index)
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
TrilinosScalar frobenius_norm() const
::ExceptionBase & ExcInternalError()
::ExceptionBase & ExcDivideByZero()
TrilinosScalar value_type
TrilinosScalar l1_norm() const
void set(const size_type i, const size_type j, const TrilinosScalar value)
std::pair< size_type, size_type > local_range() const
const Accessor< Constness > & operator*() const
bool in_local_range(const size_type index) const
void vmult_add(VectorType &dst, const VectorType &src) const