From: kronbichler Date: Thu, 21 Aug 2008 01:07:49 +0000 (+0000) Subject: Added a rudimentary implementation of Trilinos sparse matrix and vector wrappers... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=184bf0bd0bcaa9576540c42636a7ced475b75d22;p=dealii-svn.git Added a rudimentary implementation of Trilinos sparse matrix and vector wrappers for serial and distributed memory systems. It will still need a lot of testing of the features implemented, and some more comments as well. Build a parallel version of step-31 on top of that. git-svn-id: https://svn.dealii.org/trunk@16617 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h new file mode 100755 index 0000000000..dee47baf38 --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -0,0 +1,1220 @@ +//--------------------------------------------------------------------------- +// $Id: trilinos_sparse_matrix.h 16014 2008-04-30 16:18:03Z young $ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2006, 2007, 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__trilinos_sparse_matrix_h +#define __deal2__trilinos_sparse_matrix_h + + +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#ifdef DEAL_II_USE_TRILINOS + +# define TrilinosScalar double +# include +# include +# include +# ifdef DEAL_II_COMPILER_SUPPORTS_MPI +# include +# include "mpi.h" +# else +# include "Epetra_SerialComm.h" +# endif + +DEAL_II_NAMESPACE_OPEN + + +namespace TrilinosWrappers +{ + // forward declarations + class Vector; + class SparseMatrix; + + namespace MatrixIterators + { +/** + * STL conforming iterator. This class acts as an iterator walking over the + * elements of Trilinos matrices. The implementation of this class is similar + * to the one for PETSc matrices. + * + * Note that Trilinos does not give any guarantees as to the order of elements + * within each row. Note also that accessing the elements of a full matrix + * surprisingly only shows the nonzero elements of the matrix, not all + * elements. + * + * @ingroup TrilinosWrappers + * @author Martin Kronbichler, Wolfgang Bangerth, 2008 + */ + class const_iterator + { + private: + /** + * Accessor class for iterators + */ + class Accessor + { + public: + /** + * Constructor. Since we use + * accessors only for read + * access, a const matrix + * pointer is sufficient. + */ + Accessor (const SparseMatrix *matrix, + const unsigned int row, + const unsigned int index); + + /** + * Row number of the element + * represented by this + * object. + */ + unsigned int row() const; + + /** + * Index in row of the element + * represented by this + * object. + */ + unsigned int index() const; + + /** + * Column number of the + * element represented by + * this object. + */ + unsigned int column() const; + + /** + * Value of this matrix entry. + */ + TrilinosScalar value() const; + + /** + * Exception + */ + DeclException0 (ExcBeyondEndOfMatrix); + /** + * Exception + */ + DeclException3 (ExcAccessToNonlocalRow, + int, int, int, + << "You tried to access row " << arg1 + << " of a distributed matrix, but only rows " + << arg2 << " through " << arg3 + << " are stored locally and can be accessed."); + + private: + /** + * The matrix accessed. + */ + mutable SparseMatrix *matrix; + + /** + * Current row number. + */ + unsigned int a_row; + + /** + * Current index in row. + */ + unsigned int a_index; + + /** + * Cache where we store the + * column indices of the present + * row. This is necessary, since + * Trilinos makes access to the + * elements of its matrices + * rather hard, and it is much + * more efficient to copy all + * column entries of a row once + * when we enter it than + * repeatedly asking Trilinos for + * individual ones. This also + * makes some sense since it is + * likely that we will access + * them sequentially anyway. + * + * In order to make copying of + * iterators/accessor of + * acceptable performance, we + * keep a shared pointer to these + * entries so that more than one + * accessor can access this data + * if necessary. + */ + boost::shared_ptr > colnum_cache; + + /** + * Similar cache for the values + * of this row. + */ + boost::shared_ptr > value_cache; + + /** + * Discard the old row caches + * (they may still be used by + * other accessors) and generate + * new ones for the row pointed + * to presently by this accessor. + */ + void visit_present_row (); + + /** + * Make enclosing class a + * friend. + */ + friend class const_iterator; + }; + + public: + + /** + * Constructor. Create an iterator + * into the matrix @p matrix for the + * given row and the index within it. + */ + const_iterator (const SparseMatrix *matrix, + const unsigned int row, + const unsigned int index); + + /** + * Prefix increment. + */ + const_iterator& operator++ (); + + /** + * Postfix increment. + */ + const_iterator operator++ (int); + + /** + * Dereferencing operator. + */ + const Accessor& operator* () const; + + /** + * Dereferencing operator. + */ + const Accessor* operator-> () const; + + /** + * Comparison. True, if + * both iterators point to + * the same matrix + * position. + */ + bool operator == (const const_iterator&) const; + /** + * Inverse of ==. + */ + bool operator != (const const_iterator&) const; + + /** + * Comparison + * operator. Result is true + * if either the first row + * number is smaller or if + * the row numbers are + * equal and the first + * index is smaller. + */ + bool operator < (const const_iterator&) const; + + /** + * Exception + */ + DeclException2 (ExcInvalidIndexWithinRow, + int, int, + << "Attempt to access element " << arg2 + << " of row " << arg1 + << " which doesn't have that many elements."); + + private: + /** + * Store an object of the + * accessor class. + */ + Accessor accessor; + }; + + } + + +/** + * This class implements a wrapper to use the Trilinos distributed sparse matrix + * class Epetra_FECrsMatrix. This is precisely the kind of matrix we deal with + * all the time - we most likely get it from some assembly process, where also + * entries not locally owned might need to written and hence need to be + * forwarded to the owner. This class is designed to be used in a distributed + * memory architecture with an MPI compiler on the bottom, but works equally + * well also for serial processes. The only requirement for this class to + * work is that Trilinos is installed with the respective compiler as a + * basis. + * + * The interface of this class is modeled after the existing + * SparseMatrix class in deal.II. It has almost the same member + * functions, and is often exchangable. However, since Trilinos only supports a + * single scalar type (double), it is + * not templated, and only works with doubles. + * + * Note that Trilinos only guarantees that operations do what you expect if the + * functions @p GlobalAssemble has been called + * after matrix assembly. Therefore, you need to call + * SparseMatrix::compress() before you actually use the matrix. This also + * calls @p FillComplete that compresses the storage format for sparse + * matrices by discarding unused elements. Trilinos allows to continue with + * assembling the matrix after calls to these functions, but since there are + * no more free entries available after that any more, it is better to only + * call SparseMatrix::compress() once at the end of the assembly stage and + * before the matrix is actively used. + * + * @ingroup TrilinosWrappers + * @ingroup Matrix1 + * @author Martin Kronbichler, Wolfgang Bangerth, 2008 + */ + class SparseMatrix : public Subscriptor + { + public: + /** + * Declare a typedef for the iterator + * class. + */ + typedef MatrixIterators::const_iterator const_iterator; + + /** + * Declare a typedef in analogy to all + * the other container classes. + */ + typedef TrilinosScalar value_type; + + /** + * Default constructor. + */ + SparseMatrix (); + + /** + * Constructor using an Epetra_Map + * and a maximum number of nonzero + * matrix entries. + */ + SparseMatrix (const Epetra_Map &InputMap, + const unsigned int n_max_entries_per_row); + + /** + * Same as before, but now use the + * exact number of nonzeros in each + * matrix row. Since we know the + * number of elements in the matrix + * exactly in this case, we can + * already allocate the right amount + * of memory, which makes the + * creation process including the + * insertion of nonzero elements + * by the respective + * SparseMatrix::reinit call + * considerably faster. + */ + SparseMatrix (const Epetra_Map &InputMap, + const std::vector &n_entries_per_row); + + /** + * Destructor. Made virtual so that one + * can use pointers to this class. + */ + virtual ~SparseMatrix (); + /** + * This operator assigns a scalar to a + * matrix. Since this does usually not + * make much sense (should we set all + * matrix entries to this value? Only + * the nonzero entries of the sparsity + * pattern?), this operation is only + * allowed if the actual value to be + * assigned is zero. This operator only + * exists to allow for the obvious + * notation matrix=0, which + * sets all elements of the matrix to + * zero, but keeps the sparsity pattern + * previously used. + */ + SparseMatrix & + operator = (const double d); + + /** + * This function initializes the + * Trilinos matrix by attaching all + * the elements to the sparsity + * pattern provided as deal argument. + * This function uses a user- + * provided maximum number of + * elements per row. If that is + * not directly available, use one of the + * other reinit functions. + */ + void reinit (const CompressedSparsityPattern &sparsity_pattern, + const unsigned int n_max_entries_per_row); + + /** + * This function initializes the + * Trilinos matrix by attaching all + * the elements to the sparsity + * pattern provided as deal argument, + * now calculating the maximum number + * of nonzeros from the sparsity + * pattern internally. + */ + void reinit (const CompressedSparsityPattern &sparsity_pattern); + + /** + * This function is similar to the + * other initialization function above, + * but now also reassigns the matrix + * rows according to a user-supplied + * Epetra map. This might be used + * when the matrix structure changes, + * e.g. when the grid is refined. + */ + void reinit (const Epetra_Map &input_map, + const CompressedSparsityPattern &sparsity_pattern); + + /** + * Release all memory and return + * to a state just like after + * having called the default + * constructor. + */ + void clear (); + + /** + * Set the element (i,j) + * to @p value. + * + * This function + * adds a new entry to the + * matrix if it didn't exist + * before, very much in + * contrast to the SparseMatrix + * class which throws an error + * if the entry does not exist. + * If value is not a + * finite number an exception + * is thrown. + */ + void set (const unsigned int i, + const unsigned int j, + const TrilinosScalar value); + + /** + * Add @p value to the + * element (i,j). + * + * This function + * adds a new entry to the + * matrix if it didn't exist + * before, very much in + * contrast to the SparseMatrix + * class which throws an error + * if the entry does not exist. + * If value is not a + * finite number an exception + * is thrown. + */ + void add (const unsigned int i, + const unsigned int j, + const TrilinosScalar value); + + /** + * Remove all elements from + * this row by setting + * them to zero. The function + * does not modify the number + * of allocated nonzero + * entries, it only sets some + * entries to zero. It may drop + * them from the sparsity + * pattern, though (but retains + * the allocated memory in case + * new entries are again added + * later). + * + * This operation is used in + * eliminating constraints (e.g. due to + * hanging nodes) and makes sure that + * we can write this modification to + * the matrix without having to read + * entries (such as the locations of + * non-zero elements) from it -- + * without this operation, removing + * constraints on parallel matrices is + * a rather complicated procedure. + * + * The second parameter can be used to + * set the diagonal entry of this row + * to a value different from zero. The + * default is to set it to zero. + */ + // void clear_row (const unsigned int row, + // const TrilinosScalar new_diag_value = 0); + + /** + * Same as clear_row(), except that it + * works on a number of rows at once. + * + * The second parameter can be used to + * set the diagonal entries of all + * cleared rows to something different + * from zero. Note that all of these + * diagonal entries get the same value + * -- if you want different values for + * the diagonal entries, you have to + * set them by hand. + */ + // void clear_rows (const std::vector &rows, + // const TrilinosScalar new_diag_value = 0); + + /** + * Trilinos matrices store their own + * sparsity patterns. So, in analogy to + * our own SparsityPattern class, + * this function compresses the + * sparsity pattern and allows the + * resulting matrix to be used in all + * other operations where before only + * assembly functions were + * allowed. This function must + * therefore be called once you have + * assembled the matrix. + */ + void compress (); + + /** + * Return the value of the entry + * (i,j). This may be an + * expensive operation and you should + * always take care where to call this + * function. In contrast to the + * respective function in the + * @p SparseMatrix class, we don't + * throw an exception if the respective + * entry doesn't exist in the sparsity + * pattern of this class, since Trilinos + * does not transmit this information. + * On the other hand, an exception + * will be thrown when the requested + * element is not saved on the calling + * process. + * + * This function is therefore exactly + * equivalent to the el() function. + */ + TrilinosScalar operator () (const unsigned int i, + const unsigned int j) const; + + /** + * Return the value of the matrix entry + * (i,j). If this entry does not + * exist in the sparsity pattern, then + * zero is returned. While this may be + * convenient in some cases, note that + * it is simple to write algorithms + * that are slow compared to an optimal + * solution, since the sparsity of the + * matrix is not used. + */ + TrilinosScalar el (const unsigned int i, + const unsigned int j) const; + + /** + * Return the main diagonal + * element in the ith + * row. This function throws an + * error if the matrix is not + * quadratic. + * + * TODO: Trilinos can access the + * diagonal faster. Implement this! + */ + TrilinosScalar diag_element (const unsigned int i) const; + + /** + * Return the number of rows in this + * matrix. + */ + unsigned int m () const; + + /** + * Return the number of columns in this + * matrix. + */ + unsigned int n () const; + + /** + * Return the local dimension of the + * matrix, i.e. the number of rows + * stored on the present MPI + * process. For sequential matrices, + * this number is the same as m(), + * but for parallel matrices it may be + * smaller. + * + * To figure out which elements + * exactly are stored locally, + * use local_range(). + */ + unsigned int local_size () const; + + /** + * Return a pair of indices + * indicating which rows of + * this matrix are stored + * locally. The first number is + * the index of the first + * row stored, the second + * the index of the one past + * the last one that is stored + * locally. If this is a + * sequential matrix, then the + * result will be the pair + * (0,m()), otherwise it will be + * a pair (i,i+n), where + * n=local_size(). + */ + std::pair + local_range () const; + + /** + * Return whether @p index is + * in the local range or not, + * see also local_range(). + */ + bool in_local_range (const unsigned int index) const; + + /** + * Return the number of nonzero + * elements of this + * matrix. + */ + unsigned int n_nonzero_elements () const; + + /** + * Number of entries in a specific row. + */ + unsigned int row_length (const unsigned int row) const; + + /** + * Return the l1-norm of the matrix, that is + * $|M|_1=max_{all columns j}\sum_{all + * rows i} |M_ij|$, + * (max. sum of columns). + * This is the + * natural matrix norm that is compatible + * to the l1-norm for vectors, i.e. + * $|Mv|_1\leq |M|_1 |v|_1$. + * (cf. Haemmerlin-Hoffmann: + * Numerische Mathematik) + */ + TrilinosScalar l1_norm () const; + + /** + * Return the linfty-norm of the + * matrix, that is + * $|M|_infty=max_{all rows i}\sum_{all + * columns j} |M_ij|$, + * (max. sum of rows). + * This is the + * natural matrix norm that is compatible + * to the linfty-norm of vectors, i.e. + * $|Mv|_infty \leq |M|_infty |v|_infty$. + * (cf. Haemmerlin-Hoffmann: + * Numerische Mathematik) + */ + TrilinosScalar linfty_norm () const; + + /** + * Return the frobenius norm of the + * matrix, i.e. the square root of the + * sum of squares of all entries in the + * matrix. + */ + TrilinosScalar frobenius_norm () const; + + /** + * Multiply the entire matrix by a + * fixed factor. + */ + SparseMatrix & operator *= (const TrilinosScalar factor); + + /** + * Divide the entire matrix by a + * fixed factor. + */ + SparseMatrix & operator /= (const TrilinosScalar factor); + + /** + * Matrix-vector multiplication: + * let dst = M*src with + * M being this matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + void vmult (Vector &dst, + const Vector &src) const; + + /** + * Matrix-vector multiplication: let + * dst = MT*src with + * M being this matrix. This + * function does the same as vmult() + * but takes the transposed matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + void Tvmult (Vector &dst, + const Vector &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * M*src on dst + * with M being this + * matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + void vmult_add (Vector &dst, + const Vector &src) const; + + /** + * Adding Matrix-vector + * multiplication. Add + * MT*src to + * dst with M being + * this matrix. This function + * does the same as vmult_add() + * but takes the transposed + * matrix. + * + * Source and destination must + * not be the same vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + void Tvmult_add (Vector &dst, + const Vector &src) const; + + /** + * Return the square of the norm + * of the vector $v$ with respect + * to the norm induced by this + * matrix, + * i.e. $\left(v,Mv\right)$. This + * is useful, e.g. in the finite + * element context, where the + * $L_2$ norm of a function + * equals the matrix norm with + * respect to the mass matrix of + * the vector representing the + * nodal values of the finite + * element function. + * + * Obviously, the matrix needs to + * be quadratic for this operation. + * + * The implementation of this function + * is not as efficient as the one in + * the @p SparseMatrix class used in + * deal.II (i.e. the original one, not + * the Trilinos wrapper class) since Trilinos + * doesn't support this operation and + * needs a temporary vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + TrilinosScalar matrix_norm_square (const Vector &v) const; + + /** + * Compute the matrix scalar + * product $\left(u,Mv\right)$. + * + * The implementation of this function + * is not as efficient as the one in + * the @p SparseMatrix class used in + * deal.II (i.e. the original one, not + * the Trilinos wrapper class) since Trilinos + * doesn't support this operation and + * needs a temporary vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + TrilinosScalar matrix_scalar_product (const Vector &u, + const Vector &v) const; + + /** + * Compute the residual of an + * equation Mx=b, where + * the residual is defined to be + * r=b-Mx. Write the + * residual into + * @p dst. The + * l2 norm of + * the residual vector is + * returned. + * + * Source x and destination + * dst must not be the same + * vector. + * + * Note that both vectors have to be + * distributed vectors generated using + * the same Map as was used for the + * matrix in case you work on a + * distributed memory architecture, + * using the interface of the + * TrilinosWrappers::Vector class. + */ + TrilinosScalar residual (Vector &dst, + const Vector &x, + const Vector &b) const; + + /** + * STL-like iterator with the + * first entry. + */ + const_iterator begin () const; + + /** + * Final iterator. + */ + const_iterator end () const; + + /** + * STL-like iterator with the + * first entry of row @p r. + * + * Note that if the given row is empty, + * i.e. does not contain any nonzero + * entries, then the iterator returned by + * this function equals + * end(r). Note also that the + * iterator may not be dereferencable in + * that case. + */ + const_iterator begin (const unsigned int r) const; + + /** + * Final iterator of row r. It + * points to the first element past the + * end of line @p r, or past the end of + * the entire sparsity pattern. + * + * Note that the end iterator is not + * necessarily dereferencable. This is in + * particular the case if it is the end + * iterator for the last row of a matrix. + */ + const_iterator end (const unsigned int r) const; + + /** + * Make an in-place transpose of a + * matrix. + */ + void transpose (); + + /** + * Test whether a matrix is symmetric. + * Default tolerance is zero. + * TODO: Not implemented. + */ + bool is_symmetric (const double tol = 0.0); + + /** + * Test whether a matrix is Hermitian, + * i.e. it is the complex conjugate + * of its transpose. + * TODO: Not implemented. + */ + bool is_hermitian (); + + /** + * Abstract Trilinos object that helps view + * in ASCII other Trilinos objects. Currently + * this function is not implemented. + * TODO: Not implemented. + */ + void write_ascii (); + + // TODO: Write an overloading + // of the operator << for output. + // Since the underlying Trilinos + // object supports it, this should + // be very easy. + + /** + * Exception + */ + DeclException1 (ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occured while calling a Trilinos function"); + + /** + * Exception + */ + DeclException0 (ExcSourceEqualsDestination); + + /** + * Exception + */ + DeclException4 (ExcAccessToNonlocalElement, + int, int, int, int, + << "You tried to access element (" << arg1 + << "/" << arg2 << ")" + << " of a distributed matrix, but only rows " + << arg3 << " through " << arg4 + << " are stored locally and can be accessed."); + + /** + * The Epetra Trilinos mapping that + * assigns parts of the matrix to + * the individual processes. + */ + Epetra_Map map; + + /** + * A sparse matrix object in + * Trilinos to be used for + * finite element based problems + * which allows for assembling into + * non-local elements. + * The actual type, a sparse + * matrix, is set in the constructor. + */ + std::auto_ptr matrix; + + protected: + /** + * Trilinos doesn't allow to mix additions + * to matrix entries and overwriting + * them (to make synchronisation of + * parallel computations + * simpler). The way we do it is to, + * for each access operation, store + * whether it is an insertion or an + * addition. If the previous one was of + * different type, then we first have + * to flush the Trilinos buffers; + * otherwise, we can simply go on. + * Luckily, Trilinos has an object + * for this which does already all + * the parallel communications in + * such a case, so we simply use their + * model, which stores whether the + * last operation was an addition + * or an insertion. + */ + Epetra_CombineMode last_action; + }; + + + +#ifndef DOXYGEN +// -------------------------- inline and template functions ---------------------- + + + namespace MatrixIterators + { + + inline + const_iterator::Accessor:: + Accessor (const SparseMatrix *matrix, + const unsigned int row, + const unsigned int index) + : + matrix(const_cast(matrix)), + a_row(row), + a_index(index) + { + visit_present_row (); + } + + + inline + unsigned int + const_iterator::Accessor::row() const + { + Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return a_row; + } + + + inline + unsigned int + const_iterator::Accessor::column() const + { + Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return (*colnum_cache)[a_index]; + } + + + inline + unsigned int + const_iterator::Accessor::index() const + { + Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return a_index; + } + + + inline + TrilinosScalar + const_iterator::Accessor::value() const + { + Assert (a_row < matrix->m(), ExcBeyondEndOfMatrix()); + return (*value_cache)[a_index]; + } + + + inline + const_iterator:: + const_iterator(const SparseMatrix *matrix, + const unsigned int row, + const unsigned int index) + : + accessor(matrix, row, index) + {} + + + + inline + const_iterator & + const_iterator::operator++ () + { + Assert (accessor.a_row < accessor.matrix->m(), ExcIteratorPastEnd()); + + ++accessor.a_index; + + // if at end of line: do one step, then + // cycle until we find a row with a + // nonzero number of entries + if (accessor.a_index >= accessor.colnum_cache->size()) + { + accessor.a_index = 0; + ++accessor.a_row; + + while ((accessor.a_row < accessor.matrix->m()) + && + (accessor.matrix->row_length(accessor.a_row) == 0)) + ++accessor.a_row; + + accessor.visit_present_row(); + } + return *this; + } + + + inline + const_iterator + const_iterator::operator++ (int) + { + const const_iterator old_state = *this; + ++(*this); + return old_state; + } + + + inline + const const_iterator::Accessor & + const_iterator::operator* () const + { + return accessor; + } + + + inline + const const_iterator::Accessor * + const_iterator::operator-> () const + { + return &accessor; + } + + + inline + bool + const_iterator:: + operator == (const const_iterator& other) const + { + return (accessor.a_row == other.accessor.a_row && + accessor.a_index == other.accessor.a_index); + } + + + inline + bool + const_iterator:: + operator != (const const_iterator& other) const + { + return ! (*this == other); + } + + + inline + bool + const_iterator:: + operator < (const const_iterator& other) const + { + return (accessor.row() < other.accessor.row() || + (accessor.row() == other.accessor.row() && + accessor.index() < other.accessor.index())); + } + + } + + + inline + TrilinosScalar + SparseMatrix::operator() (const unsigned int i, + const unsigned int j) const + { + return el(i,j); + } + + + + inline + SparseMatrix::const_iterator + SparseMatrix::begin() const + { + return const_iterator(this, 0, 0); + } + + + inline + SparseMatrix::const_iterator + SparseMatrix::end() const + { + return const_iterator(this, m(), 0); + } + + + inline + SparseMatrix::const_iterator + SparseMatrix::begin(const unsigned int r) const + { + Assert (r < m(), ExcIndexRange(r, 0, m())); + if (row_length(r) > 0) + return const_iterator(this, r, 0); + else + return end (r); + } + + + inline + SparseMatrix::const_iterator + SparseMatrix::end(const unsigned int r) const + { + Assert (r < m(), ExcIndexRange(r, 0, m())); + + // place the iterator on the first entry + // past this line, or at the end of the + // matrix + for (unsigned int i=r+1; i 0) + return const_iterator(this, i, 0); + + // if there is no such line, then take the + // end iterator of the matrix + return end(); + } + + + + inline + bool + SparseMatrix::in_local_range (const unsigned int index) const + { + int begin, end; + begin = matrix->RowMap().MinMyGID(); + end = matrix->RowMap().MaxMyGID(); + + return ((index >= static_cast(begin)) && + (index < static_cast(end))); + } + +#endif // DOXYGEN +} + + +DEAL_II_NAMESPACE_CLOSE + + +#endif // DEAL_II_USE_TRILINOS + + +/*---------------------------- trilinos_matrix_base.h ---------------------------*/ + +#endif +/*---------------------------- trilinos_matrix_base.h ---------------------------*/ diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h new file mode 100755 index 0000000000..7998b13a13 --- /dev/null +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -0,0 +1,1013 @@ +//--------------------------------------------------------------------------- +// $Id: Trilinos_vector_base.h 15453 2008-01-29 1:00:00Z marquez $ +// Version: $Name$ +// +// Copyright (C) 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- +#ifndef __deal2__trilinos_vector_h +#define __deal2__trilinos_vector_h + + +#include +#include +#include +#include + +#include +#include +#include + +#ifdef DEAL_II_USE_TRILINOS + +#define TrilinosScalar double +# include "Epetra_ConfigDefs.h" +# ifdef DEAL_II_COMPILER_SUPPORTS_MPI // only if MPI is installed +# include "mpi.h" +# include "Epetra_MpiComm.h" +# else +# include "Epetra_SerialComm.h" +# endif +# include "Epetra_FEVector.h" +# include "Epetra_Map.h" + +DEAL_II_NAMESPACE_OPEN + + // forward declaration +template class Vector; + + +/** + * A namespace in which wrapper classes for Trilinos objects reside. + * + * @ingroup TrilinosWrappers + * @ingroup Vectors + * @see @ref SoftwareTrilinos + * @author Martin Kronbichler, Wolfgang Bangerth, 2008 + */ +namespace TrilinosWrappers +{ + // forward declaration + class Vector; + + /** + * @cond internal + */ + +/** + * A namespace for internal implementation details of the TrilinosWrapper + * members. + * @ingroup TrilinosWrappers + */ + namespace internal + { + /** + * This class implements a wrapper for + * accessing the Trilinos vector + * in the same way as we access + * deal.II objects: + * it is initialized with a vector and an + * element within it, and has a + * conversion operator to extract the + * scalar value of this element. It also + * has a variety of assignment operator + * for writing to this one element. + * @ingroup TrilinosWrappers + */ + class VectorReference + { + private: + /** + * Constructor. It is made private so + * as to only allow the actual vector + * class to create it. + */ + VectorReference (Vector &vector, + const unsigned int index); + + public: + /** + * This looks like a copy operator, + * but does something different than + * usual. In particular, it does not + * copy the member variables of this + * reference. Rather, it handles the + * situation where we have two + * vectors @p v and @p w, and assign + * elements like in + * v(i)=w(i). Here, both + * left and right hand side of the + * assignment have data type + * VectorReference, but what we + * really mean is to assign the + * vector elements represented by the + * two references. This operator + * implements this operation. Note + * also that this allows us to make + * the assignment operator const. + */ + const VectorReference & operator = (const VectorReference &r) const; + + /** + * Set the referenced element of the + * vector to s. + */ + const VectorReference & operator = (const TrilinosScalar &s) const; + + /** + * Add s to the referenced + * element of the vector-> + */ + const VectorReference & operator += (const TrilinosScalar &s) const; + + /** + * Subtract s from the + * referenced element of the vector-> + */ + const VectorReference & operator -= (const TrilinosScalar &s) const; + + /** + * Multiply the referenced element of + * the vector by s. + */ + const VectorReference & operator *= (const TrilinosScalar &s) const; + + /** + * Divide the referenced element of + * the vector by s. + */ + const VectorReference & operator /= (const TrilinosScalar &s) const; + + /** + * Convert the reference to an actual + * value, i.e. return the value of + * the referenced element of the + * vector. + */ + operator TrilinosScalar () const; + + /** + * Exception + */ + DeclException1 (ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occured while calling a Trilinos function"); + + /** + * Exception + */ + DeclException3 (ExcAccessToNonlocalElement, + int, int, int, + << "You tried to access element " << arg1 + << " of a distributed vector, but only elements " + << arg2 << " through " << arg3 + << " are stored locally and can be accessed."); + + private: + /** + * Point to the vector we are + * referencing. + */ + Vector &vector; + + /** + * Index of the referenced element + * of the vector. + */ + const unsigned int index; + + /** + * Make the vector class a friend, so + * that it can create objects of the + * present type. + */ + friend class ::dealii::TrilinosWrappers::Vector; + }; + } + /** + * @endcond + */ + + +/** + * This class implements a wrapper to use the Trilinos distributed vector + * class Epetra_FEVector. This is precisely the kind of vector we deal with + * all the time - we probably get it from some assembly process, where also + * entries not locally owned might need to written and hence need to be + * forwarded to the owner. This class is designed to be used in a distributed + * memory architecture with an MPI compiler on the bottom, but works equally + * well also for serial processes. The only requirement for this class to + * work is that Trilinos is installed with the respective compiler as a + * basis. + * + * The interface of this class is modeled after the existing Vector + * class in deal.II. It has almost the same member functions, and is often + * exchangable. However, since Trilinos only supports a single scalar type + * (double), it is not templated, and only works with that type. + * + * Note that Trilinos only guarantees that operations do what you expect if the + * function @p GlobalAssemble has been called after vector assembly in order to + * distribute the data. Therefore, you need to call Vector::compress() + * before you actually use the vectors. + * + *

