From 233b915b6307ffe7117a7b911a02859f3d75ceb4 Mon Sep 17 00:00:00 2001 From: kronbichler Date: Fri, 12 Sep 2008 09:40:38 +0000 Subject: [PATCH] Introduced a class LocalizedVector that does take a parallel vector and generates copies of all the content on all processors. Need to fix the = operator to work in parallel as well. git-svn-id: https://svn.dealii.org/trunk@16814 0785d39b-7218-0410-832d-ea1e28bc413d --- .../lac/include/lac/trilinos_sparse_matrix.h | 241 ++++++++++-------- deal.II/lac/include/lac/trilinos_vector.h | 96 ++++++- deal.II/lac/include/lac/vector.templates.h | 11 +- deal.II/lac/source/trilinos_sparse_matrix.cc | 40 ++- deal.II/lac/source/trilinos_vector.cc | 73 ++++++ 5 files changed, 334 insertions(+), 127 deletions(-) diff --git a/deal.II/lac/include/lac/trilinos_sparse_matrix.h b/deal.II/lac/include/lac/trilinos_sparse_matrix.h index d5b46b3ef5..18359e99c0 100755 --- a/deal.II/lac/include/lac/trilinos_sparse_matrix.h +++ b/deal.II/lac/include/lac/trilinos_sparse_matrix.h @@ -477,6 +477,12 @@ namespace TrilinosWrappers * holds the sparsity_pattern * structure because each * processor sets its rows. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void reinit (const SparsityPattern &sparsity_pattern); @@ -508,7 +514,14 @@ namespace TrilinosWrappers * necessary that each processor * holds the sparsity_pattern * structure because each - * processor sets its rows. + * processor sets its + * rows. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void reinit (const Epetra_Map &input_map, const SparsityPattern &sparsity_pattern); @@ -522,6 +535,12 @@ namespace TrilinosWrappers * user-supplied Epetra maps. * To be used for rectangular * matrices. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void reinit (const Epetra_Map &input_row_map, const Epetra_Map &input_col_map, @@ -532,6 +551,12 @@ namespace TrilinosWrappers * content in * sparse_matrix to * the calling matrix. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void reinit (const SparseMatrix &sparse_matrix); @@ -546,6 +571,12 @@ namespace TrilinosWrappers * threshold (so zeros in the * deal.II matrix can be * filtered away). + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void reinit (const Epetra_Map &input_map, const ::dealii::SparseMatrix &dealii_sparse_matrix, @@ -560,6 +591,12 @@ namespace TrilinosWrappers * the rows and the columns of * the matrix. Chosen for * rectangular matrices. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void reinit (const Epetra_Map &input_row_map, const Epetra_Map &input_col_map, @@ -571,6 +608,12 @@ namespace TrilinosWrappers * return to a state just like * after having called the * default constructor. + * + * This is a + * collective operation that + * needs to be called on all + * processors in order to avoid a + * dead lock. */ void clear (); @@ -620,12 +663,14 @@ namespace TrilinosWrappers * 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 + * Just as the respective call in + * deal.II SparseMatrix + * class (but in contrast to the + * situation for PETSc based + * matrices), this function + * throws an exception if an + * entry does not exist in the + * sparsity pattern. Moreover, if * value is not a finite * number an exception is thrown. */ @@ -637,16 +682,16 @@ namespace TrilinosWrappers * 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. + * Just as the respective call in + * deal.II SparseMatrix + * class (but in contrast to the + * situation for PETSc based + * matrices), this function + * throws an exception if an + * entry does not exist in the + * sparsity pattern. Moreover, if + * value is not a finite + * number an exception is thrown. */ void add (const unsigned int i, const unsigned int j, @@ -720,7 +765,7 @@ namespace TrilinosWrappers * operation and you should * always take care where to * call this function. As in - * the &p SparseMatrix + * the deal.II sparse matrix * class, we throw an exception * if the respective entry * doesn't exist in the @@ -731,10 +776,6 @@ namespace TrilinosWrappers * 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; @@ -810,7 +851,7 @@ namespace TrilinosWrappers * n=local_size(). */ std::pair - local_range () const; + local_range () const; /** * Return whether @p index is @@ -832,31 +873,34 @@ namespace TrilinosWrappers 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$. + * Return the + * l1-norm of + * the matrix, that is + * $|M|_1=\max_{\mathrm{all + * columns} j}\sum_{\mathrm{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$. + * Return the linfty-norm of the + * matrix, that is + * $|M|_\infty=\max_{\mathrm{all + * rows} i}\sum_{\mathrm{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) */ @@ -884,22 +928,20 @@ namespace TrilinosWrappers SparseMatrix & operator /= (const TrilinosScalar factor); /** - * Matrix-vector - * multiplication: let dst = - * M*src with M - * being this matrix. + * 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 + * 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 + * as was used for the matrix in + * case you work on a distributed + * memory architecture, using the + * interface in the * TrilinosWrappers::Vector * class. */ @@ -907,24 +949,23 @@ namespace TrilinosWrappers const Vector &src) const; /** - * Matrix-vector - * multiplication: let dst = + * Matrix-vector multiplication: + * let dst = * MT*src with * M being this - * matrix. This function does - * the same as vmult() but - * takes the transposed matrix. + * 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 + * 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 + * as was used for the matrix in + * case you work on a distributed + * memory architecture, using the * interface of the * TrilinosWrappers::Vector * class. @@ -942,13 +983,12 @@ namespace TrilinosWrappers * Source and destination must * not be the same vector. * - * Note that both vectors have - * to be distributed vectors + * 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 + * as was used for the matrix in + * case you work on a distributed + * memory architecture, using the * interface of the * TrilinosWrappers::Vector * class. @@ -960,22 +1000,21 @@ namespace TrilinosWrappers * 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. + * 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 + * 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 + * as was used for the matrix in + * case you work on a distributed + * memory architecture, using the * interface of the * TrilinosWrappers::Vector * class. @@ -984,23 +1023,22 @@ namespace TrilinosWrappers 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. + * 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 + * Obviously, the matrix needs to + * be quadratic for this * operation. * * The implementation of this @@ -1008,11 +1046,10 @@ namespace TrilinosWrappers * 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. + * 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 diff --git a/deal.II/lac/include/lac/trilinos_vector.h b/deal.II/lac/include/lac/trilinos_vector.h index 7ddcf36136..41ad01a263 100755 --- a/deal.II/lac/include/lac/trilinos_vector.h +++ b/deal.II/lac/include/lac/trilinos_vector.h @@ -35,6 +35,8 @@ # endif # include "Epetra_FEVector.h" # include "Epetra_Map.h" +# include "Epetra_LocalMap.h" +# include "Epetra_MultiVector.h" DEAL_II_NAMESPACE_OPEN @@ -50,6 +52,7 @@ namespace TrilinosWrappers { // forward declaration class Vector; + class LocalizedVector; /** * @cond internal @@ -322,7 +325,7 @@ namespace TrilinosWrappers * @see @ref SoftwareTrilinos * @author Martin Kronbichler, Wolfgang Bangerth, 2008 */ - class Vector + class Vector : public Subscriptor { public: /** @@ -992,6 +995,97 @@ namespace TrilinosWrappers }; +/** + * This class implements a localized version of a Trilinos vector. The + * purpose of this class is to provide a copy interface from the + * possibly parallel Vector class to a local vector on each processor. + * + * @ingroup TrilinosWrappers + * @ingroup Vectors + * @see @ref SoftwareTrilinos + * @author Martin Kronbichler, 2008 + */ + class LocalizedVector : Vector + { + public: + /** + * Default constructor that + * generates an empty (zero size) + * vector. The function + * reinit() will have to + * give the vector the correct + * size. + */ + LocalizedVector (); + + /** + * This constructor takes an + * Epetra_LocalMap that already knows + * the number of elements in the vector. + */ + LocalizedVector (const Epetra_LocalMap &InputMap); + + /** + * This constructor takes a + * (possibly parallel) Trilinos + * Vector and generates a + * localized version of the whole + * content on each processor. + */ + LocalizedVector (const Vector &V); + + /** + * This constructor generates a + * LocalizedVector based on a + * LocalizedVector input. + */ + LocalizedVector (const LocalizedVector &V); + + /** + * Reinit function. Takes the + * information of a parallel + * Trilinos vector and copies + * everything to this local + * version. + */ + void reinit (const Vector &V); + + /** + * Reinit function. Takes the + * information of a + * LocalizedVector and copies + * everything to the calling + * vector. + */ + void reinit (const LocalizedVector &V); + + /** + * Sets the left hand argument to + * the (parallel) Trilinos + * Vector. Equivalent to the @p + * reinit function. + */ + LocalizedVector & + operator = (const Vector &V); + + /** + * Copy operator. Copies both the + * dimension and the content in + * the right hand argument. + */ + LocalizedVector & + operator = (const LocalizedVector &V); + + private: + Epetra_LocalMap map; + + public: + std::auto_ptr vector; + }; + + + + // ------------------- inline and template functions -------------- diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index 1dcabdf37d..e3d336e386 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -184,10 +184,19 @@ Vector::Vector (const TrilinosWrappers::Vector &v) val = new Number[max_vec_size]; Assert (val != 0, ExcOutOfMemory()); + // Copy the distributed vector to + // a local one at all + // processors. TODO: There could + // be a better solution than + // this, but it has not yet been + // found. + TrilinosWrappers::LocalizedVector localized_vector (v); + // get a representation of the vector // and copy it TrilinosScalar **start_ptr; - int ierr = v.vector->ExtractView (&start_ptr); + + int ierr = localized_vector.vector->ExtractView (&start_ptr); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); std::copy (start_ptr[0], start_ptr[0]+vec_size, begin()); diff --git a/deal.II/lac/source/trilinos_sparse_matrix.cc b/deal.II/lac/source/trilinos_sparse_matrix.cc index 50b3cd382f..6e70110c5b 100755 --- a/deal.II/lac/source/trilinos_sparse_matrix.cc +++ b/deal.II/lac/source/trilinos_sparse_matrix.cc @@ -110,7 +110,7 @@ namespace TrilinosWrappers matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, (int*)const_cast(&(n_entries_per_row[0])), - true))) + true))) {} SparseMatrix::SparseMatrix (const Epetra_Map &InputRowMap, @@ -135,7 +135,7 @@ namespace TrilinosWrappers matrix (std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, col_map, (int*)const_cast(&(n_entries_per_row[0])), - true))) + true))) {} @@ -161,21 +161,19 @@ namespace TrilinosWrappers unsigned int n_rows = sparsity_pattern.n_rows(); // TODO: As of now, assume that the - // sparsity pattern only sits at the - // zeroth processor (completely), let - // this determine the Trilinos - // sparsity pattern on that - // processor, and only broadcast the - // pattern afterwards. + // sparsity pattern sits at the all + // processors (completely), and let + // each processor set its rows. Since + // this is wasteful, a better solution + // should be found in the future. Assert (matrix->NumGlobalRows() == (int)sparsity_pattern.n_rows(), ExcDimensionMismatch (matrix->NumGlobalRows(), sparsity_pattern.n_rows())); - // Trilinos seems to have a bug for - // rectangular matrices at this point, - // so do not check for consistent - // column numbers here. - // + // Trilinos has a bug for rectangular + // matrices at this point, so do not + // check for consistent column numbers + // here. // // this bug is filed in the Sandia // bugzilla under #4123 and should be @@ -261,7 +259,7 @@ namespace TrilinosWrappers // correct values. matrix = std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, - &n_entries_per_row[0], false)); + &n_entries_per_row[0], true)); reinit (sparsity_pattern); } @@ -330,7 +328,7 @@ namespace TrilinosWrappers // correct values. matrix = std::auto_ptr (new Epetra_FECrsMatrix(Copy, row_map, - &n_entries_per_row[0], false)); + &n_entries_per_row[0], true)); std::vector values; std::vector row_indices; @@ -490,8 +488,6 @@ namespace TrilinosWrappers const_cast(&value), &trilinos_j); - if (ierr > 0) - std::cout << ierr << " " << m() << " " << n() << std::endl; AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j)); AssertThrow (ierr == 0, ExcTrilinosError(ierr)); } @@ -938,7 +934,7 @@ namespace TrilinosWrappers { vmult (dst, x); dst -= b; - dst *= -1; + dst *= -1.; return dst.l2_norm(); } @@ -1023,9 +1019,9 @@ namespace TrilinosWrappers bool SparseMatrix::is_symmetric (const double tolerance) { - //bool truth; - if (tolerance == 0) - Assert (false, ExcNotImplemented()); + (void)tolerance; + + Assert (false, ExcNotImplemented()); return false; } @@ -1035,8 +1031,6 @@ namespace TrilinosWrappers bool SparseMatrix::is_hermitian () { - //bool truth; - Assert (false, ExcNotImplemented()); return false; } diff --git a/deal.II/lac/source/trilinos_vector.cc b/deal.II/lac/source/trilinos_vector.cc index 7a0233ab32..3a8243dae9 100755 --- a/deal.II/lac/source/trilinos_vector.cc +++ b/deal.II/lac/source/trilinos_vector.cc @@ -73,6 +73,7 @@ namespace TrilinosWrappers Vector::Vector (const Vector &v, const bool fast) : + Subscriptor(), map (v.map), last_action (Insert), vector(std::auto_ptr @@ -881,6 +882,78 @@ namespace TrilinosWrappers return 0; } + + + LocalizedVector::LocalizedVector () + : +#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_MultiVector(map,1))) + {} + + + + LocalizedVector::LocalizedVector (const Epetra_LocalMap &InputMap) + : + map (InputMap), + vector (std::auto_ptr + (new Epetra_MultiVector(map,1))) + {} + + + + LocalizedVector::LocalizedVector (const Vector &v) + : +#ifdef DEAL_II_COMPILER_SUPPORTS_MPI + map (0, 0, Epetra_MpiComm(MPI_COMM_WORLD)) +#else + map (0, 0, Epetra_SerialComm()) +#endif + { + reinit (v); + } + + + + void + LocalizedVector::reinit (const Vector &v) + { + map = Epetra_LocalMap (v.size(),0,v.vector->Comm()); + vector = std::auto_ptr (new Epetra_MultiVector (map,1,false)); + *vector = *v.vector; + } + + + + void + LocalizedVector::reinit (const LocalizedVector &v) + { + map = v.map; + *vector = *v.vector; + } + + + + LocalizedVector & + LocalizedVector::operator = (const Vector &v) + { + reinit (v); + return *this; + } + + + + LocalizedVector & + LocalizedVector::operator = (const LocalizedVector &v) + { + reinit (v); + return *this; + } + } DEAL_II_NAMESPACE_CLOSE -- 2.39.5