--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/compressed_sparsity_pattern.h>
+#include <lac/exceptions.h>
+#include <lac/trilinos_vector.h>
+
+#include <boost/shared_ptr.hpp>
+#include <vector>
+#include <cmath>
+#include <memory>
+
+#ifdef DEAL_II_USE_TRILINOS
+
+# define TrilinosScalar double
+# include <Epetra_FECrsMatrix.h>
+# include <Epetra_Map.h>
+# include <Epetra_CrsGraph.h>
+# ifdef DEAL_II_COMPILER_SUPPORTS_MPI
+# include <Epetra_MpiComm.h>
+# 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<const std::vector<unsigned int> > colnum_cache;
+
+ /**
+ * Similar cache for the values
+ * of this row.
+ */
+ boost::shared_ptr<const std::vector<TrilinosScalar> > 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 <tt>==</tt>.
+ */
+ 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<unsigned int> &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 <tt>matrix=0</tt>, 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>i,j</i>)
+ * 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 <tt>value</tt> 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>i,j</i>).
+ *
+ * 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 <tt>value</tt> 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 <tt>row</tt> 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<unsigned int> &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>i,j</i>). 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 <tt>el()</tt> function.
+ */
+ TrilinosScalar operator () (const unsigned int i,
+ const unsigned int j) const;
+
+ /**
+ * Return the value of the matrix entry
+ * (<i>i,j</i>). 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 <i>i</i>th
+ * 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
+ * <tt>n=local_size()</tt>.
+ */
+ std::pair<unsigned int, unsigned int>
+ 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 <i>dst = M*src</i> with
+ * <i>M</i> 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
+ * <i>dst = M<sup>T</sup>*src</i> with
+ * <i>M</i> 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
+ * <i>M*src</i> on <i>dst</i>
+ * with <i>M</i> 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
+ * <i>M<sup>T</sup>*src</i> to
+ * <i>dst</i> with <i>M</i> 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 <i>Mx=b</i>, where
+ * the residual is defined to be
+ * <i>r=b-Mx</i>. Write the
+ * residual into
+ * @p dst. The
+ * <i>l<sub>2</sub></i> norm of
+ * the residual vector is
+ * returned.
+ *
+ * Source <i>x</i> and destination
+ * <i>dst</i> 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
+ * <tt>end(r)</tt>. Note also that the
+ * iterator may not be dereferencable in
+ * that case.
+ */
+ const_iterator begin (const unsigned int r) const;
+
+ /**
+ * Final iterator of row <tt>r</tt>. 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<Epetra_FECrsMatrix> 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<SparseMatrix*>(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<m(); ++i)
+ if (row_length(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<unsigned int>(begin)) &&
+ (index < static_cast<unsigned int>(end)));
+ }
+
+#endif // DOXYGEN
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+
+#endif // DEAL_II_USE_TRILINOS
+
+
+/*---------------------------- trilinos_matrix_base.h ---------------------------*/
+
+#endif
+/*---------------------------- trilinos_matrix_base.h ---------------------------*/
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <base/config.h>
+#include <base/subscriptor.h>
+#include <lac/exceptions.h>
+#include <lac/vector.h>
+
+#include <vector>
+#include <utility>
+#include <memory>
+
+#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 <typename number> 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
+ * <tt>v(i)=w(i)</tt>. 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 <tt>s</tt>.
+ */
+ const VectorReference & operator = (const TrilinosScalar &s) const;
+
+ /**
+ * Add <tt>s</tt> to the referenced
+ * element of the vector->
+ */
+ const VectorReference & operator += (const TrilinosScalar &s) const;
+
+ /**
+ * Subtract <tt>s</tt> from the
+ * referenced element of the vector->
+ */
+ const VectorReference & operator -= (const TrilinosScalar &s) const;
+
+ /**
+ * Multiply the referenced element of
+ * the vector by <tt>s</tt>.
+ */
+ const VectorReference & operator *= (const TrilinosScalar &s) const;
+
+ /**
+ * Divide the referenced element of
+ * the vector by <tt>s</tt>.
+ */
+ 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.
+ *
+ * <h3>Parallel communication model</h3>
+ *
+ * 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.
+ *
+ * <h3>Accessing individual elements of a vector</h3>
+ *
+ * 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 <tt>d=vec(i)</tt>. 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, <tt>vec(i)=d</tt> or <tt>vec(i)+=d</tt>, 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; i<vector->size(); ++i)
+ * vector(i) = i;
+ *
+ * // do some additions to vector elements, but
+ * // only for some elements
+ * for (unsigned int i=0; i<vector->size(); ++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 <tt>C</tt>
+ * standard libraries <tt>vector<...></tt>
+ * 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 <tt>v=0</tt>. 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
+ * <tt>n=local_size()</tt>.
+ */
+ std::pair<unsigned int, unsigned int>
+ 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 <code>el()</code> command.
+ */
+ TrilinosScalar
+ operator () (const unsigned int index) const;
+
+ /**
+ * Return the value of the vector entry
+ * <i>i</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<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &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<unsigned int> &indices,
+ const dealii::Vector<TrilinosScalar> &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<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &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
+ * <i>p</i>th root of the sum of the
+ * <i>p</i>th
+ * 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
+ * <tt>operator =</tt>.
+ */
+ void add (const Vector &V);
+
+ /**
+ * Simple addition of a multiple of a
+ * vector, i.e. <tt>*this = a*V</tt>.
+ */
+ void add (const TrilinosScalar a, const Vector &V);
+
+ /**
+ * Multiple addition of scaled vectors,
+ * i.e. <tt>*this = a*V + b*W</tt>.
+ */
+ void add (const TrilinosScalar a, const Vector &V,
+ const TrilinosScalar b, const Vector &W);
+
+ /**
+ * Scaling and simple vector addition,
+ * i.e.
+ * <tt>*this = s*(*this) + V</tt>.
+ */
+ void sadd (const TrilinosScalar s,
+ const Vector &V);
+
+ /**
+ * Scaling and simple addition, i.e.
+ * <tt>*this = s*(*this) + a*V</tt>.
+ */
+ 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.
+ * <tt>*this = s*(*this) + a*V + b*W + c*X</tt>.
+ */
+ 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 <tt>*this = a*V</tt>.
+ */
+ void equ (const TrilinosScalar a, const Vector &V);
+
+ /**
+ * Assignment <tt>*this = a*V + b*W</tt>.
+ */
+ 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
+ * <tt>this[i] = a[i]/b[i]</tt>. 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 <tt>b[i]</tt> 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
+ * <tt>swap(u,v)</tt> that simply calls
+ * <tt>u.swap(v)</tt>, 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<Epetra_FEVector> 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<TrilinosScalar> (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<TrilinosScalar>(*this) * value;
+ vector.set (1, &index, &new_value);
+ return *this;
+ }
+
+
+
+ inline
+ const VectorReference &
+ VectorReference::operator /= (const TrilinosScalar &value) const
+ {
+ TrilinosScalar new_value = static_cast<TrilinosScalar>(*this) / value;
+ vector.set (1, &index, &new_value);
+ return *this;
+ }
+ }
+
+
+
+ inline
+ bool
+ Vector::in_local_range (const unsigned int index) const
+ {
+ std::pair<unsigned int, unsigned int> 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 ---------------------------*/
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <lac/trilinos_vector.h>
+#include <lac/trilinos_sparse_matrix.h>
+
+#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<unsigned int> (colnums,
+ colnums+ncols));
+ value_cache.reset (new std::vector<TrilinosScalar> (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<Epetra_FECrsMatrix>
+ (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<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, map,
+ int(n_max_entries_per_row), false))),
+ last_action (Insert)
+ {}
+
+ SparseMatrix::SparseMatrix (const Epetra_Map &InputMap,
+ const std::vector<unsigned int> &n_entries_per_row)
+ :
+ map (InputMap),
+ matrix (std::auto_ptr<Epetra_FECrsMatrix>
+ (new Epetra_FECrsMatrix(Copy, map,
+ (int*)const_cast<unsigned int*>(&(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<double> values(n_max_entries_per_row, 0.);
+ std::vector<int> row_indices(n_max_entries_per_row);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ {
+ const int row_length = sparsity_pattern.row_length(row);
+ row_indices.resize (row_length, 0);
+ values.resize (row_length, 0.);
+
+ for (int col=0; col< row_length; ++col)
+ row_indices[col] = sparsity_pattern.column_number (row, col);
+
+ matrix->InsertGlobalValues(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<int> n_entries_per_row(n_rows);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+
+ 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::reinit (const Epetra_Map &input_map,
+ const CompressedSparsityPattern &sparsity_pattern)
+ {
+
+ unsigned int n_rows = sparsity_pattern.n_rows();
+ matrix.reset();
+ map = input_map;
+ Assert (input_map.NumGlobalElements() == (int)sparsity_pattern.n_rows(),
+ ExcDimensionMismatch (input_map.NumGlobalElements(),
+ sparsity_pattern.n_rows()));
+
+ std::vector<int> n_entries_per_row(n_rows);
+
+ for (unsigned int row=0; row<n_rows; ++row)
+ n_entries_per_row[(int)row] = sparsity_pattern.row_length(row);
+
+ matrix = std::auto_ptr<Epetra_FECrsMatrix>
+ (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<Epetra_FECrsMatrix>
+ (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<double*>(&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<double*>(&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<unsigned int, unsigned int>
+ 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<unsigned int>(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
--- /dev/null
+//---------------------------------------------------------------------------
+// $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 <lac/trilinos_vector.h>
+
+#include <cmath>
+
+#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<signed int>(index) >= vector.map.MinMyGID()) &&
+ (static_cast<signed int>(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<Epetra_FEVector>
+ (new Epetra_FEVector(map))),
+ last_action (Insert)
+ {}
+
+ Vector::Vector (unsigned int GlobalSize, Epetra_Comm &Comm)
+ :
+ map (GlobalSize, 0, Comm),
+ vector (std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(map))),
+ last_action (Insert)
+ {}
+
+
+ Vector::Vector (const Epetra_Map &InputMap)
+ :
+ map (InputMap),
+ vector (std::auto_ptr<Epetra_FEVector>
+ (new Epetra_FEVector(map))),
+ last_action (Insert)
+ {}
+
+
+ Vector::Vector (const Vector &v)
+ :
+ map (v.map),
+ vector(std::auto_ptr<Epetra_FEVector>
+ (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<Epetra_FEVector> (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<Epetra_FEVector> (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; i<size(); i++)
+ if ((*(v.vector))[i]!=(*vector)[i]) return false;
+ return true;
+ }
+
+
+
+ bool
+ Vector::operator != (const Vector &v) const
+ {
+ Assert (size() == v.size(),
+ ExcDimensionMismatch(size(), v.size()));
+
+ return (!(*this==v));
+ }
+
+
+
+ unsigned int
+ Vector::size () const
+ {
+ return (unsigned int) vector->Map().NumGlobalElements();
+ }
+
+
+
+ unsigned int
+ Vector::local_size () const
+ {
+ return (unsigned int) vector->Map().NumMyElements();
+ }
+
+
+
+ std::pair<unsigned int, unsigned int>
+ 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<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &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<unsigned int*>(indices)),
+ const_cast<TrilinosScalar*>(values));
+ AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+ }
+
+
+
+ void
+ Vector::add (const std::vector<unsigned int> &indices,
+ const std::vector<TrilinosScalar> &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<unsigned int*>(indices)),
+ const_cast<TrilinosScalar*>(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<TrilinosScalar>(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<TrilinosScalar> list (size(), s);
+
+ int* index = new int[size()];
+ for (unsigned int i=0; i<size(); i++)
+ index[i]=i;
+
+ const int ierr = vector->SumIntoGlobalValues(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<size(); ++i)
+ out << static_cast<double>(val[i]) << ' ';
+ else
+ for (unsigned int i=0; i<size(); ++i)
+ out << static_cast<double>(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<Epetra_FEVector> tmp;
+ tmp = v.vector;
+ v.vector = vector;
+ vector = tmp;
+ }
+
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_USE_TRILINOS