Parallel communication model

+ * + * The parallel functionality of Trilinos is built on top of the Message Passing + * Interface (MPI). MPI's communication model is built on collective + * communications: if one process wants something from another, that other + * process has to be willing to accept this communication. A process cannot + * query data from another process by calling a remote function, without that + * other process expecting such a transaction. The consequence is that most of + * the operations in the base class of this class have to be called + * collectively. For example, if you want to compute the l2 norm of a parallel + * vector, @em all processes across which this vector is shared have to call + * the @p l2_norm function. If you don't do this, but instead only call the @p + * l2_norm function on one process, then the following happens: This one + * process will call one of the collective MPI functions and wait for all the + * other processes to join in on this. Since the other processes don't call + * this function, you will either get a time-out on the first process, or, + * worse, by the time the next a callto a Trilinos function generates an MPI + * message on the other processes , you will get a cryptic message that only a + * subset of processes attempted a communication. These bugs can be very hard + * to figure out, unless you are well-acquainted with the communication model + * of MPI, and know which functions may generate MPI messages. + * + * One particular case, where an MPI message may be generated unexpectedly is + * discussed below. + * + *

Accessing individual elements of a vector

+ * + * Trilinos does allow read access to individual elements of a vector, but in the + * distributed case only to elements that are stored locally. We implement + * this through calls like d=vec(i). However, if you access an + * element outside the locally stored range, an exception is generated. + * + * In contrast to read access, Trilinos (and the respective deal.II wrapper + * classes) allow to write (or add) to individual elements of vectors, even if + * they are stored on a different process. You can do this writing, for + * example, vec(i)=d or vec(i)+=d, or similar + * operations. There is one catch, however, that may lead to very confusing + * error messages: Trilinos requires application programs to call the compress() + * function when they switch from adding, to elements to writing to + * elements. The reasoning is that all processes might accumulate addition + * operations to elements, even if multiple processes write to the same + * elements. By the time we call compress() the next time, all these additions + * are executed. However, if one process adds to an element, and another + * overwrites to it, the order of execution would yield non-deterministic + * behavior if we don't make sure that a synchronisation with compress() + * happens in between. + * + * In order to make sure these calls to compress() happen at the appropriate + * time, the deal.II wrappers keep a state variable that store which is the + * presently allowed operation: additions or writes. If it encounters an + * operation of the opposite kind, it calls compress() and flips the + * state. This can sometimes lead to very confusing behavior, in code that may + * for example look like this: + * @verbatim + * TrilinosWrappers::Vector vector; + * ... + * // do some write operations on the vector + * for (unsigned int i=0; isize(); ++i) + * vector(i) = i; + * + * // do some additions to vector elements, but + * // only for some elements + * for (unsigned int i=0; isize(); ++i) + * if (some_condition(i) == true) + * vector(i) += 1; + * + * // do another collective operation + * const double norm = vector->l2_norm(); + * @endverbatim + * + * This code can run into trouble: by the time we see the first addition + * operation, we need to flush the overwrite buffers for the vector, and the + * deal.II library will do so by calling compress(). However, it will only do + * so for all processes that actually do an addition -- if the condition is + * never true for one of the processes, then this one will not get to the + * actual compress() call, whereas all the other ones do. This gets us into + * trouble, since all the other processes hang in the call to flush the write + * buffers, while the one other process advances to the call to compute the l2 + * norm. At this time, you will get an error that some operation was attempted + * by only a subset of processes. This behavior may seem surprising, unless + * you know that write/addition operations on single elements may trigger this + * behavior. + * + * The problem described here may be avoided by placing additional calls to + * compress(), or making sure that all processes do the same type of + * operations at the same time, for example by placing zero additions if + * necessary. + * + * @ingroup TrilinosWrappers + * @author Martin Kronbichler, Wolfgang Bangerth, 2008 + */ + class Vector + { + public: + /** + * Declare some of the standard types + * used in all containers. These types + * parallel those in the C + * standard libraries vector<...> + * class. + */ + typedef TrilinosScalar value_type; + typedef TrilinosScalar real_type; + typedef size_t size_type; + typedef internal::VectorReference reference; + typedef const internal::VectorReference const_reference; + + /** + * Default constructor. It doesn't do + * anything, derived classes will have + * to initialize the data. + */ + Vector (); + /** + * One of the constructors that + * actually build a vector. This + * one requires prior knowledge + * of the size of the vector and + * a communicator from + * Epetra_CommSerial or Epetra_CommMpi, + * depending on whether we use a + * serial or parallel MPI-based program. + * This command distributes the + * vector linearly among the processes, + * from the beginning to the end, + * so you might want to use some + * more advanced mapping and the + * third constructor with argument + * Epetra_Map. + */ + Vector (unsigned int GlobalSize, Epetra_Comm &Comm); + + /** + * Third constructor. It takes an + * Epetra_Map that already knows how + * to distribute the individual + * components among the MPI processors. + */ + Vector (const Epetra_Map &InputMap); + + /** + * Copy constructor. Sets the dimension + * to that of the given vector and uses + * the map of that vector, and + * copies all elements. + */ + Vector (const Vector &v); + + /** + * Destructor + */ + virtual ~Vector (); + + /** + * Reinit functionality. This function + * destroys the old vector content + * and generates a new one based on + * the input map. + */ + void reinit (const Epetra_Map &input_map); + + /** + * Release all memory and return + * to a state just like after + * having called the default + * constructor. + */ + void clear (); + + /** + * Compress the underlying + * representation of the Trilinos object, + * i.e. flush the buffers of the vector + * object if it has any. This function + * is necessary after writing into a + * vector element-by-element and before + * anything else can be done on it. + */ + void compress (); + + /** + * Set all components of the vector to + * the given number @p s. Simply pass + * this down to the Trilinos Epetra + * object, but we still need to declare + * this function to make the example + * given in the discussion about making + * the constructor explicit work. + * + * Since the semantics of assigning a + * scalar to a vector are not + * immediately clear, this operator + * should really only be used if you + * want to set the entire vector to + * zero. This allows the intuitive + * notation v=0. Assigning + * other values is deprecated and may + * be disallowed in the future. + */ + Vector & operator = (const TrilinosScalar s); + + /** + * Test for equality. This function + * assumes that the present vector and + * the one to compare with have the same + * size already, since comparing vectors + * of different sizes makes not much + * sense anyway. + */ + bool operator == (const Vector &v) const; + + /** + * Test for inequality. This function + * assumes that the present vector and + * the one to compare with have the same + * size already, since comparing vectors + * of different sizes makes not much + * sense anyway. + */ + bool operator != (const Vector &v) const; + + /** + * Return the global dimension of the + * vector. + */ + unsigned int size () const; + + /** + * Return the local dimension of the + * vector, i.e. the number of elements + * stored on the present MPI + * process. For sequential vectors, + * this number is the same as size(), + * but for parallel vectors it may be + * smaller. + * + * To figure out which elements + * exactly are stored locally, + * use local_range(). + */ + unsigned int local_size () const; + + /** + * Return a pair of indices + * indicating which elements of + * this vector are stored + * locally. The first number is + * the index of the first + * element stored, the second + * the index of the one past + * the last one that is stored + * locally. If this is a + * sequential vector, then the + * result will be the pair + * (0,N), otherwise it will be + * a pair (i,i+n), where + * n=local_size(). + */ + std::pair + local_range () const; + + /** + * Return whether @p index is + * in the local range or not, + * see also local_range(). + */ + bool in_local_range (const unsigned int index) const; + + /** + * Provide access to a given element, + * both read and write. + */ + reference + operator () (const unsigned int index); + + /** + * Provide read-only access to an + * element. This is equivalent to + * the el() command. + */ + TrilinosScalar + operator () (const unsigned int index) const; + + /** + * Return the value of the vector entry + * i. Note that this function + * does only work properly when + * we request a data stored on the + * local processor. The function will + * throw an exception in case the + * elements sits on another process. + */ + TrilinosScalar el (const unsigned int index) const; + + /** + * A collective set operation: instead + * of setting individual elements of a + * vector, this function allows to set + * a whole set of elements at once. The + * indices of the elements to be set + * are stated in the first argument, + * the corresponding values in the + * second. + */ + void set (const std::vector &indices, + const std::vector &values); + + /** + * This is a second collective set + * operation. As a difference, this + * function takes a deal.II vector + * of values. + */ + void set (const std::vector &indices, + const dealii::Vector &values); + + /** + * This collective set operation is + * of lower level and can handle + * anything else – the only + * thing you have to provide is + * an address where all the indices + * are stored and the number of + * elements to be set. + */ + void set (const unsigned int n_elements, + const unsigned int *indices, + const TrilinosScalar *values); + + /** + * A collective add operation: This + * function adds a whole set of values + * stored in @p values to the vector + * components specified by @p indices. + */ + void add (const std::vector &indices, + const std::vector &values); + + void add (const unsigned int n_elements, + const unsigned int *indices, + const TrilinosScalar *values); + + /** + * Return the scalar/inner product of two + * vectors. The vectors must have the + * same size. + */ + TrilinosScalar operator * (const Vector &vec) const; + + /** + * Return square of the $l_2$-norm. + */ + real_type norm_sqr () const; + + /** + * Mean value of the elements of + * this vector. + */ + TrilinosScalar mean_value () const; + + /** + * $l_1$-norm of the vector. + * The sum of the absolute values. + */ + real_type l1_norm () const; + + /** + * $l_2$-norm of the vector. The + * square root of the sum of the + * squares of the elements. + */ + real_type l2_norm () const; + + /** + * $l_p$-norm of the vector. The + * pth root of the sum of the + * pth + * powers of the absolute values + * of the elements. + */ + real_type lp_norm (const TrilinosScalar p) const; + + /** + * Maximum absolute value of the + * elements. + */ + real_type linfty_norm () const; + + /** + * Return whether the vector contains + * only elements with value zero. This + * function is mainly for internal + * consistency checks and should + * seldomly be used when not in debug + * mode since it uses quite some time. + */ + bool all_zero () const; + + /** + * Return @p true if the vector has no + * negative entries, i.e. all entries + * are zero or positive. This function + * is used, for example, to check + * whether refinement indicators are + * really all positive (or zero). + */ + bool is_non_negative () const; + + /** + * Multiply the entire vector by a + * fixed factor. + */ + Vector & operator *= (const TrilinosScalar factor); + + /** + * Divide the entire vector by a + * fixed factor. + */ + Vector & operator /= (const TrilinosScalar factor); + + /** + * Add the given vector to the present + * one. + */ + Vector & operator += (const Vector &V); + + /** + * Subtract the given vector from the + * present one. + */ + Vector & operator -= (const Vector &V); + + /** + * Addition of @p s to all + * components. Note that @p s is a + * scalar and not a vector. + */ + void add (const TrilinosScalar s); + + /** + * Simple vector addition, equal to the + * operator =. + */ + void add (const Vector &V); + + /** + * Simple addition of a multiple of a + * vector, i.e. *this = a*V. + */ + void add (const TrilinosScalar a, const Vector &V); + + /** + * Multiple addition of scaled vectors, + * i.e. *this = a*V + b*W. + */ + void add (const TrilinosScalar a, const Vector &V, + const TrilinosScalar b, const Vector &W); + + /** + * Scaling and simple vector addition, + * i.e. + * *this = s*(*this) + V. + */ + void sadd (const TrilinosScalar s, + const Vector &V); + + /** + * Scaling and simple addition, i.e. + * *this = s*(*this) + a*V. + */ + void sadd (const TrilinosScalar s, + const TrilinosScalar a, + const Vector &V); + + /** + * Scaling and multiple addition. + */ + void sadd (const TrilinosScalar s, + const TrilinosScalar a, + const Vector &V, + const TrilinosScalar b, + const Vector &W); + + /** + * Scaling and multiple addition. + * *this = s*(*this) + a*V + b*W + c*X. + */ + void sadd (const TrilinosScalar s, + const TrilinosScalar a, + const Vector &V, + const TrilinosScalar b, + const Vector &W, + const TrilinosScalar c, + const Vector &X); + + /** + * Scale each element of this + * vector by the corresponding + * element in the argument. This + * function is mostly meant to + * simulate multiplication (and + * immediate re-assignment) by a + * diagonal scaling matrix. + */ + void scale (const Vector &scaling_factors); + + /** + * Assignment *this = a*V. + */ + void equ (const TrilinosScalar a, const Vector &V); + + /** + * Assignment *this = a*V + b*W. + */ + void equ (const TrilinosScalar a, const Vector &V, + const TrilinosScalar b, const Vector &W); + + /** + * Compute the elementwise ratio of the + * two given vectors, that is let + * this[i] = a[i]/b[i]. This is + * useful for example if you want to + * compute the cellwise ratio of true to + * estimated error. + * + * This vector is appropriately + * scaled to hold the result. + * + * If any of the b[i] is + * zero, the result is + * undefined. No attempt is made + * to catch such situations. + */ + void ratio (const Vector &a, + const Vector &b); + + /** + * Print to a + * stream. @p precision denotes + * the desired precision with + * which values shall be printed, + * @p scientific whether + * scientific notation shall be + * used. If @p across is + * @p true then the vector is + * printed in a line, while if + * @p false then the elements + * are printed on a separate line + * each. + */ + void print (std::ostream &out, + const unsigned int precision = 3, + const bool scientific = true, + const bool across = true) const; + + /** + * Swap the contents of this + * vector and the other vector + * @p v. One could do this + * operation with a temporary + * variable and copying over the + * data elements, but this + * function is significantly more + * efficient since it only swaps + * the pointers to the data of + * the two vectors and therefore + * does not need to allocate + * temporary storage and move + * data around. Note that the + * vectors need to be of the same + * size and base on the same + * map. + * + * This function is analog to the + * the @p swap function of all C + * standard containers. Also, + * there is a global function + * swap(u,v) that simply calls + * u.swap(v), again in analogy + * to standard functions. + */ + void swap (Vector &v); + + /** + * Exception + */ + DeclException1 (ExcTrilinosError, + int, + << "An error with error number " << arg1 + << " occured while calling a Trilinos function"); + + /** + * Exception + */ + DeclException3 (ExcAccessToNonlocalElement, + int, int, int, + << "You tried to access element " << arg1 + << " of a distributed vector, but only entries " + << arg2 << " through " << arg3 + << " are stored locally and can be accessed."); + + + /** + * An Epetra map used to map vector data + * accross multiple processes. This is + * the communicator and data distribution + * object common to all + * Trilinos objects used by deal.II. + * TODO: we probably only need a pointer + * to the map, since the information + * is provided from outside and there + * is no need to copy the map. + * Especially not when we have many + * vectors based on the same map. + */ + Epetra_Map map; + + /** + * An Epetra distibuted vector type. + * Requires an existing Epetra_Map for + * storing data. + */ + std::auto_ptr vector; + + + /** + * Trilinos doesn't allow to mix additions + * to matrix entries and overwriting + * them (to make synchronisation of + * parallel computations + * simpler). The way we do it is to, + * for each access operation, store + * whether it is an insertion or an + * addition. If the previous one was of + * different type, then we first have + * to flush the Trilinos buffers; + * otherwise, we can simply go on. + * Luckily, Trilinos has an object + * for this which does already all + * the parallel communications in + * such a case, so we simply use their + * model, which stores whether the + * last operation was an addition + * or an insertion. + */ + Epetra_CombineMode last_action; + + /** + * Make the reference class a friend. + */ + friend class internal::VectorReference; + + }; + + + +// ------------------- inline and template functions -------------- + +/** + * Global function @p swap which overloads the default implementation + * of the C standard library which uses a temporary object. The + * function simply exchanges the data of the two vectors. + * + * @ingroup TrilinosWrappers + * @relates TrilinosWrappers::Vector + * @author Wolfgang Bangerth, 2004 + */ + inline + void swap (Vector &u, Vector &v) + { + u.swap (v); + } + +#ifndef DOXYGEN + namespace internal + { + inline + VectorReference::VectorReference (Vector &vector, + const unsigned int index) + : + vector (vector), + index (index) + {} + + + inline + const VectorReference & + VectorReference::operator = (const VectorReference &r) const + { + // as explained in the class + // documentation, this is not the copy + // operator. so simply pass on to the + // "correct" assignment operator + *this = static_cast (r); + + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator = (const TrilinosScalar &value) const + { + vector.set (1, &index, &value); + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator += (const TrilinosScalar &value) const + { + vector.add (1, &index, &value); + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator -= (const TrilinosScalar &value) const + { + TrilinosScalar new_value = -value; + vector.add (1, &index, &new_value); + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator *= (const TrilinosScalar &value) const + { + TrilinosScalar new_value = static_cast(*this) * value; + vector.set (1, &index, &new_value); + return *this; + } + + + + inline + const VectorReference & + VectorReference::operator /= (const TrilinosScalar &value) const + { + TrilinosScalar new_value = static_cast(*this) / value; + vector.set (1, &index, &new_value); + return *this; + } + } + + + + inline + bool + Vector::in_local_range (const unsigned int index) const + { + std::pair range = local_range(); + + return ((index >= range.first) && (index < range.second)); + } + + + + inline + internal::VectorReference + Vector::operator () (const unsigned int index) + { + return internal::VectorReference (*this, index); + } + + + + inline + TrilinosScalar + Vector::operator () (const unsigned int index) const + { + TrilinosScalar value = el(index); + + return value; + } + + +#endif // DOXYGEN +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS + +/*---------------------------- trilinos_vector_base.h ---------------------------*/ + +#endif +/*---------------------------- trilinos_vector_base.h ---------------------------*/ diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc new file mode 100755 index 0000000000..5ce2e21e09 --- /dev/null +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -0,0 +1,687 @@ +//--------------------------------------------------------------------------- +// $Id: trilinos_matrix_base.cc 16014 2008-04-30 16:18:03Z young $ +// Version: $Name$ +// +// Copyright (C) 2004, 2005, 2006, 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + +#include +#include + +#ifdef DEAL_II_USE_TRILINOS + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + namespace MatrixIterators + { + void + SparseMatrix::const_iterator::Accessor:: + visit_present_row () + { + // if we are asked to visit the + // past-the-end line, then simply + // release all our caches and go on + // with life + if (this->a_row == matrix->m()) + { + colnum_cache.reset (); + value_cache.reset (); + + return; + } + + // otherwise first flush Trilinos caches + matrix->compress (); + + // get a representation of the present + // row + int ncols; + int colnums = matrix->n(); + TrilinosScalar *values = new TrilinosScalar(colnums); + + int ierr; + ierr = matrix->matrix->ExtractGlobalRowCopy((int)this->a_row, colnums, + ncols, &(values[0])); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + // copy it into our caches if the line + // isn't empty. if it is, then we've + // done something wrong, since we + // shouldn't have initialized an + // iterator for an empty line (what + // would it point to?) + Assert (ncols != 0, ExcInternalError()); + colnum_cache.reset (new std::vector (colnums, + colnums+ncols)); + value_cache.reset (new std::vector (values, values+ncols)); + } + } + + + // The constructor is actually the + // only point where we have to check + // whether we build a serial or + // a parallel Trilinos matrix. + // In the end, it does not even + // matter how many threads there + // are, but only if we use an + // MPI compiler or a standard + // compiler. So, one thread on + // an MPI compiler will still get + // a parallel interface. + SparseMatrix::SparseMatrix () + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), +#else + map (0,0,Epetra_SerialComm()), +#endif + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, map, 0))), + last_action (Insert) + {} + + + SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, + const unsigned int n_max_entries_per_row) + : + map (InputMap), + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, map, + int(n_max_entries_per_row), false))), + last_action (Insert) + {} + + SparseMatrix::SparseMatrix (const Epetra_Map &InputMap, + const std::vector &n_entries_per_row) + : + map (InputMap), + matrix (std::auto_ptr + (new Epetra_FECrsMatrix(Copy, map, + (int*)const_cast(&(n_entries_per_row[0])), + true))), + last_action (Insert) + {} + + + SparseMatrix::~SparseMatrix () + { + } + + + + void + SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern, + const unsigned int n_max_entries_per_row) + { + + unsigned int n_rows = sparsity_pattern.n_rows(); + + Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (matrix->NumGlobalRows(), + sparsity_pattern.n_rows())); + + std::vector values(n_max_entries_per_row, 0.); + std::vector row_indices(n_max_entries_per_row); + + for (unsigned int row=0; rowInsertGlobalValues(row, row_length, + &values[0], &row_indices[0]); + } + + // In the end, the matrix is to + // be compressed in order to be + // really ready. However, that is + // a collective operation, so it + // has to be called on all processes + // by the user, whereas this function + // should only be used on one processor + // since our sparsity pattern data + // types are all serial. + } + + + + void + SparseMatrix::reinit (const CompressedSparsityPattern &sparsity_pattern) + { + unsigned int n_rows = sparsity_pattern.n_rows(); + + Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), + ExcDimensionMismatch (matrix->NumGlobalRows(), + sparsity_pattern.n_rows())); + + std::vector n_entries_per_row(n_rows); + + for (unsigned int row=0; row n_entries_per_row(n_rows); + + for (unsigned int row=0; row + (new Epetra_FECrsMatrix(Copy, map, &n_entries_per_row[0], true)); + + const unsigned int n_max_entries_per_row = *std::max_element ( + &n_entries_per_row[0], &n_entries_per_row[n_rows-1]); + + reinit (sparsity_pattern, n_max_entries_per_row); + } + + + + void + SparseMatrix::clear () + { + // When we clear the matrix, + // reset the pointer and + // generate an empty matrix. + matrix.reset(); +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), +#else + map = Epetra_Map (0,0,Epetra_SerialComm()), +#endif + + matrix = std::auto_ptr + (new Epetra_FECrsMatrix(Copy, map, 0)); + } + + + + void + SparseMatrix::compress () + { + // flush buffers + int ierr; + ierr = matrix->GlobalAssemble (true); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + ierr = matrix->OptimizeStorage (); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + SparseMatrix & + SparseMatrix::operator = (const double d) + { + Assert (d==0, ExcScalarAssignmentOnlyForZeroValue()); + + compress (); + + const int ierr = matrix->PutScalar(d); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + + void + SparseMatrix::set (const unsigned int i, + const unsigned int j, + const TrilinosScalar value) + { + + Assert (numbers::is_finite(value), + ExcMessage("The given value is not finite but either " + "infinite or Not A Number (NaN)")); + + if (last_action == Add) + { + int ierr; + ierr = matrix->GlobalAssemble(false); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + last_action = Insert; + } + + int trilinos_i = i; + int trilinos_j = j; + + const int ierr = matrix->ReplaceGlobalValues (trilinos_i, 1, + const_cast(&value), + &trilinos_j); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + SparseMatrix::add (const unsigned int i, + const unsigned int j, + const TrilinosScalar value) + { + + Assert (numbers::is_finite(value), + ExcMessage("The given value is not finite but either " + "infinite or Not A Number (NaN)")); + + if (last_action == Insert) + { + int ierr; + ierr = matrix->GlobalAssemble(false); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + last_action = Add; + } + + // we have to do above actions in any + // case to be consistent with the MPI + // communication model (see the + // comments in the documentation of + // TrilinosWrappers::MPI::Vector), but we + // can save some work if the addend is + // zero + if (value == 0) + return; + + int trilinos_i = i; + int trilinos_j = j; + + const int ierr = matrix->SumIntoGlobalValues (trilinos_i, 1, + const_cast(&value), + &trilinos_j); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + TrilinosScalar + SparseMatrix::el (const unsigned int i, + const unsigned int j) const + { + // Extract local indices in + // the matrix. + int trilinos_i = matrix->LRID(i), trilinos_j = matrix->LRID(j); + TrilinosScalar value = 0.; + + // If the data is not on the + // present processor, we can't + // continue. + if ((trilinos_i == -1 ) || (trilinos_j == -1)) + { + Assert (false, ExcAccessToNonlocalElement(i, j, local_range().first, + local_range().second)); + } + else + { + // Check whether the matrix + // already is transformed to + // local indices. + if (!matrix->Filled()) + matrix->FillComplete(true); + + // Prepare pointers for extraction + // of a view of the row. + int nnz_present = matrix->NumMyEntries(trilinos_i); + int nnz_extracted; + int *col_indices; + TrilinosScalar *values; + + // Generate the view and make + // sure that we have not generated + // an error. + int ierr = matrix->ExtractMyRowView(trilinos_i, nnz_extracted, + values, col_indices); + Assert (ierr==0, ExcTrilinosError(ierr)); + + Assert (nnz_present == nnz_extracted, + ExcDimensionMismatch(nnz_present, nnz_extracted)); + + // Search the index where we + // look for the value, and then + // finally get it. + int* index = std::find(&col_indices[0],&col_indices[0] + nnz_present, + trilinos_j); + + int position; + if (!index) + value = 0; + else + { + position = (int)(index - &(col_indices[0])); + value = values[position]; + } + } + + return value; + } + + + + TrilinosScalar + SparseMatrix::diag_element (const unsigned int i) const + { + Assert (m() == n(), ExcNotQuadratic()); + + // this doesn't seem to work any + // different than any other element + return el(i,i); + } + + + + unsigned int + SparseMatrix::m () const + { + int n_rows = matrix->NumGlobalRows(); + + return n_rows; + } + + + + unsigned int + SparseMatrix::n () const + { + int n_cols = matrix -> NumGlobalCols(); + return n_cols; + } + + + + unsigned int + SparseMatrix::local_size () const + { + int n_rows = matrix -> NumMyRows(); + + return n_rows; + } + + + + std::pair + SparseMatrix::local_range () const + { + int begin, end; + begin = matrix->RowMap().MinMyGID(); + end = matrix->RowMap().MaxMyGID(); + + return std::make_pair (begin, end); + } + + + + unsigned int + SparseMatrix::n_nonzero_elements () const + { + int nnz = matrix->NumGlobalNonzeros(); + + return static_cast(nnz); + } + + + + unsigned int + SparseMatrix::row_length (const unsigned int row) const + { + Assert (row < m(), ExcInternalError()); + + // get a representation of the present + // row + int ncols = -1; + int local_row = matrix->RowMap().LID(row); + + // on the processor who owns this + // row, we'll have a non-negative + // value. + if (local_row >= 0) + { + int ierr = matrix->NumMyRowEntries (local_row, ncols); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + return ncols; + } + + + TrilinosScalar + SparseMatrix::l1_norm () const + { + if (!matrix->Filled()) + matrix->FillComplete(); + + TrilinosScalar result = matrix->NormOne(); + + return result; + } + + + + TrilinosScalar + SparseMatrix::linfty_norm () const + { + if (!matrix->Filled()) + matrix->FillComplete(); + + TrilinosScalar result = matrix->NormInf(); + + return result; + } + + + + TrilinosScalar + SparseMatrix::frobenius_norm () const + { + if (!matrix->Filled()) + matrix->FillComplete(); + + TrilinosScalar result = matrix->NormFrobenius(); + + return result; + } + + + + SparseMatrix & + SparseMatrix::operator *= (const TrilinosScalar a) + { + const int ierr = matrix->Scale (a); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + + SparseMatrix & + SparseMatrix::operator /= (const TrilinosScalar a) + { + const TrilinosScalar factor = 1./a; + + const int ierr = matrix->Scale (factor); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + void + SparseMatrix::vmult (Vector &dst, + const Vector &src) const + { + Assert (&src != &dst, ExcSourceEqualsDestination()); + + if (!matrix->Filled()) + matrix->FillComplete(); + + const int ierr = matrix->Multiply (false, *(src.vector), *(dst.vector)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + SparseMatrix::Tvmult (Vector &dst, + const Vector &src) const + { + Assert (&src != &dst, ExcSourceEqualsDestination()); + + if (!matrix->Filled()) + matrix->FillComplete(); + + const int ierr = matrix->Multiply (true, *(src.vector), *(dst.vector)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + SparseMatrix::vmult_add (Vector &dst, + const Vector &src) const + { + Assert (&src != &dst, ExcSourceEqualsDestination()); + + Vector tmp = dst; + vmult (dst, src); + dst += tmp; + } + + + + void + SparseMatrix::Tvmult_add (Vector &dst, + const Vector &src) const + { + Assert (&src != &dst, ExcSourceEqualsDestination()); + + Vector tmp = dst; + vmult (dst, src); + dst += tmp; + } + + + + TrilinosScalar + SparseMatrix::matrix_norm_square (const Vector &v) const + { + Vector tmp(v.map); + vmult (tmp, v); + return tmp*v; + } + + + + TrilinosScalar + SparseMatrix::matrix_scalar_product (const Vector &u, + const Vector &v) const + { + Vector tmp(v.map); + vmult (tmp, v); + return u*tmp; + } + + + + TrilinosScalar + SparseMatrix::residual (Vector &dst, + const Vector &x, + const Vector &b) const + { + vmult (dst, x); + dst -= b; + dst *= -1; + + return dst.l2_norm(); + } + + + + // TODO: Currently this only flips + // a flag that tells Trilinos that + // any application should be done with + // the transpose. However, the matrix + // structure is not reset. + void + SparseMatrix::transpose () + { + int ierr; + + if (!matrix->UseTranspose()) + { + ierr = matrix->SetUseTranspose (true); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + else + { + ierr = matrix->SetUseTranspose (false); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + } + + + + bool + SparseMatrix::is_symmetric (const double tolerance) + { + //bool truth; + if (tolerance == 0) + Assert (false, ExcNotImplemented()); + + return false; + } + + + + bool + SparseMatrix::is_hermitian () + { + //bool truth; + + Assert (false, ExcNotImplemented()); + return false; + } + + void + SparseMatrix::write_ascii () + { + Assert (false, ExcNotImplemented()); + } + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc new file mode 100755 index 0000000000..18f8afc3a8 --- /dev/null +++ b/deal.II/lac/source/trilinos_vector.cc @@ -0,0 +1,786 @@ +//--------------------------------------------------------------------------- +// $Id: trilinos_vector_base.cc 15631 2008-01-29 01:00:00Z marquez $ +// Version: $Name$ +// +// Copyright (C) 2008 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//--------------------------------------------------------------------------- + + +#include + +#include + +#ifdef DEAL_II_USE_TRILINOS + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + namespace internal + { + VectorReference::operator TrilinosScalar () const // i believe useless with + // trilinos + { + Assert (index < vector.size(), + ExcIndexRange (index, 0, vector.size())); + + // Trilinos allows for vectors to be + // referenced by the [] or () operators + // but only () checks index bounds + // Also, can only get local values + + AssertThrow ((static_cast(index) >= vector.map.MinMyGID()) && + (static_cast(index) <= vector.map.MaxMyGID()), + ExcAccessToNonlocalElement (index, vector.map.MinMyGID(), + vector.map.MaxMyGID()-1)); + return *(*(vector.vector))[index]; + } + } + + Vector::Vector () + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), +#else + map (0,0,Epetra_SerialComm()), +#endif + vector(std::auto_ptr + (new Epetra_FEVector(map))), + last_action (Insert) + {} + + Vector::Vector (unsigned int GlobalSize, Epetra_Comm &Comm) + : + map (GlobalSize, 0, Comm), + vector (std::auto_ptr + (new Epetra_FEVector(map))), + last_action (Insert) + {} + + + Vector::Vector (const Epetra_Map &InputMap) + : + map (InputMap), + vector (std::auto_ptr + (new Epetra_FEVector(map))), + last_action (Insert) + {} + + + Vector::Vector (const Vector &v) + : + map (v.map), + vector(std::auto_ptr + (new Epetra_FEVector(*(v.vector)))), + last_action (Insert) + {} + + + + Vector::~Vector () + {} + + + + void + Vector::reinit (const Epetra_Map &input_map) + { + vector.reset(); + map = input_map; + + vector = std::auto_ptr (new Epetra_FEVector(input_map)); + last_action = Insert; + } + + + + void + Vector::clear () + { + // When we clear the matrix, + // reset the pointer and + // generate an empty matrix. + vector.reset(); +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + map = Epetra_Map (0,0,Epetra_MpiComm(MPI_COMM_WORLD)), +#else + map = Epetra_Map (0,0,Epetra_SerialComm()), +#endif + + vector = std::auto_ptr (new Epetra_FEVector(map)); + } + + + + void + Vector::compress () + { + // Now pass over the information + // about what we did last to Trilinos. + const int ierr = vector->GlobalAssemble(last_action); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + Vector & + Vector::operator = (const TrilinosScalar s) + { + + Assert (numbers::is_finite(s), + ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); + + const int ierr = vector->PutScalar(s); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + + bool + Vector::operator == (const Vector &v) const + { + Assert (size() == v.size(), + ExcDimensionMismatch(size(), v.size())); + unsigned int i; + for (i=0; iMap().NumGlobalElements(); + } + + + + unsigned int + Vector::local_size () const + { + return (unsigned int) vector->Map().NumMyElements(); + } + + + + std::pair + Vector::local_range () const + { + int begin, end; + begin = vector->Map().MinMyGID(); + end = vector->Map().MaxMyGID(); + return std::make_pair (begin, end); + } + + + + TrilinosScalar + Vector::el (const unsigned int index) const + { + // Extract local indices in + // the vector. + int trilinos_i = map.LID(index); + TrilinosScalar value = 0.; + if (trilinos_i == -1 ) + { + Assert (false, ExcAccessToNonlocalElement(index, local_range().first, + local_range().second)); + } + else + value = (*vector)[0][trilinos_i]; + + return value; + } + + + void + Vector::set (const std::vector &indices, + const std::vector &values) + { + Assert (indices.size() == values.size(), + ExcMessage ("Function called with arguments of different sizes")); + + set (indices.size(), &indices[0], &values[0]); + } + + + void + Vector::set (const unsigned int n_elements, + const unsigned int *indices, + const TrilinosScalar *values) + { + if (last_action == Add) + { + vector->GlobalAssemble(Add); + last_action = Insert; + } + + const int ierr= vector->ReplaceGlobalValues (n_elements, + (int*)(const_cast(indices)), + const_cast(values)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::add (const std::vector &indices, + const std::vector &values) + { + Assert (indices.size() == values.size(), + ExcMessage ("Function called with arguments of different sizes")); + + add (indices.size(), &indices[0], &values[0]); + } + + + + void + Vector::add (const unsigned int n_elements, + const unsigned int *indices, + const TrilinosScalar *values) + { + if (last_action == Insert) + { + vector->GlobalAssemble(Insert); + last_action = Add; + } + + const int ierr= vector->SumIntoGlobalValues (n_elements, + (int*)(const_cast(indices)), + const_cast(values)); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + TrilinosScalar + Vector::operator * (const Vector &vec) const + { + Assert (size() == vec.size(), + ExcDimensionMismatch(size(), vec.size())); + + TrilinosScalar result; + + const int ierr = vector->Dot(*(vec.vector), &result); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return result; + } + + + + Vector::real_type + Vector::norm_sqr () const + { + const TrilinosScalar d = l2_norm(); + return d*d; + } + + + + TrilinosScalar + Vector::mean_value () const + { + + TrilinosScalar mean; + int ierr; + + ierr = vector->MeanValue (&mean); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return mean; + } + + + + Vector::real_type + Vector::l1_norm () const + { + TrilinosScalar d; + + const int ierr = vector->Norm1 (&d); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return d; + } + + + + Vector::real_type + Vector::l2_norm () const + { + TrilinosScalar d; + + const int ierr = vector->Norm2 (&d); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return d; + } + + + + Vector::real_type + Vector::lp_norm (const TrilinosScalar p) const + { + // get a representation of the vector and + // loop over all the elements + TrilinosScalar *start_ptr; + int leading_dimension; + int ierr = vector->ExtractView (&start_ptr, &leading_dimension); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + TrilinosScalar norm = 0; + TrilinosScalar sum=0; + + const TrilinosScalar * ptr = start_ptr; + // add up elements + while (ptr != start_ptr+size()) + sum += std::pow(std::fabs(*ptr++), p); + + norm = std::pow(sum, static_cast(1./p)); + + return norm; + } + + + + Vector::real_type + Vector::linfty_norm () const + { + TrilinosScalar d; + + const int ierr = vector->NormInf (&d); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return d; + } + + + + bool + Vector::all_zero () const + { + // get a representation of the vector and + // loop over all the elements + TrilinosScalar *start_ptr; + int leading_dimension; + int ierr = vector->ExtractView (&start_ptr, &leading_dimension); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + const TrilinosScalar *ptr = start_ptr, + *eptr = start_ptr + size(); + bool flag = true; + while (ptr != eptr) + { + if (*ptr != 0) + { + flag = false; + break; + } + ++ptr; + } + + return flag; + } + + + + bool + Vector::is_non_negative () const + { + // get a representation of the vector and + // loop over all the elements + TrilinosScalar *start_ptr; + int leading_dimension; + int ierr = vector->ExtractView (&start_ptr, &leading_dimension); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + const TrilinosScalar *ptr = start_ptr, + *eptr = start_ptr + size(); + bool flag = true; + while (ptr != eptr) + { + if (*ptr < 0.0) + { + flag = false; + break; + } + ++ptr; + } + + return flag; + } + + + + Vector & + Vector::operator *= (const TrilinosScalar a) + { + + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Scale(a); + Assert (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + + Vector & + Vector::operator /= (const TrilinosScalar a) + { + + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const TrilinosScalar factor = 1./a; + + Assert (numbers::is_finite(factor), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Scale(factor); + Assert (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + + Vector & + Vector::operator += (const Vector &v) + { + const int ierr = vector->Update (1.0, *(v.vector), 1.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + return *this; + } + + + + Vector & + Vector::operator -= (const Vector &v) + { + const int ierr = vector->Update (-1.0, *(v.vector), 1.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + return *this; + } + + + + void + Vector::add (const TrilinosScalar s) + { + + Assert (numbers::is_finite(s), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + std::vector list (size(), s); + + int* index = new int[size()]; + for (unsigned int i=0; iSumIntoGlobalValues(size(), index, &list[0]); + + delete[] index; + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::add (const Vector &v) + { + *this += v; + } + + + + void + Vector::add (const TrilinosScalar a, + const Vector &v) + { + + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Update(a, *(v.vector), 1); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::add (const TrilinosScalar a, + const Vector &v, + const TrilinosScalar b, + const Vector &w) + { + + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(b), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.0); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::sadd (const TrilinosScalar s, + const Vector &v) + { + + Assert (numbers::is_finite(s), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Update(1.0, *(v.vector), s); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::sadd (const TrilinosScalar s, + const TrilinosScalar a, + const Vector &v) + { + + Assert (numbers::is_finite(s), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Update(a, *(v.vector), s); + + AssertThrow(ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::sadd (const TrilinosScalar s, + const TrilinosScalar a, + const Vector &v, + const TrilinosScalar b, + const Vector &w) + { + + Assert (numbers::is_finite(s), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(b), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::sadd (const TrilinosScalar s, + const TrilinosScalar a, + const Vector &v, + const TrilinosScalar b, + const Vector &w, + const TrilinosScalar c, + const Vector &x) + { + + Assert (numbers::is_finite(s), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(b), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(c), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + // Update member can only input two other vectors so + // do it in two steps + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), s); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + const int jerr = vector->Update(c, *(x.vector), 1.0); + AssertThrow (jerr == 0, ExcTrilinosError(jerr)); + } + + + + void + Vector::scale (const Vector &factors) + { + const int ierr = vector->Multiply (1.0, *(factors.vector), *vector, 0.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::equ (const TrilinosScalar a, + const Vector &v) + { + + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); + + const int ierr = vector->Update(a, *(v.vector), 0.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::equ (const TrilinosScalar a, + const Vector &v, + const TrilinosScalar b, + const Vector &w) + { + + Assert (numbers::is_finite(a), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + Assert (numbers::is_finite(b), + ExcMessage("The given value is not finite but " + "either infinite or Not A Number (NaN)")); + + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); + + Assert (size() == w.size(), + ExcDimensionMismatch (size(), w.size())); + + const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 0.0); + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + + } + + + + void + Vector::ratio (const Vector &v, + const Vector &w) + { + Assert (size() == v.size(), + ExcDimensionMismatch (size(), v.size())); + + Assert (size() == w.size(), + ExcDimensionMismatch (size(), w.size())); + + const int ierr = vector->ReciprocalMultiply(1.0, *(w.vector), + *(v.vector), 0.0); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + } + + + + void + Vector::print (std::ostream &out, + const unsigned int precision, + const bool scientific, + const bool across) const + { + AssertThrow (out, ExcIO()); + + // get a representation of the vector and + // loop over all the elements + // TODO: up to now only local data + // printed out! + TrilinosScalar *val; + int leading_dimension; + int ierr = vector->ExtractView (&val, &leading_dimension); + + AssertThrow (ierr == 0, ExcTrilinosError(ierr)); + out.precision (precision); + if (scientific) + out.setf (std::ios::scientific, std::ios::floatfield); + else + out.setf (std::ios::fixed, std::ios::floatfield); + + if (across) + for (unsigned int i=0; i(val[i]) << ' '; + else + for (unsigned int i=0; i(val[i]) << std::endl; + out << std::endl; + + // restore the representation of the + // vector + AssertThrow (out, ExcIO()); + } + + + + void + Vector::swap (Vector &v) + { + // Just swap the pointers to the + // two Epetra vectors that hold all + // the data. + std::auto_ptr tmp; + tmp = v.vector; + v.vector = vector; + vector = tmp; + } + +} + +DEAL_II_NAMESPACE_CLOSE + +#endif // DEAL_II_USE_TRILINOS