]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Corrected some minor bugs in Trilinos wrappers. Added clear_row() functionality to...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Aug 2008 15:36:33 +0000 (15:36 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 26 Aug 2008 15:36:33 +0000 (15:36 +0000)
git-svn-id: https://svn.dealii.org/trunk@16669 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/source/numerics/matrices.all_dimensions.cc
deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/include/lac/trilinos_vector.h
deal.II/lac/source/trilinos_precondition_amg.cc
deal.II/lac/source/trilinos_sparse_matrix.cc
deal.II/lac/source/trilinos_vector.cc

index 78da5349246314fb3f83c1ad502c24f494454dd6..3034f8123651288e26e30297cdbde01880a1ae3a 100644 (file)
@@ -746,18 +746,11 @@ namespace TrilinosWrappers
 
                                      // then eliminate these rows and set
                                      // their diagonal entry to what we have
-                                     // determined above. note that for trilinos
-                                     // matrices interleaving read with write
-                                     // operations is very expensive. thus, we
-                                     // here always replace the diagonal
-                                     // element, rather than first checking
-                                     // whether it is nonzero and in that case
-                                     // preserving it. this is different from
-                                     // the case of deal.II sparse matrices
-                                     // treated in the other functions.
-//TODO: clear_row is not currently implemented for Trilinos    
-    Assert (false, ExcInternalError());
-//    matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
+                                     // determined above. if the value already
+                                     // is nonzero, it will be preserved, 
+                                     // in accordance with the basic 
+                                     // matrix classes in deal.II.
+    matrix.clear_rows (constrained_rows, average_nonzero_diagonal_entry);
 
                                      // the next thing is to set right hand
                                      // side to the wanted value. there's one
index 3285cc6a638a4c4071a85059fdceef8d41b36dd4..8d3bdfd5d4c7f78372c03a85c8bf51e0bdbb50eb 100755 (executable)
@@ -477,6 +477,21 @@ namespace TrilinosWrappers
                                         */
       void clear ();
 
+                                       /**
+                                        * 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 ();
+
                                        /**
                                         * This operator assigns a scalar to a
                                         * matrix. Since this does usually not
@@ -545,7 +560,10 @@ namespace TrilinosWrappers
                                         * pattern, though (but retains
                                         * the allocated memory in case
                                         * new entries are again added
-                                        * later).
+                                        * later). Note that this
+                                       * is a global operation, so this
+                                       * needs to be done on all 
+                                       * MPI processes.
                                         *
                                         * This operation is used in
                                         * eliminating constraints (e.g. due to
@@ -554,7 +572,7 @@ namespace TrilinosWrappers
                                         * the matrix without having to read
                                         * entries (such as the locations of
                                         * non-zero elements) from it --
-                                        * without this operation, removing
+                                        * without this o peration, removing
                                         * constraints on parallel matrices is
                                         * a rather complicated procedure.
                                         *
@@ -563,8 +581,8 @@ namespace TrilinosWrappers
                                         * 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);
+      void clear_row (const unsigned int   row,
+                      const TrilinosScalar new_diag_value = 0);
 
                                        /**
                                         * Same as clear_row(), except that it
@@ -579,24 +597,9 @@ namespace TrilinosWrappers
                                         * 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 ();
-      
+      void clear_rows (const std::vector<unsigned int> &rows,
+                       const TrilinosScalar             new_diag_value = 0);
+
                                        /**
                                         * Return the value of the entry
                                         * (<i>i,j</i>).  This may be an
@@ -998,6 +1001,15 @@ namespace TrilinosWrappers
                                         */
       void write_ascii ();
       
+                                    /**
+                                     * Print the matrix to the given
+                                     * stream, using the format
+                                     * <tt>(line,col) value</tt>,
+                                     * i.e. one nonzero entry of the
+                                     * matrix per line.
+                                     */
+    void print (std::ostream &out) const;
+    
                                         // TODO: Write an overloading
                                         // of the operator << for output.
                                         // Since the underlying Trilinos 
index 7b9e151a6d11fdb7afb689600a0ec609554b30d8..33af1a0c38035bba0cdca20bfb9277ce2f82df97 100755 (executable)
@@ -770,6 +770,13 @@ namespace TrilinosWrappers
       void ratio (const Vector &a,
                   const Vector &b);
 
+                                    /**
+                                     *  Output of vector in user-defined
+                                     *  format in analogy to the 
+                                     *  dealii::Vector<number> class.
+                                     */
+    void print (const char* format = 0) const;
+
                                        /**
                                         * Print to a
                                         * stream. @p precision denotes
@@ -900,7 +907,7 @@ namespace TrilinosWrappers
  *
  * @ingroup TrilinosWrappers
  * @relates TrilinosWrappers::Vector
- * @author Wolfgang Bangerth, 2004
+ * @author Martin Kronbichler, Wolfgang Bangerth, 2008
  */
   inline
   void swap (Vector &u, Vector &v)
index 5791c8a2fce8a9097fceda02488cfa7243df9b8e..7d287c230d621730529ee01a95276ec11ae7fbb7 100755 (executable)
@@ -175,7 +175,7 @@ namespace TrilinosWrappers
     
     Epetra_Vector LHS (View, multigrid_operator->OperatorDomainMap(),
                       dst.begin());
-    Epetra_Vector RHS (View, multigrid_operator->OperatorDomainMap(),
+    Epetra_Vector RHS (View, multigrid_operator->OperatorRangeMap(),
                       const_cast<double*>(src.begin()));
   
     const int res = multigrid_operator->ApplyInverse (RHS, LHS);
index d0b2054fd53a3f55d4fe51b7b3df95c27c24b985..3af6a7a2cfd98c6156ff37ac219e7e5167c1ab61 100755 (executable)
@@ -195,14 +195,10 @@ namespace TrilinosWrappers
            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.
-                                    //
-                                    // this bug is filed in the Sandia
-                                    // bugzilla under #4123 and should be
-                                    // fixed for version 9.0
+                                // Trilinos seems to have a bug for
+                                // rectangular matrices at this point,
+                                // so do not check for consistent 
+                                // column numbers here.
 //    Assert (matrix->NumGlobalCols() == (int)sparsity_pattern.n_cols(),
 //         ExcDimensionMismatch (matrix->NumGlobalCols(),
 //                               sparsity_pattern.n_cols()));
@@ -331,7 +327,7 @@ namespace TrilinosWrappers
        for (dealii::SparseMatrix<double>::const_iterator  
              p  = deal_ii_sparse_matrix.begin(row);
              p != deal_ii_sparse_matrix.end(row); ++p)
-         if (std::abs(p->value()) > drop_tolerance)
+         if (std::fabs(p->value()) > drop_tolerance)
            {
              row_indices[index] = p->column();
              values[index]      = p->value();
@@ -461,7 +457,7 @@ namespace TrilinosWrappers
                                      // case to be consistent with the MPI
                                      // communication model (see the
                                      // comments in the documentation of
-                                     // TrilinosWrappers::MPI::Vector), but we
+                                     // TrilinosWrappers::Vector), but we
                                      // can save some work if the addend is
                                      // zero
     if (value == 0)
@@ -480,6 +476,52 @@ namespace TrilinosWrappers
 
 
 
+  void
+  SparseMatrix::clear_row (const unsigned int   row,
+                          const TrilinosScalar new_diag_value)
+  {
+    Assert (matrix->Filled()==true,
+           ExcMessage("Matrix must be compressed before invoking clear_row."));
+
+                                    // Only do this on the rows
+                                    // owned locally on this processor.
+    int local_row = matrix->LRID(row);
+    if (local_row >= 0)
+      {
+       TrilinosScalar *values;
+       int *col_indices;
+       int num_entries;
+       const int ierr = matrix->ExtractMyRowView(local_row, num_entries,
+                                                 values, col_indices);
+       
+       Assert (ierr == 0,
+               ExcTrilinosError(ierr));
+       
+       int* diag_find = std::find(col_indices,col_indices+num_entries, 
+                                  local_row);
+       int diag_index = (int)(diag_find - col_indices);
+
+       for (int j=0; j<num_entries; ++j)
+         if (diag_index != col_indices[j])
+          values[j] = 0.;
+
+       if (diag_find && std::fabs(values[diag_index]) > 0.)
+         values[diag_index] = new_diag_value;
+      }
+  }
+
+
+
+  void
+  SparseMatrix::clear_rows (const std::vector<unsigned int> &rows,
+                           const TrilinosScalar             new_diag_value)
+  {
+    for (unsigned int row=0; row<rows.size(); ++row)
+      clear_row(rows[row], new_diag_value);
+  }
+
+
+
   TrilinosScalar
   SparseMatrix::el (const unsigned int i,
                    const unsigned int j) const
@@ -488,7 +530,7 @@ namespace TrilinosWrappers
                                       // 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.
@@ -525,16 +567,17 @@ namespace TrilinosWrappers
                                      // 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* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+                              trilinos_j);
+      
+      int el_index = (int)(el_find - col_indices);
 
-      int position;
-      if (!index)
+      if (!el_find)
        value = 0;
       else
        {
-         position = (int)(index - &(col_indices[0]));
-         value = values[position];
+         value = values[el_index];
        }
     }
 
@@ -548,8 +591,6 @@ namespace TrilinosWrappers
   {
     Assert (m() == n(), ExcNotQuadratic());
 
-                                     // this doesn't seem to work any
-                                     // different than any other element
     return el(i,i);
   }
 
@@ -589,7 +630,7 @@ namespace TrilinosWrappers
   {
     int begin, end;
     begin = matrix->RowMap().MinMyGID();
-    end = matrix->RowMap().MaxMyGID();
+    end = matrix->RowMap().MaxMyGID()+1;
     
     return std::make_pair (begin, end);
   }
@@ -827,12 +868,36 @@ namespace TrilinosWrappers
     return false;
   }  
 
+
+
   void
   SparseMatrix::write_ascii ()
   {
     Assert (false, ExcNotImplemented());
   }
 
+
+
+                                  // As of now, no particularly neat
+                                  // ouput is generated in case of 
+                                  // multiple processors.
+  void SparseMatrix::print (std::ostream &out) const
+  {
+    double * values;
+    int * indices;
+    int num_entries;
+  
+    for (int i=0; i<matrix->NumMyRows(); ++i)
+      {
+        matrix->ExtractMyRowView (i, num_entries, values, indices);
+       for (int j=0; j<num_entries; ++j)
+         out << "(" << i << "," << indices[matrix->GRID(j)] << ") " 
+             << values[j] << std::endl;
+      }
+  
+    AssertThrow (out, ExcIO());
+  }
+
 }
 
 DEAL_II_NAMESPACE_CLOSE
