include/deal.II/lac/trilinos_sparse_matrix.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: trilinos_sparse_matrix.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2008, 2009, 2010, 2011, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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                                  // forward declarations
00047 template <typename MatrixType> class BlockMatrixBase;
00048 
00049 template <typename number> class SparseMatrix;
00050 class SparsityPattern;
00051 
00052 namespace TrilinosWrappers
00053 {
00054                                    // forward declarations
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   &parallel_partitioning,
00522                     const unsigned int  n_max_entries_per_row = 0);
00523 
00537       SparseMatrix (const Epetra_Map                &parallel_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    &parallel_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                     &parallel_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     &parallel_partitioning,
00768                     const MPI_Comm     &communicator = MPI_COMM_WORLD,
00769                     const unsigned int  n_max_entries_per_row = 0);
00770 
00785       SparseMatrix (const IndexSet                  &parallel_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      &parallel_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                       &parallel_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 // -------------------------- inline and template functions ----------------------
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                                         // If at end of line: do one
02431                                         // step, then cycle until we
02432                                         // find a row with a nonzero
02433                                         // number of entries.
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                                      // place the iterator on the first entry
02553                                      // past this line, or at the end of the
02554                                      // matrix
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                                      // if there is no such line, then take the
02560                                      // end iterator of the matrix
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                                   // flush buffers
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                                         // Inline the set() and add()
02626                                         // functions, since they will be
02627                                         // called frequently, and the
02628                                         // compiler can optimize away
02629                                         // some unnecessary loops when
02630                                         // the sizes are given at
02631                                         // compile time.
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                                    // If we don't elide zeros, the pointers
02722                                    // are already available...
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                                    // Otherwise, extract nonzero values in
02732                                    // each row and get the respective
02733                                    // indices.
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                                    // If the calling matrix owns the row to
02761                                    // which we want to insert values, we
02762                                    // can directly call the Epetra_CrsMatrix
02763                                    // input function, which is much faster
02764                                    // than the Epetra_FECrsMatrix
02765                                    // function. We distinguish between two
02766                                    // cases: the first one is when the matrix
02767                                    // is not filled (i.e., it is possible to
02768                                    // add new elements to the sparsity pattern),
02769                                    // and the second one is when the pattern is
02770                                    // already fixed. In the former case, we
02771                                    // add the possibility to insert new values,
02772                                    // and in the second we just replace
02773                                    // data.
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                                         // When inserting elements, we do
02783                                         // not want to create exceptions in
02784                                         // the case when inserting non-local
02785                                         // data (since that's what we want
02786                                         // to do right now).
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                                    // When we're at off-processor data, we
02798                                    // have to stick with the standard
02799                                    // Insert/ReplaceGlobalValues
02800                                    // function. Nevertheless, the way we
02801                                    // call it is the fastest one (any other
02802                                    // will lead to repeated allocation and
02803                                    // deallocation of memory in order to
02804                                    // call the function we already use,
02805                                    // which is very unefficient if writing
02806                                    // one element at a time).
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                                   // we have to do checkings on Insert/Add
02842                                   // in any case
02843                                   // to be consistent with the MPI
02844                                   // communication model (see the comments
02845                                   // in the documentation of
02846                                   // TrilinosWrappers::Vector), but we can
02847                                   // save some work if the addend is
02848                                   // zero. However, these actions are done
02849                                   // in case we pass on to the other
02850                                   // function.
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            /*col_indices_are_sorted*/)
02930   {
02931     int ierr;
02932     if (last_action == Insert)
02933       {
02934         // TODO: this could lead to a dead lock when only one processor
02935         // calls GlobalAssemble.
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                                    // If we don't elide zeros, the pointers
02949                                    // are already available...
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                                    // Otherwise, extract nonzero values in
02959                                    // each row and the corresponding index.
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                                    // If the calling matrix owns the row to
02986                                    // which we want to add values, we
02987                                    // can directly call the Epetra_CrsMatrix
02988                                    // input function, which is much faster
02989                                    // than the Epetra_FECrsMatrix function.
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                                    // When we're at off-processor data, we
02999                                    // have to stick with the standard
03000                                    // SumIntoGlobalValues
03001                                    // function. Nevertheless, the way we
03002                                    // call it is the fastest one (any other
03003                                    // will lead to repeated allocation and
03004                                    // deallocation of memory in order to
03005                                    // call the function we already use,
03006                                    // which is very unefficient if writing
03007                                    // one element at a time).
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                                    // inline "simple" functions that are
03045                                    // called frequently and do only involve
03046                                    // a call to some Trilinos function.
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      &parallel_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      &parallel_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                                    // Choose to reinit the vector with fast
03310                                    // argument set, which does not overwrite
03311                                    // the content -- this is what we need
03312                                    // since we're going to overwrite that
03313                                    // anyway in the vmult operation.
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                                    //nothing to do here
03441   }
03442 
03443 
03444 
03445   inline
03446   void
03447   SparseMatrix::prepare_set()
03448   {
03449                                    //nothing to do here
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 /*-----------------------   trilinos_sparse_matrix.h     --------------------*/
03466 
03467 #endif
03468 /*-----------------------   trilinos_sparse_matrix.h     --------------------*/
03469 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 12:03:30 by doxygen 1.7.3