00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__trilinos_sparse_matrix_h
00013 #define __deal2__trilinos_sparse_matrix_h
00014
00015
00016 #include <deal.II/base/config.h>
00017
00018 #ifdef DEAL_II_USE_TRILINOS
00019
00020 # include <deal.II/base/std_cxx1x/shared_ptr.h>
00021 # include <deal.II/base/subscriptor.h>
00022 # include <deal.II/base/index_set.h>
00023 # include <deal.II/lac/full_matrix.h>
00024 # include <deal.II/lac/exceptions.h>
00025 # include <deal.II/lac/trilinos_vector_base.h>
00026 # include <deal.II/lac/parallel_vector.h>
00027
00028 # include <vector>
00029 # include <cmath>
00030 # include <memory>
00031
00032 # define TrilinosScalar double
00033 # include <Epetra_FECrsMatrix.h>
00034 # include <Epetra_Map.h>
00035 # include <Epetra_CrsGraph.h>
00036 # include <Epetra_Vector.h>
00037 # ifdef DEAL_II_COMPILER_SUPPORTS_MPI
00038 # include <Epetra_MpiComm.h>
00039 # include "mpi.h"
00040 # else
00041 # include "Epetra_SerialComm.h"
00042 # endif
00043
00044 DEAL_II_NAMESPACE_OPEN
00045
00046
00047 template <typename MatrixType> class BlockMatrixBase;
00048
00049 template <typename number> class SparseMatrix;
00050 class SparsityPattern;
00051
00052 namespace TrilinosWrappers
00053 {
00054
00055 class VectorBase;
00056 class SparseMatrix;
00057 class SparsityPattern;
00058
00059 namespace MatrixIterators
00060 {
00075 class const_iterator
00076 {
00077 private:
00081 class Accessor
00082 {
00083 public:
00090 Accessor (const SparseMatrix *matrix,
00091 const unsigned int row,
00092 const unsigned int index);
00093
00098 unsigned int row() const;
00099
00104 unsigned int index() const;
00105
00110 unsigned int column() const;
00111
00115 TrilinosScalar value() const;
00116
00120 DeclException0 (ExcBeyondEndOfMatrix);
00121
00125 DeclException3 (ExcAccessToNonlocalRow,
00126 int, int, int,
00127 << "You tried to access row " << arg1
00128 << " of a distributed matrix, but only rows "
00129 << arg2 << " through " << arg3
00130 << " are stored locally and can be accessed.");
00131
00132 private:
00136 mutable SparseMatrix *matrix;
00137
00141 unsigned int a_row;
00142
00146 unsigned int a_index;
00147
00173 std_cxx1x::shared_ptr<std::vector<unsigned int> > colnum_cache;
00174
00179 std_cxx1x::shared_ptr<std::vector<TrilinosScalar> > value_cache;
00180
00189 void visit_present_row ();
00190
00195 friend class const_iterator;
00196 };
00197
00198 public:
00199
00206 const_iterator (const SparseMatrix *matrix,
00207 const unsigned int row,
00208 const unsigned int index);
00209
00213 const_iterator& operator++ ();
00214
00218 const_iterator operator++ (int);
00219
00223 const Accessor& operator* () const;
00224
00228 const Accessor* operator-> () const;
00229
00235 bool operator == (const const_iterator&) const;
00236
00240 bool operator != (const const_iterator&) const;
00241
00250 bool operator < (const const_iterator&) const;
00251
00255 DeclException2 (ExcInvalidIndexWithinRow,
00256 int, int,
00257 << "Attempt to access element " << arg2
00258 << " of row " << arg1
00259 << " which doesn't have that many elements.");
00260
00261 private:
00266 Accessor accessor;
00267 };
00268
00269 }
00270
00271
00302 class SparseMatrix : public Subscriptor
00303 {
00304 public:
00317 struct Traits
00318 {
00324 static const bool zero_addition_can_be_elided = true;
00325 };
00326
00331 typedef MatrixIterators::const_iterator const_iterator;
00332
00338 typedef TrilinosScalar value_type;
00339
00348 SparseMatrix ();
00349
00359 SparseMatrix (const unsigned int m,
00360 const unsigned int n,
00361 const unsigned int n_max_entries_per_row);
00362
00373 SparseMatrix (const unsigned int m,
00374 const unsigned int n,
00375 const std::vector<unsigned int> &n_entries_per_row);
00376
00381 SparseMatrix (const SparsityPattern &InputSparsityPattern);
00382
00390 SparseMatrix (const SparseMatrix &InputMatrix);
00391
00397 virtual ~SparseMatrix ();
00398
00427 template<typename SparsityType>
00428 void reinit (const SparsityType &sparsity_pattern);
00429
00441 void reinit (const SparsityPattern &sparsity_pattern);
00442
00453 void reinit (const SparseMatrix &sparse_matrix);
00454
00477 template <typename number>
00478 void reinit (const ::SparseMatrix<number> &dealii_sparse_matrix,
00479 const double drop_tolerance=1e-13,
00480 const bool copy_values=true,
00481 const ::SparsityPattern *use_this_sparsity=0);
00482
00490 void reinit (const Epetra_CrsMatrix &input_matrix,
00491 const bool copy_values = true);
00493
00521 SparseMatrix (const Epetra_Map ¶llel_partitioning,
00522 const unsigned int n_max_entries_per_row = 0);
00523
00537 SparseMatrix (const Epetra_Map ¶llel_partitioning,
00538 const std::vector<unsigned int> &n_entries_per_row);
00539
00569 SparseMatrix (const Epetra_Map &row_parallel_partitioning,
00570 const Epetra_Map &col_parallel_partitioning,
00571 const unsigned int n_max_entries_per_row = 0);
00572
00602 SparseMatrix (const Epetra_Map &row_parallel_partitioning,
00603 const Epetra_Map &col_parallel_partitioning,
00604 const std::vector<unsigned int> &n_entries_per_row);
00605
00643 template<typename SparsityType>
00644 void reinit (const Epetra_Map ¶llel_partitioning,
00645 const SparsityType &sparsity_pattern,
00646 const bool exchange_data = false);
00647
00669 template<typename SparsityType>
00670 void reinit (const Epetra_Map &row_parallel_partitioning,
00671 const Epetra_Map &col_parallel_partitioning,
00672 const SparsityType &sparsity_pattern,
00673 const bool exchange_data = false);
00674
00703 template <typename number>
00704 void reinit (const Epetra_Map ¶llel_partitioning,
00705 const ::SparseMatrix<number> &dealii_sparse_matrix,
00706 const double drop_tolerance=1e-13,
00707 const bool copy_values=true,
00708 const ::SparsityPattern *use_this_sparsity=0);
00709
00731 template <typename number>
00732 void reinit (const Epetra_Map &row_parallel_partitioning,
00733 const Epetra_Map &col_parallel_partitioning,
00734 const ::SparseMatrix<number> &dealii_sparse_matrix,
00735 const double drop_tolerance=1e-13,
00736 const bool copy_values=true,
00737 const ::SparsityPattern *use_this_sparsity=0);
00739
00767 SparseMatrix (const IndexSet ¶llel_partitioning,
00768 const MPI_Comm &communicator = MPI_COMM_WORLD,
00769 const unsigned int n_max_entries_per_row = 0);
00770
00785 SparseMatrix (const IndexSet ¶llel_partitioning,
00786 const MPI_Comm &communicator,
00787 const std::vector<unsigned int> &n_entries_per_row);
00788
00816 SparseMatrix (const IndexSet &row_parallel_partitioning,
00817 const IndexSet &col_parallel_partitioning,
00818 const MPI_Comm &communicator = MPI_COMM_WORLD,
00819 const unsigned int n_max_entries_per_row = 0);
00820
00850 SparseMatrix (const IndexSet &row_parallel_partitioning,
00851 const IndexSet &col_parallel_partitioning,
00852 const MPI_Comm &communicator,
00853 const std::vector<unsigned int> &n_entries_per_row);
00854
00894 template<typename SparsityType>
00895 void reinit (const IndexSet ¶llel_partitioning,
00896 const SparsityType &sparsity_pattern,
00897 const MPI_Comm &communicator = MPI_COMM_WORLD,
00898 const bool exchange_data = false);
00899
00921 template<typename SparsityType>
00922 void reinit (const IndexSet &row_parallel_partitioning,
00923 const IndexSet &col_parallel_partitioning,
00924 const SparsityType &sparsity_pattern,
00925 const MPI_Comm &communicator = MPI_COMM_WORLD,
00926 const bool exchange_data = false);
00927
00956 template <typename number>
00957 void reinit (const IndexSet ¶llel_partitioning,
00958 const ::SparseMatrix<number> &dealii_sparse_matrix,
00959 const MPI_Comm &communicator = MPI_COMM_WORLD,
00960 const double drop_tolerance=1e-13,
00961 const bool copy_values=true,
00962 const ::SparsityPattern *use_this_sparsity=0);
00963
00985 template <typename number>
00986 void reinit (const IndexSet &row_parallel_partitioning,
00987 const IndexSet &col_parallel_partitioning,
00988 const ::SparseMatrix<number> &dealii_sparse_matrix,
00989 const MPI_Comm &communicator = MPI_COMM_WORLD,
00990 const double drop_tolerance=1e-13,
00991 const bool copy_values=true,
00992 const ::SparsityPattern *use_this_sparsity=0);
00994
00998
01003 unsigned int m () const;
01004
01009 unsigned int n () const;
01010
01025 unsigned int local_size () const;
01026
01043 std::pair<unsigned int, unsigned int>
01044 local_range () const;
01045
01051 bool in_local_range (const unsigned int index) const;
01052
01057 unsigned int n_nonzero_elements () const;
01058
01063 unsigned int row_length (const unsigned int row) const;
01064
01074 bool is_compressed () const;
01075
01084 std::size_t memory_consumption () const;
01085
01087
01091
01108 SparseMatrix &
01109 operator = (const double d);
01110
01121 void clear ();
01122
01166 void compress ();
01167
01194 void set (const unsigned int i,
01195 const unsigned int j,
01196 const TrilinosScalar value);
01197
01234 void set (const std::vector<unsigned int> &indices,
01235 const FullMatrix<TrilinosScalar> &full_matrix,
01236 const bool elide_zero_values = false);
01237
01245 void set (const std::vector<unsigned int> &row_indices,
01246 const std::vector<unsigned int> &col_indices,
01247 const FullMatrix<TrilinosScalar> &full_matrix,
01248 const bool elide_zero_values = false);
01249
01277 void set (const unsigned int row,
01278 const std::vector<unsigned int> &col_indices,
01279 const std::vector<TrilinosScalar> &values,
01280 const bool elide_zero_values = false);
01281
01309 void set (const unsigned int row,
01310 const unsigned int n_cols,
01311 const unsigned int *col_indices,
01312 const TrilinosScalar *values,
01313 const bool elide_zero_values = false);
01314
01330 void add (const unsigned int i,
01331 const unsigned int j,
01332 const TrilinosScalar value);
01333
01370 void add (const std::vector<unsigned int> &indices,
01371 const FullMatrix<TrilinosScalar> &full_matrix,
01372 const bool elide_zero_values = true);
01373
01381 void add (const std::vector<unsigned int> &row_indices,
01382 const std::vector<unsigned int> &col_indices,
01383 const FullMatrix<TrilinosScalar> &full_matrix,
01384 const bool elide_zero_values = true);
01385
01412 void add (const unsigned int row,
01413 const std::vector<unsigned int> &col_indices,
01414 const std::vector<TrilinosScalar> &values,
01415 const bool elide_zero_values = true);
01416
01442 void add (const unsigned int row,
01443 const unsigned int n_cols,
01444 const unsigned int *col_indices,
01445 const TrilinosScalar *values,
01446 const bool elide_zero_values = true,
01447 const bool col_indices_are_sorted = false);
01448
01453 SparseMatrix & operator *= (const TrilinosScalar factor);
01454
01459 SparseMatrix & operator /= (const TrilinosScalar factor);
01460
01465 void copy_from (const SparseMatrix &source);
01466
01479 void add (const TrilinosScalar factor,
01480 const SparseMatrix &matrix);
01481
01521 void clear_row (const unsigned int row,
01522 const TrilinosScalar new_diag_value = 0);
01523
01540 void clear_rows (const std::vector<unsigned int> &rows,
01541 const TrilinosScalar new_diag_value = 0);
01542
01547 void transpose ();
01548
01550
01554
01574 TrilinosScalar operator () (const unsigned int i,
01575 const unsigned int j) const;
01576
01614 TrilinosScalar el (const unsigned int i,
01615 const unsigned int j) const;
01616
01628 TrilinosScalar diag_element (const unsigned int i) const;
01629
01631
01635
01664 void vmult (VectorBase &dst,
01665 const VectorBase &src) const;
01666
01672 void vmult (parallel::distributed::Vector<TrilinosScalar> &dst,
01673 const parallel::distributed::Vector<TrilinosScalar> &src) const;
01674
01707 void Tvmult (VectorBase &dst,
01708 const VectorBase &src) const;
01709
01715 void Tvmult (parallel::distributed::Vector<TrilinosScalar> &dst,
01716 const parallel::distributed::Vector<TrilinosScalar> &src) const;
01717
01748 void vmult_add (VectorBase &dst,
01749 const VectorBase &src) const;
01750
01784 void Tvmult_add (VectorBase &dst,
01785 const VectorBase &src) const;
01786
01836 TrilinosScalar matrix_norm_square (const VectorBase &v) const;
01837
01873 TrilinosScalar matrix_scalar_product (const VectorBase &u,
01874 const VectorBase &v) const;
01875
01910 TrilinosScalar residual (VectorBase &dst,
01911 const VectorBase &x,
01912 const VectorBase &b) const;
01913
01938 void mmult (SparseMatrix &C,
01939 const SparseMatrix &B,
01940 const VectorBase &V = VectorBase()) const;
01941
01942
01969 void Tmmult (SparseMatrix &C,
01970 const SparseMatrix &B,
01971 const VectorBase &V = VectorBase()) const;
01972
01974
01978
01995 TrilinosScalar l1_norm () const;
01996
02012 TrilinosScalar linfty_norm () const;
02013
02021 TrilinosScalar frobenius_norm () const;
02022
02024
02028
02034 const Epetra_CrsMatrix & trilinos_matrix () const;
02035
02043 const Epetra_CrsGraph & trilinos_sparsity_pattern () const;
02044
02054 const Epetra_Map & domain_partitioner () const;
02055
02065 const Epetra_Map & range_partitioner () const;
02066
02074 const Epetra_Map & row_partitioner () const;
02075
02085 const Epetra_Map & col_partitioner () const;
02087
02091
02096 const_iterator begin () const;
02097
02101 const_iterator end () const;
02102
02116 const_iterator begin (const unsigned int r) const;
02117
02133 const_iterator end (const unsigned int r) const;
02134
02136
02140
02150 void write_ascii ();
02151
02164 void print (std::ostream &out,
02165 const bool write_extended_trilinos_info = false) const;
02166
02168
02175 DeclException1 (ExcTrilinosError,
02176 int,
02177 << "An error with error number " << arg1
02178 << " occured while calling a Trilinos function");
02179
02183 DeclException2 (ExcInvalidIndex,
02184 int, int,
02185 << "The entry with index <" << arg1 << ',' << arg2
02186 << "> does not exist.");
02187
02191 DeclException0 (ExcSourceEqualsDestination);
02192
02196 DeclException0 (ExcMatrixNotCompressed);
02197
02201 DeclException4 (ExcAccessToNonLocalElement,
02202 int, int, int, int,
02203 << "You tried to access element (" << arg1
02204 << "/" << arg2 << ")"
02205 << " of a distributed matrix, but only rows "
02206 << arg3 << " through " << arg4
02207 << " are stored locally and can be accessed.");
02208
02212 DeclException2 (ExcAccessToNonPresentElement,
02213 int, int,
02214 << "You tried to access element (" << arg1
02215 << "/" << arg2 << ")"
02216 << " of a sparse matrix, but it appears to not"
02217 << " exist in the Trilinos sparsity pattern.");
02219
02220
02221
02222 protected:
02223
02246 void prepare_add();
02247
02255 void prepare_set();
02256
02257
02258
02259 private:
02260
02268 std_cxx1x::shared_ptr<Epetra_Map> column_space_map;
02269
02280 std_cxx1x::shared_ptr<Epetra_FECrsMatrix> matrix;
02281
02303 Epetra_CombineMode last_action;
02304
02310 bool compressed;
02311
02319 mutable VectorBase temp_vector;
02320
02328 std::vector<unsigned int> column_indices;
02329
02337 std::vector<TrilinosScalar> column_values;
02338
02344 friend class BlockMatrixBase<SparseMatrix>;
02345 };
02346
02347
02348
02349
02350
02351
02352 #ifndef DOXYGEN
02353
02354 namespace MatrixIterators
02355 {
02356
02357 inline
02358 const_iterator::Accessor::
02359 Accessor (const SparseMatrix *matrix,
02360 const unsigned int row,
02361 const unsigned int index)
02362 :
02363 matrix(const_cast<SparseMatrix*>(matrix)),
02364 a_row(row),
02365 a_index(index)
02366 {
02367 visit_present_row ();
02368 }
02369
02370
02371 inline
02372 unsigned int
02373 const_iterator::Accessor::row() const
02374 {
02375 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
02376 return a_row;
02377 }
02378
02379
02380
02381 inline
02382 unsigned int
02383 const_iterator::Accessor::column() const
02384 {
02385 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
02386 return (*colnum_cache)[a_index];
02387 }
02388
02389
02390
02391 inline
02392 unsigned int
02393 const_iterator::Accessor::index() const
02394 {
02395 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
02396 return a_index;
02397 }
02398
02399
02400
02401 inline
02402 TrilinosScalar
02403 const_iterator::Accessor::value() const
02404 {
02405 Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix());
02406 return (*value_cache)[a_index];
02407 }
02408
02409
02410
02411 inline
02412 const_iterator::
02413 const_iterator(const SparseMatrix *matrix,
02414 const unsigned int row,
02415 const unsigned int index)
02416 :
02417 accessor(matrix, row, index)
02418 {}
02419
02420
02421
02422 inline
02423 const_iterator &
02424 const_iterator::operator++ ()
02425 {
02426 Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd());
02427
02428 ++accessor.a_index;
02429
02430
02431
02432
02433
02434 if (accessor.a_index >= accessor.colnum_cache->size())
02435 {
02436 accessor.a_index = 0;
02437 ++accessor.a_row;
02438
02439 while ((accessor.a_row < accessor.matrix->m())
02440 &&
02441 (accessor.matrix->row_length(accessor.a_row) == 0))
02442 ++accessor.a_row;
02443
02444 accessor.visit_present_row();
02445 }
02446 return *this;
02447 }
02448
02449
02450
02451 inline
02452 const_iterator
02453 const_iterator::operator++ (int)
02454 {
02455 const const_iterator old_state = *this;
02456 ++(*this);
02457 return old_state;
02458 }
02459
02460
02461
02462 inline
02463 const const_iterator::Accessor &
02464 const_iterator::operator* () const
02465 {
02466 return accessor;
02467 }
02468
02469
02470
02471 inline
02472 const const_iterator::Accessor *
02473 const_iterator::operator-> () const
02474 {
02475 return &accessor;
02476 }
02477
02478
02479
02480 inline
02481 bool
02482 const_iterator::
02483 operator == (const const_iterator& other) const
02484 {
02485 return (accessor.a_row == other.accessor.a_row &&
02486 accessor.a_index == other.accessor.a_index);
02487 }
02488
02489
02490
02491 inline
02492 bool
02493 const_iterator::
02494 operator != (const const_iterator& other) const
02495 {
02496 return ! (*this == other);
02497 }
02498
02499
02500
02501 inline
02502 bool
02503 const_iterator::
02504 operator < (const const_iterator& other) const
02505 {
02506 return (accessor.row() < other.accessor.row() ||
02507 (accessor.row() == other.accessor.row() &&
02508 accessor.index() < other.accessor.index()));
02509 }
02510
02511 }
02512
02513
02514
02515 inline
02516 SparseMatrix::const_iterator
02517 SparseMatrix::begin() const
02518 {
02519 return const_iterator(this, 0, 0);
02520 }
02521
02522
02523
02524 inline
02525 SparseMatrix::const_iterator
02526 SparseMatrix::end() const
02527 {
02528 return const_iterator(this, m(), 0);
02529 }
02530
02531
02532
02533 inline
02534 SparseMatrix::const_iterator
02535 SparseMatrix::begin(const unsigned int r) const
02536 {
02537 Assert (r < m(), ExcIndexRange(r, 0, m()));
02538 if (row_length(r) > 0)
02539 return const_iterator(this, r, 0);
02540 else
02541 return end (r);
02542 }
02543
02544
02545
02546 inline
02547 SparseMatrix::const_iterator
02548 SparseMatrix::end(const unsigned int r) const
02549 {
02550 Assert (r < m(), ExcIndexRange(r, 0, m()));
02551
02552
02553
02554
02555 for (unsigned int i=r+1; i<m(); ++i)
02556 if (row_length(i) > 0)
02557 return const_iterator(this, i, 0);
02558
02559
02560
02561 return end();
02562 }
02563
02564
02565
02566 inline
02567 bool
02568 SparseMatrix::in_local_range (const unsigned int index) const
02569 {
02570 int begin, end;
02571 begin = matrix->RowMap().MinMyGID();
02572 end = matrix->RowMap().MaxMyGID()+1;
02573
02574 return ((index >= static_cast<unsigned int>(begin)) &&
02575 (index < static_cast<unsigned int>(end)));
02576 }
02577
02578
02579
02580 inline
02581 bool
02582 SparseMatrix::is_compressed () const
02583 {
02584 return compressed;
02585 }
02586
02587
02588
02589 inline
02590 void
02591 SparseMatrix::compress ()
02592 {
02593
02594 int ierr;
02595 ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(),
02596 true);
02597
02598 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
02599
02600 ierr = matrix->OptimizeStorage ();
02601 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
02602
02603 last_action = Zero;
02604
02605 compressed = true;
02606 }
02607
02608
02609
02610 inline
02611 SparseMatrix &
02612 SparseMatrix::operator = (const double d)
02613 {
02614 Assert (d==0, ExcScalarAssignmentOnlyForZeroValue());
02615 compress ();
02616
02617 const int ierr = matrix->PutScalar(d);
02618 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
02619
02620 return *this;
02621 }
02622
02623
02624
02625
02626
02627
02628
02629
02630
02631
02632 inline
02633 void
02634 SparseMatrix::set (const unsigned int i,
02635 const unsigned int j,
02636 const TrilinosScalar value)
02637 {
02638
02639 Assert (numbers::is_finite(value), ExcNumberNotFinite());
02640
02641 set (i, 1, &j, &value, false);
02642 }
02643
02644
02645
02646 inline
02647 void
02648 SparseMatrix::set (const std::vector<unsigned int> &indices,
02649 const FullMatrix<TrilinosScalar> &values,
02650 const bool elide_zero_values)
02651 {
02652 Assert (indices.size() == values.m(),
02653 ExcDimensionMismatch(indices.size(), values.m()));
02654 Assert (values.m() == values.n(), ExcNotQuadratic());
02655
02656 for (unsigned int i=0; i<indices.size(); ++i)
02657 set (indices[i], indices.size(), &indices[0], &values(i,0),
02658 elide_zero_values);
02659 }
02660
02661
02662
02663 inline
02664 void
02665 SparseMatrix::set (const std::vector<unsigned int> &row_indices,
02666 const std::vector<unsigned int> &col_indices,
02667 const FullMatrix<TrilinosScalar> &values,
02668 const bool elide_zero_values)
02669 {
02670 Assert (row_indices.size() == values.m(),
02671 ExcDimensionMismatch(row_indices.size(), values.m()));
02672 Assert (col_indices.size() == values.n(),
02673 ExcDimensionMismatch(col_indices.size(), values.n()));
02674
02675 for (unsigned int i=0; i<row_indices.size(); ++i)
02676 set (row_indices[i], col_indices.size(), &col_indices[0], &values(i,0),
02677 elide_zero_values);
02678 }
02679
02680
02681
02682 inline
02683 void
02684 SparseMatrix::set (const unsigned int row,
02685 const std::vector<unsigned int> &col_indices,
02686 const std::vector<TrilinosScalar> &values,
02687 const bool elide_zero_values)
02688 {
02689 Assert (col_indices.size() == values.size(),
02690 ExcDimensionMismatch(col_indices.size(), values.size()));
02691
02692 set (row, col_indices.size(), &col_indices[0], &values[0],
02693 elide_zero_values);
02694 }
02695
02696
02697
02698 inline
02699 void
02700 SparseMatrix::set (const unsigned int row,
02701 const unsigned int n_cols,
02702 const unsigned int *col_indices,
02703 const TrilinosScalar *values,
02704 const bool elide_zero_values)
02705 {
02706 int ierr;
02707 if (last_action == Add)
02708 {
02709 ierr = matrix->GlobalAssemble (*column_space_map, matrix->RowMap(),
02710 true);
02711
02712 Assert (ierr == 0, ExcTrilinosError(ierr));
02713 }
02714
02715 last_action = Insert;
02716
02717 int * col_index_ptr;
02718 TrilinosScalar const* col_value_ptr;
02719 int n_columns;
02720
02721
02722
02723 if (elide_zero_values == false)
02724 {
02725 col_index_ptr = (int*)col_indices;
02726 col_value_ptr = values;
02727 n_columns = n_cols;
02728 }
02729 else
02730 {
02731
02732
02733
02734 if (column_indices.size() < n_cols)
02735 {
02736 column_indices.resize(n_cols);
02737 column_values.resize(n_cols);
02738 }
02739
02740 n_columns = 0;
02741 for (unsigned int j=0; j<n_cols; ++j)
02742 {
02743 const double value = values[j];
02744 Assert (numbers::is_finite(value), ExcNumberNotFinite());
02745 if (value != 0)
02746 {
02747 column_indices[n_columns] = col_indices[j];
02748 column_values[n_columns] = value;
02749 n_columns++;
02750 }
02751 }
02752
02753 Assert(n_columns <= (int)n_cols, ExcInternalError());
02754
02755 col_index_ptr = (int*)&column_indices[0];
02756 col_value_ptr = &column_values[0];
02757 }
02758
02759
02760
02761
02762
02763
02764
02765
02766
02767
02768
02769
02770
02771
02772
02773
02774 if (row_partitioner().MyGID(row) == true)
02775 {
02776 if (matrix->Filled() == false)
02777 {
02778 ierr = matrix->Epetra_CrsMatrix::InsertGlobalValues(row, n_columns,
02779 const_cast<double*>(col_value_ptr),
02780 col_index_ptr);
02781
02782
02783
02784
02785
02786
02787 if (ierr > 0)
02788 ierr = 0;
02789 }
02790 else
02791 ierr = matrix->Epetra_CrsMatrix::ReplaceGlobalValues(row, n_columns,
02792 const_cast<double*>(col_value_ptr),
02793 col_index_ptr);
02794 }
02795 else
02796 {
02797
02798
02799
02800
02801
02802
02803
02804
02805
02806
02807 compressed = false;
02808
02809 if (matrix->Filled() == false)
02810 {
02811 ierr = matrix->InsertGlobalValues (1, (int*)&row,
02812 n_columns, col_index_ptr,
02813 &col_value_ptr,
02814 Epetra_FECrsMatrix::ROW_MAJOR);
02815 if (ierr > 0)
02816 ierr = 0;
02817 }
02818 else
02819 ierr = matrix->ReplaceGlobalValues (1, (int*)&row,
02820 n_columns, col_index_ptr,
02821 &col_value_ptr,
02822 Epetra_FECrsMatrix::ROW_MAJOR);
02823 }
02824
02825 Assert (ierr <= 0, ExcAccessToNonPresentElement(row, col_index_ptr[0]));
02826 AssertThrow (ierr >= 0, ExcTrilinosError(ierr));
02827 }
02828
02829
02830
02831 inline
02832 void
02833 SparseMatrix::add (const unsigned int i,
02834 const unsigned int j,
02835 const TrilinosScalar value)
02836 {
02837 Assert (numbers::is_finite(value), ExcNumberNotFinite());
02838
02839 if (value == 0)
02840 {
02841
02842
02843
02844
02845
02846
02847
02848
02849
02850
02851 if (last_action == Insert)
02852 {
02853 int ierr;
02854 ierr = matrix->GlobalAssemble(*column_space_map,
02855 row_partitioner(), false);
02856
02857 Assert (ierr == 0, ExcTrilinosError(ierr));
02858 }
02859
02860 last_action = Add;
02861
02862 return;
02863 }
02864 else
02865 add (i, 1, &j, &value, false);
02866 }
02867
02868
02869
02870 inline
02871 void
02872 SparseMatrix::add (const std::vector<unsigned int> &indices,
02873 const FullMatrix<TrilinosScalar> &values,
02874 const bool elide_zero_values)
02875 {
02876 Assert (indices.size() == values.m(),
02877 ExcDimensionMismatch(indices.size(), values.m()));
02878 Assert (values.m() == values.n(), ExcNotQuadratic());
02879
02880 for (unsigned int i=0; i<indices.size(); ++i)
02881 add (indices[i], indices.size(), &indices[0], &values(i,0),
02882 elide_zero_values);
02883 }
02884
02885
02886
02887 inline
02888 void
02889 SparseMatrix::add (const std::vector<unsigned int> &row_indices,
02890 const std::vector<unsigned int> &col_indices,
02891 const FullMatrix<TrilinosScalar> &values,
02892 const bool elide_zero_values)
02893 {
02894 Assert (row_indices.size() == values.m(),
02895 ExcDimensionMismatch(row_indices.size(), values.m()));
02896 Assert (col_indices.size() == values.n(),
02897 ExcDimensionMismatch(col_indices.size(), values.n()));
02898
02899 for (unsigned int i=0; i<row_indices.size(); ++i)
02900 add (row_indices[i], col_indices.size(), &col_indices[0],
02901 &values(i,0), elide_zero_values);
02902 }
02903
02904
02905
02906 inline
02907 void
02908 SparseMatrix::add (const unsigned int row,
02909 const std::vector<unsigned int> &col_indices,
02910 const std::vector<TrilinosScalar> &values,
02911 const bool elide_zero_values)
02912 {
02913 Assert (col_indices.size() == values.size(),
02914 ExcDimensionMismatch(col_indices.size(), values.size()));
02915
02916 add (row, col_indices.size(), &col_indices[0], &values[0],
02917 elide_zero_values);
02918 }
02919
02920
02921
02922 inline
02923 void
02924 SparseMatrix::add (const unsigned int row,
02925 const unsigned int n_cols,
02926 const unsigned int *col_indices,
02927 const TrilinosScalar *values,
02928 const bool elide_zero_values,
02929 const bool )
02930 {
02931 int ierr;
02932 if (last_action == Insert)
02933 {
02934
02935
02936 ierr = matrix->GlobalAssemble(*column_space_map,
02937 row_partitioner(), false);
02938
02939 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
02940 }
02941
02942 last_action = Add;
02943
02944 int * col_index_ptr;
02945 TrilinosScalar const* col_value_ptr;
02946 int n_columns;
02947
02948
02949
02950 if (elide_zero_values == false)
02951 {
02952 col_index_ptr = (int*)col_indices;
02953 col_value_ptr = values;
02954 n_columns = n_cols;
02955 }
02956 else
02957 {
02958
02959
02960 if (column_indices.size() < n_cols)
02961 {
02962 column_indices.resize(n_cols);
02963 column_values.resize(n_cols);
02964 }
02965
02966 n_columns = 0;
02967 for (unsigned int j=0; j<n_cols; ++j)
02968 {
02969 const double value = values[j];
02970 Assert (numbers::is_finite(value), ExcNumberNotFinite());
02971 if (value != 0)
02972 {
02973 column_indices[n_columns] = col_indices[j];
02974 column_values[n_columns] = value;
02975 n_columns++;
02976 }
02977 }
02978
02979 Assert(n_columns <= (int)n_cols, ExcInternalError());
02980
02981 col_index_ptr = (int*)&column_indices[0];
02982 col_value_ptr = &column_values[0];
02983 }
02984
02985
02986
02987
02988
02989
02990 if (row_partitioner().MyGID(row) == true)
02991 {
02992 ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues(row, n_columns,
02993 const_cast<double*>(col_value_ptr),
02994 col_index_ptr);
02995 }
02996 else
02997 {
02998
02999
03000
03001
03002
03003
03004
03005
03006
03007
03008 compressed = false;
03009
03010 ierr = matrix->SumIntoGlobalValues (1, (int*)&row, n_columns,
03011 col_index_ptr,
03012 &col_value_ptr,
03013 Epetra_FECrsMatrix::ROW_MAJOR);
03014 }
03015
03016 #ifdef DEBUG
03017 if (ierr > 0)
03018 {
03019 std::cout << "------------------------------------------"
03020 << std::endl;
03021 std::cout << "Got error " << ierr << " in row " << row
03022 << " of proc " << row_partitioner().Comm().MyPID()
03023 << " when trying to add the columns:" << std::endl;
03024 for (int i=0; i<n_columns; ++i)
03025 std::cout << col_index_ptr[i] << " ";
03026 std::cout << std::endl << std::endl;
03027 std::cout << "Matrix row has the following indices:" << std::endl;
03028 int n_indices, *indices;
03029 trilinos_sparsity_pattern().ExtractMyRowView(row_partitioner().LID(row),
03030 n_indices,
03031 indices);
03032 for (int i=0; i<n_indices; ++i)
03033 std::cout << indices[i] << " ";
03034 std::cout << endl << std::endl;
03035 Assert (ierr <= 0,
03036 ExcAccessToNonPresentElement(row, col_index_ptr[0]));
03037 }
03038 #endif
03039 Assert (ierr >= 0, ExcTrilinosError(ierr));
03040 }
03041
03042
03043
03044
03045
03046
03047 inline
03048 unsigned int
03049 SparseMatrix::m () const
03050 {
03051 return matrix -> NumGlobalRows();
03052 }
03053
03054
03055
03056 inline
03057 unsigned int
03058 SparseMatrix::n () const
03059 {
03060 return matrix -> NumGlobalCols();
03061 }
03062
03063
03064
03065 inline
03066 unsigned int
03067 SparseMatrix::local_size () const
03068 {
03069 return matrix -> NumMyRows();
03070 }
03071
03072
03073
03074 inline
03075 std::pair<unsigned int, unsigned int>
03076 SparseMatrix::local_range () const
03077 {
03078 unsigned int begin, end;
03079 begin = matrix -> RowMap().MinMyGID();
03080 end = matrix -> RowMap().MaxMyGID()+1;
03081
03082 return std::make_pair (begin, end);
03083 }
03084
03085
03086
03087 inline
03088 unsigned int
03089 SparseMatrix::n_nonzero_elements () const
03090 {
03091 return matrix->NumGlobalNonzeros();
03092 }
03093
03094
03095
03096 template <typename SparsityType>
03097 inline
03098 void SparseMatrix::reinit (const IndexSet ¶llel_partitioning,
03099 const SparsityType &sparsity_pattern,
03100 const MPI_Comm &communicator,
03101 const bool exchange_data)
03102 {
03103 Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator, false);
03104 reinit (map, map, sparsity_pattern, exchange_data);
03105 }
03106
03107
03108
03109 template <typename SparsityType>
03110 inline
03111 void SparseMatrix::reinit (const IndexSet &row_parallel_partitioning,
03112 const IndexSet &col_parallel_partitioning,
03113 const SparsityType &sparsity_pattern,
03114 const MPI_Comm &communicator,
03115 const bool exchange_data)
03116 {
03117 Epetra_Map row_map =
03118 row_parallel_partitioning.make_trilinos_map (communicator, false);
03119 Epetra_Map col_map =
03120 col_parallel_partitioning.make_trilinos_map (communicator, false);
03121 reinit (row_map, col_map, sparsity_pattern, exchange_data);
03122 }
03123
03124
03125
03126 template <typename number>
03127 inline
03128 void SparseMatrix::reinit (const IndexSet ¶llel_partitioning,
03129 const ::SparseMatrix<number> &sparse_matrix,
03130 const MPI_Comm &communicator,
03131 const double drop_tolerance,
03132 const bool copy_values,
03133 const ::SparsityPattern *use_this_sparsity)
03134 {
03135 Epetra_Map map = parallel_partitioning.make_trilinos_map (communicator, false);
03136 reinit (map, map, sparse_matrix, drop_tolerance, copy_values,
03137 use_this_sparsity);
03138 }
03139
03140
03141
03142 template <typename number>
03143 inline
03144 void SparseMatrix::reinit (const IndexSet &row_parallel_partitioning,
03145 const IndexSet &col_parallel_partitioning,
03146 const ::SparseMatrix<number> &sparse_matrix,
03147 const MPI_Comm &communicator,
03148 const double drop_tolerance,
03149 const bool copy_values,
03150 const ::SparsityPattern *use_this_sparsity)
03151 {
03152 Epetra_Map row_map =
03153 row_parallel_partitioning.make_trilinos_map (communicator, false);
03154 Epetra_Map col_map =
03155 col_parallel_partitioning.make_trilinos_map (communicator, false);
03156 reinit (row_map, col_map, sparse_matrix, drop_tolerance, copy_values,
03157 use_this_sparsity);
03158 }
03159
03160
03161
03162 inline
03163 TrilinosScalar
03164 SparseMatrix::l1_norm () const
03165 {
03166 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03167 return matrix->NormOne();
03168 }
03169
03170
03171
03172 inline
03173 TrilinosScalar
03174 SparseMatrix::linfty_norm () const
03175 {
03176 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03177 return matrix->NormInf();
03178 }
03179
03180
03181
03182 inline
03183 TrilinosScalar
03184 SparseMatrix::frobenius_norm () const
03185 {
03186 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03187 return matrix->NormFrobenius();
03188 }
03189
03190
03191
03192 inline
03193 SparseMatrix &
03194 SparseMatrix::operator *= (const TrilinosScalar a)
03195 {
03196 const int ierr = matrix->Scale (a);
03197 Assert (ierr == 0, ExcTrilinosError(ierr));
03198
03199 return *this;
03200 }
03201
03202
03203
03204 inline
03205 SparseMatrix &
03206 SparseMatrix::operator /= (const TrilinosScalar a)
03207 {
03208 Assert (a !=0, ExcDivideByZero());
03209
03210 const TrilinosScalar factor = 1./a;
03211
03212 const int ierr = matrix->Scale (factor);
03213 Assert (ierr == 0, ExcTrilinosError(ierr));
03214
03215 return *this;
03216 }
03217
03218
03219
03220 inline
03221 void
03222 SparseMatrix::vmult (VectorBase &dst,
03223 const VectorBase &src) const
03224 {
03225 Assert (&src != &dst, ExcSourceEqualsDestination());
03226 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03227
03228 Assert (src.vector_partitioner().SameAs(matrix->DomainMap()) == true,
03229 ExcMessage ("Column map of matrix does not fit with vector map!"));
03230 Assert (dst.vector_partitioner().SameAs(matrix->RangeMap()) == true,
03231 ExcMessage ("Row map of matrix does not fit with vector map!"));
03232
03233 const int ierr = matrix->Multiply (false, src.trilinos_vector(),
03234 dst.trilinos_vector());
03235 Assert (ierr == 0, ExcTrilinosError(ierr));
03236 }
03237
03238
03239
03240 inline
03241 void
03242 SparseMatrix::vmult (parallel::distributed::Vector<TrilinosScalar> &dst,
03243 const parallel::distributed::Vector<TrilinosScalar> &src) const
03244 {
03245 Assert (&src != &dst, ExcSourceEqualsDestination());
03246 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03247
03248 AssertDimension (dst.local_size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
03249 AssertDimension (src.local_size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
03250
03251 Epetra_Vector tril_dst (View, matrix->RangeMap(), dst.begin());
03252 Epetra_Vector tril_src (View, matrix->DomainMap(),
03253 const_cast<double*>(src.begin()));
03254
03255 const int ierr = matrix->Multiply (false, tril_src, tril_dst);
03256 Assert (ierr == 0, ExcTrilinosError(ierr));
03257 }
03258
03259
03260
03261 inline
03262 void
03263 SparseMatrix::Tvmult (VectorBase &dst,
03264 const VectorBase &src) const
03265 {
03266 Assert (&src != &dst, ExcSourceEqualsDestination());
03267 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03268
03269 Assert (src.vector_partitioner().SameAs(matrix->RangeMap()) == true,
03270 ExcMessage ("Column map of matrix does not fit with vector map!"));
03271 Assert (dst.vector_partitioner().SameAs(matrix->DomainMap()) == true,
03272 ExcMessage ("Row map of matrix does not fit with vector map!"));
03273
03274 const int ierr = matrix->Multiply (true, src.trilinos_vector(),
03275 dst.trilinos_vector());
03276 Assert (ierr == 0, ExcTrilinosError(ierr));
03277 }
03278
03279
03280
03281 inline
03282 void
03283 SparseMatrix::Tvmult (parallel::distributed::Vector<TrilinosScalar> &dst,
03284 const parallel::distributed::Vector<TrilinosScalar>&src) const
03285 {
03286 Assert (&src != &dst, ExcSourceEqualsDestination());
03287 Assert (matrix->Filled(), ExcMatrixNotCompressed());
03288
03289 AssertDimension (dst.local_size(), static_cast<unsigned int>(matrix->DomainMap().NumMyElements()));
03290 AssertDimension (src.local_size(), static_cast<unsigned int>(matrix->RangeMap().NumMyElements()));
03291
03292 Epetra_Vector tril_dst (View, matrix->DomainMap(), dst.begin());
03293 Epetra_Vector tril_src (View, matrix->RangeMap(),
03294 const_cast<double*>(src.begin()));
03295
03296 const int ierr = matrix->Multiply (true, tril_src, tril_dst);
03297 Assert (ierr == 0, ExcTrilinosError(ierr));
03298 }
03299
03300
03301
03302 inline
03303 void
03304 SparseMatrix::vmult_add (VectorBase &dst,
03305 const VectorBase &src) const
03306 {
03307 Assert (&src != &dst, ExcSourceEqualsDestination());
03308
03309
03310
03311
03312
03313
03314 temp_vector.reinit(dst, true);
03315
03316 vmult (temp_vector, src);
03317 dst += temp_vector;
03318 }
03319
03320
03321
03322 inline
03323 void
03324 SparseMatrix::Tvmult_add (VectorBase &dst,
03325 const VectorBase &src) const
03326 {
03327 Assert (&src != &dst, ExcSourceEqualsDestination());
03328
03329 temp_vector.reinit(dst, true);
03330
03331 Tvmult (temp_vector, src);
03332 dst += temp_vector;
03333 }
03334
03335
03336
03337 inline
03338 TrilinosScalar
03339 SparseMatrix::matrix_norm_square (const VectorBase &v) const
03340 {
03341 Assert (row_partitioner().SameAs(domain_partitioner()),
03342 ExcNotQuadratic());
03343
03344 temp_vector.reinit(v);
03345
03346 vmult (temp_vector, v);
03347 return temp_vector*v;
03348 }
03349
03350
03351
03352 inline
03353 TrilinosScalar
03354 SparseMatrix::matrix_scalar_product (const VectorBase &u,
03355 const VectorBase &v) const
03356 {
03357 Assert (row_partitioner().SameAs(domain_partitioner()),
03358 ExcNotQuadratic());
03359
03360 temp_vector.reinit(v);
03361
03362 vmult (temp_vector, v);
03363 return u*temp_vector;
03364 }
03365
03366
03367
03368 inline
03369 TrilinosScalar
03370 SparseMatrix::residual (VectorBase &dst,
03371 const VectorBase &x,
03372 const VectorBase &b) const
03373 {
03374 vmult (dst, x);
03375 dst -= b;
03376 dst *= -1.;
03377
03378 return dst.l2_norm();
03379 }
03380
03381
03382 inline
03383 const Epetra_CrsMatrix &
03384 SparseMatrix::trilinos_matrix () const
03385 {
03386 return static_cast<const Epetra_CrsMatrix&>(*matrix);
03387 }
03388
03389
03390
03391 inline
03392 const Epetra_CrsGraph &
03393 SparseMatrix::trilinos_sparsity_pattern () const
03394 {
03395 return matrix->Graph();
03396 }
03397
03398
03399
03400 inline
03401 const Epetra_Map &
03402 SparseMatrix::domain_partitioner () const
03403 {
03404 return matrix->DomainMap();
03405 }
03406
03407
03408
03409 inline
03410 const Epetra_Map &
03411 SparseMatrix::range_partitioner () const
03412 {
03413 return matrix->RangeMap();
03414 }
03415
03416
03417
03418 inline
03419 const Epetra_Map &
03420 SparseMatrix::row_partitioner () const
03421 {
03422 return matrix->RowMap();
03423 }
03424
03425
03426
03427 inline
03428 const Epetra_Map &
03429 SparseMatrix::col_partitioner () const
03430 {
03431 return matrix->ColMap();
03432 }
03433
03434
03435
03436 inline
03437 void
03438 SparseMatrix::prepare_add()
03439 {
03440
03441 }
03442
03443
03444
03445 inline
03446 void
03447 SparseMatrix::prepare_set()
03448 {
03449
03450 }
03451
03452
03453
03454 #endif // DOXYGEN
03455
03456 }
03457
03458
03459 DEAL_II_NAMESPACE_CLOSE
03460
03461
03462 #endif // DEAL_II_USE_TRILINOS
03463
03464
03465
03466
03467 #endif
03468
03469