index 0c952be3531f75de262c9cff893d5c9f8726ebb6..09d6bfd169d0cbe5732c18cbfef01f9d4dc6f384 100755 (executable)
@@ -207,7 +207,7 @@ namespace TrilinosWrappers
   {
     int begin, end;
     begin = vector->Map().MinMyGID();
-    end = vector->Map().MaxMyGID();
+    end = vector->Map().MaxMyGID()+1;
     return std::make_pair (begin, end);
   }
 
@@ -520,17 +520,14 @@ namespace TrilinosWrappers
            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]);
+    unsigned int n_local = local_size();
+    int ierr;
 
-    delete[] index;
-    
-    AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+    for (unsigned int i=0; i<n_local; i++)
+      {
+       ierr = vector->SumIntoMyValue(i,0,s);
+        AssertThrow (ierr == 0, ExcTrilinosError(ierr));
+      }
   }
 
 
@@ -545,14 +542,14 @@ namespace TrilinosWrappers
 
   void
   Vector::add (const TrilinosScalar a,
-                   const Vector     &v)
+              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);
+    const int ierr = vector->Update(a, *(v.vector), 1.);
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
 
@@ -560,9 +557,9 @@ namespace TrilinosWrappers
 
   void
   Vector::add (const TrilinosScalar a,
-                   const Vector &v,
-                   const TrilinosScalar b,
-                   const Vector &w)
+              const Vector        &v,
+              const TrilinosScalar b,
+              const Vector        &w)
   {
 
     Assert (numbers::is_finite(a),
@@ -572,7 +569,7 @@ namespace TrilinosWrappers
            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);
+    const int ierr = vector->Update(a, *(v.vector), b, *(w.vector), 1.);
 
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
@@ -581,14 +578,14 @@ namespace TrilinosWrappers
 
   void
   Vector::sadd (const TrilinosScalar s,
-                    const Vector &v)
+               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);
+    const int ierr = vector->Update(1., *(v.vector), s);
 
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
   }
