]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
The Trilinos sparse matrix class can now insert elements using set(i,j) as long as...
authorkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 6 Nov 2008 09:52:15 +0000 (09:52 +0000)
committerkronbichler <kronbichler@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 6 Nov 2008 09:52:15 +0000 (09:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@17483 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 6520591907e89b2f5a4b1dcc8e52bad9a392da11..b557e10e7666234142a08c9066e694e339b8054a 100755 (executable)
@@ -836,15 +836,17 @@ namespace TrilinosWrappers
                                         * Set the element (<i>i,j</i>)
                                         * to @p value.
                                        *
-                                       * Just as the respective call in
-                                       * deal.II SparseMatrix<Number>
-                                       * 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
-                                       * <tt>value</tt> is not a finite
+                                       * This function is able to insert
+                                       * new elements into the matrix as
+                                       * long as compress() has not been
+                                       * called, so the sparsity pattern
+                                       * will be extended. When compress()
+                                       * is called for the first time, then
+                                       * this is no longer possible and an
+                                       * insertion of elements at positions
+                                       * which have not been initialized
+                                       * will throw an exception. Moreover,
+                                       * if <tt>value</tt> is not a finite
                                        * number an exception is thrown.
                                        */
       void set (const unsigned int i,
index 8feeac441b159e691436d6539e6992f02e7bac8b..dfc1be2f26b0ba0021f4596e3e0dc3ddfccd20a3 100755 (executable)
@@ -630,9 +630,21 @@ namespace TrilinosWrappers
     int trilinos_i = i;
     int trilinos_j = j;
 
-    const int ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
-                                                 const_cast<double*>(&value), 
-                                                 &trilinos_j);
+    int ierr;
+
+                                       // If the matrix is not yet filled,
+                                       // we can insert new entries into
+                                       // the matrix using this
+                                       // command. Otherwise, we're just
+                                       // able to replace existing entries.
+    if (matrix->Filled() == false)
+      ierr = matrix->InsertGlobalValues (trilinos_i, 1,
+                                        const_cast<double*>(&value), 
+                                        &trilinos_j);
+    else
+      ierr = matrix->ReplaceGlobalValues (trilinos_i, 1,
+                                         const_cast<double*>(&value), 
+                                         &trilinos_j);
 
     AssertThrow (ierr <= 0, ExcAccessToNonPresentElement(i,j));
     AssertThrow (ierr == 0, ExcTrilinosError(ierr));
@@ -762,7 +774,7 @@ namespace TrilinosWrappers
                                      // already is transformed to
                                      // local indices.
        if (matrix->Filled() == false)
-         matrix->FillComplete(col_map, row_map, true);
+         matrix->GlobalAssemble(col_map, row_map, true);
 
                                      // Prepare pointers for extraction
                                      // of a view of the row.
@@ -785,10 +797,10 @@ namespace TrilinosWrappers
                                      // look for the value, and then
                                      // finally get it.
       
-       int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
+       int* el_find = std::find(col_indices, col_indices + nnz_present,
                                 trilinos_j);
-      
-       int el_index = (int)(el_find - col_indices);
+
+       int local_col_index = (int)(el_find - col_indices);
 
                                        // This is actually the only
                                        // difference to the el(i,j)
@@ -798,13 +810,12 @@ namespace TrilinosWrappers
                                        // returning zero for an
                                        // element that is not present
                                        // in the sparsity pattern.
-       if (!el_find)
+       if (local_col_index == nnz_present)
          {
            Assert (false, ExcInvalidIndex (i,j));
          }
        else
-         value = values[el_index];
-
+         value = values[local_col_index];
       }
 
     return value;
@@ -841,7 +852,7 @@ namespace TrilinosWrappers
                                      // already is transformed to
                                      // local indices.
       if (!matrix->Filled())
-       matrix->FillComplete(col_map, row_map, true);
+       matrix->GlobalAssemble(col_map, row_map, true);
 
                                      // Prepare pointers for extraction
                                      // of a view of the row.
@@ -863,18 +874,23 @@ namespace TrilinosWrappers
                                      // Search the index where we
                                      // look for the value, and then
                                      // finally get it.
-      
-      int* el_find = std::find(&col_indices[0],&col_indices[0] + nnz_present,
-                               trilinos_j);
-      
-      int el_index = (int)(el_find - col_indices);
+      int* el_find = std::find(col_indices, col_indices + nnz_present,
+                              trilinos_j);
+
+      int local_col_index = (int)(el_find - col_indices);
 
-      if (!el_find)
+
+                                       // This is actually the only
+                                       // difference to the () function
+                                       // querying (i,j), where we throw an
+                                       // exception instead of just
+                                       // returning zero for an element
+                                       // that is not present in the
+                                       // sparsity pattern.
+      if (local_col_index == nnz_present)
        value = 0;
       else
-       {
-         value = values[el_index];
-       }
+       value = values[local_col_index];
     }
 
     return value;

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.