]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Improved the function that adds a Trilinos sparse matrix to another.
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 2 Dec 2008 12:37:53 +0000 (12:37 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 2 Dec 2008 12:37:53 +0000 (12:37 +0000)
git-svn-id: https://svn.dealii.org/trunk@17810 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/trilinos_sparse_matrix.h
deal.II/lac/source/trilinos_sparse_matrix.cc

index 2d25253f119a9b30995a6f2733c86e6a60132edc..21e9cf65e22e72996cdd9137d5c161e8918fef44 100755 (executable)
@@ -1140,11 +1140,16 @@ namespace TrilinosWrappers
        copy_from (const SparseMatrix &source);
 
                                       /**
-                                       * Add <tt>matrix</tt> scaled
-                                       * by <tt>factor</tt> to this
-                                       * matrix, i.e. the matrix
-                                       * <tt>factor*matrix</tt> is
-                                       * added to <tt>this</tt>.
+                                       * Add <tt>matrix</tt> scaled by
+                                       * <tt>factor</tt> to this matrix,
+                                       * i.e. the matrix
+                                       * <tt>factor*matrix</tt> is added to
+                                       * <tt>this</tt>. If the sparsity
+                                       * pattern of the calling matrix does
+                                       * not contain all the elements in
+                                       * the sparsity pattern of the input
+                                       * matrix, this function will throw
+                                       * an exception.
                                        */
       void add (const TrilinosScalar  factor,
                const SparseMatrix   &matrix);
index a71a854cf96fb5706d2fb14272068e53d51b6f78..3d4f2dedaf7912b3615ad6db62f9aeb1d4fea1be 100755 (executable)
@@ -1165,15 +1165,6 @@ namespace TrilinosWrappers
     Assert (rhs.m() == m(), ExcDimensionMismatch (rhs.m(), m()));
     Assert (rhs.n() == n(), ExcDimensionMismatch (rhs.n(), n()));
     
-                                 // I bet that there must be a
-                                  // better way to do this but it
-                                  // has not been found: currently,
-                                 // we simply go through each row
-                                 // of the argument matrix, copy
-                                 // it, scale it, and add it to
-                                 // the current matrix. that's
-                                 // probably not the most
-                                 // efficient way to do things.
     const std::pair<unsigned int, unsigned int>
       local_range = rhs.local_range();
 
@@ -1186,29 +1177,98 @@ namespace TrilinosWrappers
     
     std::vector<int>            column_indices (max_row_length);
     std::vector<TrilinosScalar> values (max_row_length);
+    int ierr;
 
                                   // If both matrices have been transformed
                                   // to local index space (in Trilinos
-                                  // speak: they are filled), we can
-                                  // extract MyRows instead of GlobalRows,
-                                  // which is faster.
-    if (matrix->Filled() == true && rhs.matrix->Filled() == true)
+                                  // speak: they are filled), we're having
+                                  // matrices based on the same indices
+                                  // with the same number of nonzeros
+                                  // (actually, we'd need sparsity pattern,
+                                  // but that is too expensive to check),
+                                  // we can extract views of the column
+                                  // data on both matrices and simply
+                                  // manipulate the values that are
+                                  // addressed by the pointers.
+    if (matrix->Filled() == true && 
+       rhs.matrix->Filled() == true &&
+       this->local_range() == local_range && 
+       matrix->NumMyNonzeros() == rhs.matrix->NumMyNonzeros())
+      for (unsigned int row=local_range.first;
+          row < local_range.second; ++row)
+       {
+         Assert (matrix->NumGlobalEntries(row) == 
+                 rhs.matrix->NumGlobalEntries(row),
+                 ExcDimensionMismatch(matrix->NumGlobalEntries(row),
+                                      rhs.matrix->NumGlobalEntries(row)));
+
+         const int row_local = matrix->RowMap().LID(row);
+         int n_entries, rhs_n_entries;
+         TrilinosScalar *value_ptr, *rhs_value_ptr;
+
+                                  // In debug mode, we want to check
+                                  // whether the indices really are the
+                                  // same in the calling matrix and the
+                                  // input matrix. The reason for doing
+                                  // this only in debug mode is that both
+                                  // extracting indices and comparing
+                                  // indices is relatively slow compared to
+                                  // just working with the values.
+#ifdef DEBUG
+         int *index_ptr, *rhs_index_ptr;
+         ierr = rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries, 
+                                              rhs_value_ptr, rhs_index_ptr);
+         Assert (ierr == 0, ExcTrilinosError(ierr));
+
+         ierr = matrix->ExtractMyRowView (row_local, n_entries, value_ptr,
+                                          index_ptr);
+         Assert (ierr == 0, ExcTrilinosError(ierr));
+#else
+         rhs.matrix->ExtractMyRowView (row_local, rhs_n_entries,rhs_value_ptr);
+         matrix->ExtractMyRowView (row_local, n_entries, value_ptr);
+#endif
+
+         AssertThrow (n_entries == rhs_n_entries,
+                      ExcDimensionMismatch (n_entries, rhs_n_entries));
+
+         for (int i=0; i<n_entries; ++i)
+           {
+             *value_ptr++ += *rhs_value_ptr++ * factor;
+#ifdef DEBUG
+             Assert (*index_ptr++ == *rhs_index_ptr++,
+                     ExcInternalError());
+#endif
+           }
+       }
+                                  // If we have different sparsity patterns
+                                  // (expressed by a different number of
+                                  // nonzero elements), we have to be more
+                                  // careful and extract a copy of the row
+                                  // data, multiply it by the factor and
+                                  // then add it to the matrix using the
+                                  // respective add() function.
+    else if (matrix->Filled() == true && rhs.matrix->Filled() == true &&
+            this->local_range() == local_range)
       for (unsigned int row=local_range.first;
           row < local_range.second; ++row)
        {
          const int row_local = matrix->RowMap().LID(row);
          int n_entries;
-         rhs.matrix->ExtractMyRowCopy (row_local, max_row_length,
-                                       n_entries,
-                                       &values[0],
-                                       &column_indices[0]);
+
+         ierr = rhs.matrix->ExtractMyRowCopy (row_local, max_row_length,
+                                              n_entries,
+                                              &values[0],
+                                              &column_indices[0]);
+         Assert (ierr == 0, ExcTrilinosError(ierr));
+
          for (int i=0; i<n_entries; ++i)
            values[i] *= factor;
 
          TrilinosScalar *value_ptr = &values[0];
 
-         matrix->SumIntoMyValues (row_local, n_entries, value_ptr,
-                                  &column_indices[0]);
+         ierr = matrix->SumIntoMyValues (row_local, n_entries, value_ptr,
+                                         &column_indices[0]);
+         Assert (ierr == 0, ExcTrilinosError(ierr));
        }
     else
       {
@@ -1216,18 +1276,16 @@ namespace TrilinosWrappers
             row < local_range.second; ++row)
          {
            int n_entries;
-           rhs.matrix->Epetra_CrsMatrix::ExtractGlobalRowCopy ((int)row, 
-                                                               max_row_length,
-                                                               n_entries,
-                                                               &values[0],
-                                                               &column_indices[0]);
+           ierr = rhs.matrix->Epetra_CrsMatrix::ExtractGlobalRowCopy 
+             ((int)row, max_row_length, n_entries, &values[0], &column_indices[0]);
+           Assert (ierr == 0, ExcTrilinosError(ierr));
+
            for (int i=0; i<n_entries; ++i)
              values[i] *= factor;
 
-           matrix->Epetra_CrsMatrix::SumIntoGlobalValues ((int)row, 
-                                                          n_entries, 
-                                                          &values[0],
-                                                          &column_indices[0]);
+           ierr = matrix->Epetra_CrsMatrix::SumIntoGlobalValues 
+             ((int)row, n_entries, &values[0], &column_indices[0]);
+           Assert (ierr == 0, ExcTrilinosError(ierr));
          }
        compress ();
       }

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.