@@ -597,8 +594,8 @@ namespace TrilinosWrappers
 
   void
   Vector::sadd (const TrilinosScalar s,
-                    const TrilinosScalar a,
-                    const Vector     &v)
+               const TrilinosScalar a,
+               const Vector        &v)
   {
 
     Assert (numbers::is_finite(s),
@@ -617,10 +614,10 @@ namespace TrilinosWrappers
 
   void
   Vector::sadd (const TrilinosScalar s,
-                    const TrilinosScalar a,
-                    const Vector     &v,
-                    const TrilinosScalar b,
-                    const Vector     &w)
+               const TrilinosScalar a,
+               const Vector        &v,
+               const TrilinosScalar b,
+               const Vector        &w)
   {
 
     Assert (numbers::is_finite(s),
@@ -642,12 +639,12 @@ namespace TrilinosWrappers
 
   void
   Vector::sadd (const TrilinosScalar s,
-                    const TrilinosScalar a,
-                    const Vector     &v,
-                    const TrilinosScalar b,
-                    const Vector     &w,
-                    const TrilinosScalar c,
-                    const Vector     &x)
+               const TrilinosScalar a,
+               const Vector        &v,
+               const TrilinosScalar b,
+               const Vector        &w,
+               const TrilinosScalar c,
+               const Vector        &x)
   {
 
     Assert (numbers::is_finite(s),
@@ -663,12 +660,13 @@ namespace TrilinosWrappers
            ExcMessage("The given value is not finite but "
                       "either infinite or Not A Number (NaN)"));
 
-                                     // Update member can only input two other vectors so
+                                     // 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);
+
+    const int jerr = vector->Update(c, *(x.vector), 1.);
     AssertThrow (jerr == 0, ExcTrilinosError(jerr));
   }
 
@@ -692,11 +690,10 @@ namespace TrilinosWrappers
            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));
