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"
49 DEAL_II_NAMESPACE_OPEN
57 namespace TrilinosWrappers
66 namespace SparseMatrixIterators
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.");
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>
413 template <
bool Other>
456 bool operator < (const Iterator<Constness> &)
const;
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;
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>
1127 const SparsityType &sparsity_pattern,
1128 const MPI_Comm &communicator = MPI_COMM_WORLD,
1129 const bool exchange_data =
false);
1152 template<
typename SparsityType>
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>
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>
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);
1274 std::pair<size_type, size_type>
1397 void compress (::VectorOperation::values operation);
1432 const TrilinosScalar value);
1471 const
FullMatrix<TrilinosScalar> &full_matrix,
1472 const
bool elide_zero_values = false);
1482 const std::vector<
size_type> &col_indices,
1483 const
FullMatrix<TrilinosScalar> &full_matrix,
1484 const
bool elide_zero_values = false);
1514 const std::vector<
size_type> &col_indices,
1515 const std::vector<TrilinosScalar> &values,
1516 const
bool elide_zero_values = false);
1548 const TrilinosScalar *values,
1549 const
bool elide_zero_values = false);
1568 const TrilinosScalar value);
1607 const
FullMatrix<TrilinosScalar> &full_matrix,
1608 const
bool elide_zero_values = true);
1618 const std::vector<
size_type> &col_indices,
1619 const
FullMatrix<TrilinosScalar> &full_matrix,
1620 const
bool elide_zero_values = true);
1649 const std::vector<
size_type> &col_indices,
1650 const std::vector<TrilinosScalar> &values,
1651 const
bool elide_zero_values = true);
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,
1758 const TrilinosScalar new_diag_value = 0);
1777 const TrilinosScalar new_diag_value = 0);
1810 TrilinosScalar operator () (const
size_type i,
2185 TrilinosScalar
l1_norm () const;
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");
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.");
2513 std_cxx1x::shared_ptr<Epetra_FECrsMatrix> matrix;
2560 namespace SparseMatrixIterators
2563 AccessorBase::AccessorBase(SparseMatrix *matrix, size_type row, size_type index)
2569 visit_present_row ();
2574 AccessorBase::size_type
2575 AccessorBase::row()
const
2577 Assert (a_row < matrix->
m(), ExcBeyondEndOfMatrix());
2583 AccessorBase::size_type
2584 AccessorBase::column()
const
2586 Assert (a_row < matrix->
m(), ExcBeyondEndOfMatrix());
2587 return (*colnum_cache)[a_index];
2592 AccessorBase::size_type
2593 AccessorBase::index()
const
2595 Assert (a_row < matrix->
m(), ExcBeyondEndOfMatrix());
2601 Accessor<true>::Accessor (MatrixType *matrix,
2602 const size_type row,
2603 const size_type index)
2605 AccessorBase(const_cast<SparseMatrix *>(matrix), row, index)
2609 template <
bool Other>
2611 Accessor<true>::Accessor(
const Accessor<Other> &other)
2619 Accessor<true>::value()
const
2621 Assert (a_row < matrix->
m(), ExcBeyondEndOfMatrix());
2622 return (*value_cache)[a_index];
2627 Accessor<false>::Reference::Reference (
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];
2641 const Accessor<false>::Reference &
2642 Accessor<false>::Reference::operator = (
const TrilinosScalar
n)
const
2644 (*accessor.value_cache)[accessor.a_index] = n;
2645 accessor.matrix->set(accessor.row(), accessor.column(),
2646 static_cast<TrilinosScalar
>(*this));
2652 const Accessor<false>::Reference &
2653 Accessor<false>::Reference::operator += (
const TrilinosScalar n)
const
2655 (*accessor.value_cache)[accessor.a_index] += n;
2656 accessor.matrix->set(accessor.row(), accessor.column(),
2657 static_cast<TrilinosScalar
>(*this));
2663 const Accessor<false>::Reference &
2664 Accessor<false>::Reference::operator -= (
const TrilinosScalar n)
const
2666 (*accessor.value_cache)[accessor.a_index] -= n;
2667 accessor.matrix->set(accessor.row(), accessor.column(),
2668 static_cast<TrilinosScalar
>(*this));
2674 const Accessor<false>::Reference &
2675 Accessor<false>::Reference::operator *= (
const TrilinosScalar n)
const
2677 (*accessor.value_cache)[accessor.a_index] *= n;
2678 accessor.matrix->set(accessor.row(), accessor.column(),
2679 static_cast<TrilinosScalar
>(*this));
2685 const Accessor<false>::Reference &
2686 Accessor<false>::Reference::operator /= (
const TrilinosScalar n)
const
2688 (*accessor.value_cache)[accessor.a_index] /= n;
2689 accessor.matrix->set(accessor.row(), accessor.column(),
2690 static_cast<TrilinosScalar
>(*this));
2696 Accessor<false>::Accessor (MatrixType *matrix,
2697 const size_type row,
2698 const size_type index)
2700 AccessorBase(matrix, row, index)
2705 Accessor<false>::Reference
2706 Accessor<false>::value()
const
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> &
2736 Iterator<Constness>::operator++ ()
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>
2765 Iterator<Constness>::operator++ (
int)
2767 const Iterator<Constness> old_state = *
this;
2774 template <
bool Constness>
2776 const Accessor<Constness> &
2777 Iterator<Constness>::operator* ()
const
2784 template <
bool Constness>
2786 const Accessor<Constness> *
2787 Iterator<Constness>::operator-> ()
const
2794 template <
bool Constness>
2797 Iterator<Constness>::operator == (
const Iterator<Constness> &other)
const
2799 return (accessor.a_row == other.accessor.a_row &&
2800 accessor.a_index == other.accessor.a_index);
2805 template <
bool Constness>
2808 Iterator<Constness>::operator != (
const Iterator<Constness> &other)
const
2810 return ! (*
this == other);
2815 template <
bool Constness>
2818 Iterator<Constness>::operator < (
const Iterator<Constness> &other)
const
2820 return (accessor.row() < other.accessor.row() ||
2821 (accessor.row() == other.accessor.row() &&
2822 accessor.index() < other.accessor.index()));
2826 template <
bool Constness>
2829 Iterator<Constness>::operator > (
const Iterator<Constness> &other)
const
2831 return (other < *
this);
2878 for (size_type i=r+1; i<
m(); ++i)
2929 for (size_type i=r+1; i<
m(); ++i)
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)));
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,
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,
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
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;
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>
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>
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>
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>
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>
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;
3800 VectorBase temp_vector;
3801 temp_vector.
reinit(v,
true);
3803 vmult (temp_vector, v);
3804 return temp_vector*v;
3812 const VectorBase &v)
const
3817 VectorBase temp_vector;
3818 temp_vector.reinit(v,
true);
3820 vmult (temp_vector, v);
3821 return u*temp_vector;
3829 const VectorBase &x,
3830 const VectorBase &b)
const
3836 return dst.l2_norm();
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
static const bool zero_addition_can_be_elided
DeclException0(ExcSourceEqualsDestination)
bool operator!=(const Iterator< Constness > &) const
#define AssertDimension(dim1, dim2)
SparseMatrixIterators::Iterator< false > iterator
::types::global_dof_index size_type
std_cxx1x::shared_ptr< Epetra_FECrsMatrix > matrix
const Epetra_CrsMatrix & trilinos_matrix() const
std_cxx1x::shared_ptr< Epetra_Map > column_space_map
void Tvmult(VectorType &dst, const VectorType &src) const
Epetra_CombineMode last_action
void clear_rows(const std::vector< size_type > &rows, const TrilinosScalar new_diag_value=0)
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
void Tmmult(SparseMatrix &C, const SparseMatrix &B, const VectorBase &V=VectorBase()) 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
DeclException2(ExcInvalidIndexWithinRow, size_type, size_type,<< "Attempt to access element "<< arg2<< " of row "<< arg1<< " which doesn't have that many elements.")
DeclException4(ExcAccessToNonLocalElement, size_type, size_type, size_type, size_type,<< "You tried to access element ("<< arg1<< "/"<< arg2<< ")"<< " of a distributed matrix, but only rows "<< arg3<< " through "<< arg4<< " are stored locally and can be accessed.")
TrilinosScalar value() const
void vmult(VectorType &dst, const VectorType &src) const
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
virtual void reinit(const size_type N, const bool fast=false)
DeclException1(ExcTrilinosError, int,<< "An error with error number "<< arg1<< " occurred while calling a Trilinos function")
DeclException2(ExcInvalidIndex, size_type, size_type,<< "The entry with index <"<< arg1<< ','<< arg2<< "> does not exist.")
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
void copy_from(const SparseMatrix &source)
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
TrilinosScalar el(const size_type i, const size_type j) const
types::global_dof_index size_type
::ExceptionBase & ExcNumberNotFinite()
size_type n_nonzero_elements() 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
TrilinosScalar diag_element(const size_type i) const
::types::global_dof_index size_type
Accessor< Constness >::MatrixType MatrixType
const_iterator end() const
::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
::ExceptionBase & ExcScalarAssignmentOnlyForZeroValue()
TrilinosScalar frobenius_norm() const
void mmult(SparseMatrix &C, const SparseMatrix &B, const VectorBase &V=VectorBase()) 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
size_type memory_consumption() const
unsigned int row_length(const size_type row) const
void clear_row(const size_type row, const TrilinosScalar new_diag_value=0)
const Accessor< Constness > & operator*() const
bool in_local_range(const size_type index) const
void print(std::ostream &out, const bool write_extended_trilinos_info=false) const
void vmult_add(VectorType &dst, const VectorType &src) const