+    *vector = *v.vector;
+    map = v.map;
+    
+    *this *= a;
   }
 
 
@@ -715,15 +712,12 @@ namespace TrilinosWrappers
            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));
+    Assert (v.size() == w.size(),
+            ExcDimensionMismatch (v.size(), w.size()));
 
+    *vector = *v.vector;
+    map = v.map;
+    sadd (a, b, w);
   }
 
 
@@ -732,8 +726,8 @@ namespace TrilinosWrappers
   Vector::ratio (const Vector &v,
                 const Vector &w)
   {
-    Assert (size() == v.size(),
-            ExcDimensionMismatch (size(), v.size()));
+    Assert (v.size() == w.size(),
+            ExcDimensionMismatch (v.size(), w.size()));
 
     Assert (size() == w.size(),
             ExcDimensionMismatch (size(), w.size()));
@@ -745,6 +739,28 @@ namespace TrilinosWrappers
   }
 
 
+  
+                                    // TODO: up to now only local data
+                                     // printed out! Find a way to neatly
+                                    // output distributed data...
+  void
+  Vector::print (const char *format) const
+  {
+    Assert (vector->GlobalLength()!=0, ExcEmptyObject());
+    
+    for (unsigned int j=0; j<size(); ++j)
+      {
+        double t = (*vector)[0][j];
+
+       if (format != 0)
+          std::printf (format, t);
+        else
+          std::printf (" %5.2f", double(t));
+      }
+    std::printf ("\n");
+  }
+
+
 
   void
   Vector::print (std::ostream      &out,
@@ -757,7 +773,8 @@ namespace TrilinosWrappers
                                      // get a representation of the vector and
                                      // loop over all the elements 
                                      // TODO: up to now only local data
-                                     // printed out!
+                                     // printed out! Find a way to neatly
+                                    // output distributed data...
     TrilinosScalar *val;
     int leading_dimension;
     int ierr = vector->ExtractView (&val, &leading_dimension);
@@ -790,10 +807,10 @@ namespace TrilinosWrappers
                                     // 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;
+    Vector *p_v = &v, *p_this = this;
+    Vector* tmp = p_v;
+    p_v = p_this;
+    p_this = tmp;
   }
 
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.