]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Implement apply_boundary_values also for parallel petsc matrices.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 May 2005 22:04:24 +0000 (22:04 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 2 May 2005 22:04:24 +0000 (22:04 +0000)
git-svn-id: https://svn.dealii.org/trunk@10596 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/matrices.h
deal.II/deal.II/source/numerics/matrices.all_dimensions.cc

index 34498542ae9ef75b7aada9b35e646ddbdf333154..4c0f8fe6120a18e321c7fc58702b6e94c5e22b0c 100644 (file)
@@ -41,6 +41,11 @@ namespace PETScWrappers
 {
   class SparseMatrix;
   class Vector;
+  namespace MPI
+  {
+    class SparseMatrix;
+    class Vector;
+  }
 }
 #endif
 
@@ -834,6 +839,13 @@ class MatrixTools : public MatrixCreator
                           PETScWrappers::Vector  &solution,
                           PETScWrappers::Vector  &right_hand_side,
                           const bool             eliminate_columns = true);
+
+    static void
+    apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                          PETScWrappers::MPI::SparseMatrix  &matrix,
+                          PETScWrappers::MPI::Vector  &solution,
+                          PETScWrappers::MPI::Vector  &right_hand_side,
+                          const bool             eliminate_columns = true);
 #endif
 
                                      /**
index d5227f37283767daf6e21de7a25ad5b54ae26182..69463f6190720397eabdaaa3594b139752d79245 100644 (file)
@@ -20,7 +20,9 @@
 #include <lac/block_sparse_matrix.h>
 
 #ifdef DEAL_II_USE_PETSC
+#  include <lac/petsc_parallel_sparse_matrix.h>
 #  include <lac/petsc_sparse_matrix.h>
+#  include <lac/petsc_parallel_vector.h>
 #  include <lac/petsc_vector.h>
 #endif
 
@@ -556,7 +558,7 @@ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
       {
        first_nonzero_diagonal_entry = matrix.diag_element(i);
        break;
-      };
+      }
 
   
   std::map<unsigned int,double>::const_iterator
@@ -592,19 +594,9 @@ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
       }
       
 
-                                      // set right hand side to
-                                      // wanted value: if main diagonal
-                                      // entry nonzero, don't touch it
-                                      // and scale rhs accordingly. If
-                                      // zero, take the first main
-                                      // diagonal entry we can find, or
-                                      // one if no nonzero main diagonal
-                                      // element exists. Normally, however,
-                                      // the main diagonal entry should
-                                      // not be zero.
-                                      //
-                                      // store the new rhs entry to make the
-                                      // gauss step (when preserving the
+                                      // set right hand side to wanted value;
+                                      // also store the new rhs entry to make
+                                      // the gauss step (when preserving the
                                       // symmetry of the matrix) more
                                       // efficient
                                        //
@@ -639,8 +631,139 @@ apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
                                       // preset solution vector
       solution(dof_number) = dof->second;
     }
+
+  matrix.compress ();
+  solution.compress ();
+  right_hand_side.compress ();
 }
 
+
+
+void
+MatrixTools::
+apply_boundary_values (const std::map<unsigned int,double> &boundary_values,
+                       PETScWrappers::MPI::SparseMatrix   &matrix,
+                       PETScWrappers::MPI::Vector   &solution,
+                       PETScWrappers::MPI::Vector   &right_hand_side,
+                       const bool        preserve_symmetry)
+{
+                                   // this function works almost exactly as
+                                   // the one above, with the only exception
+                                   // that we have to make sure that we only
+                                   // access elements that belong to the slice
+                                   // of the matrix that is available locally
+  Assert (matrix.n() == right_hand_side.size(),
+         ExcDimensionMismatch(matrix.n(), right_hand_side.size()));
+  Assert (matrix.n() == solution.size(),
+         ExcDimensionMismatch(matrix.n(), solution.size()));
+
+                                   // if no boundary values are to be applied
+                                  // simply return
+  if (boundary_values.size() == 0)
+    return;
+
+
+                                   // we have to read and write from this
+                                   // matrix (in this order). this will only
+                                   // work if we compress the matrix first,
+                                   // done here:
+  matrix.compress ();
+  
+  const unsigned int n_dofs = matrix.m();
+
+                                  // determine the first nonzero diagonal
+                                  // entry from within the part of the matrix
+                                  // that we can see. if we can't find such
+                                  // an entry, take one
+  PetscScalar first_nonzero_diagonal_entry = 1;
+  for (unsigned int i=matrix.local_range().first;
+       i<matrix.local_range().second; ++i)
+    if (matrix.diag_element(i) != 0)
+      {
+       first_nonzero_diagonal_entry = std::fabs(matrix.diag_element(i));
+       break;
+      }
+
+                                   // iterate over all fixed degrees of
+                                   // freedom that are within the local range
+                                   // of this matrix
+  std::map<unsigned int,double>::const_iterator
+    dof  = boundary_values.begin(),
+    endd = boundary_values.end();
+  for (; dof != endd; ++dof)
+    if ((dof->first >= matrix.local_range().first) &&
+        (dof->first < matrix.local_range().second))
+      {
+        const unsigned int dof_number = dof->first;
+      
+                                         // for each constrained dof:
+      
+                                         // set entries of this line
+                                         // to zero except for the diagonal
+                                         // entry.
+        {
+          PETScWrappers::MatrixBase::const_iterator
+            p = matrix.begin(dof_number),
+            e = matrix.end(dof_number);
+
+                                           // iterate over all elements of this
+                                           // row and set elements to zero
+                                           // except for the diagonal
+                                           // element. note that this is not
+                                           // exactly clean programming, since
+                                           // we change the matrix underneath,
+                                           // while we still keep working with
+                                           // the iterators into it
+          for (; p!=e; ++p)
+            if (p->column() != dof_number)
+              matrix.set (dof_number, p->column(), 0.);
+        }
+      
+
+                                      // set right hand side to wanted value;
+                                      // also store the new rhs entry to make
+                                      // the gauss step (when preserving the
+                                      // symmetry of the matrix) more
+                                      // efficient
+                                       //
+                                       // note that for petsc 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.
+      PetscScalar new_rhs;
+      matrix.set (dof_number, dof_number,
+                  first_nonzero_diagonal_entry);
+      new_rhs = dof->second * first_nonzero_diagonal_entry;
+      right_hand_side(dof_number) = new_rhs;
+
+                                      // if the user wants to have
+                                      // the symmetry of the matrix
+                                      // preserved, and if the
+                                      // sparsity pattern is
+                                      // symmetric, then do a Gauss
+                                      // elimination step with the
+                                      // present row
+      if (preserve_symmetry)
+       {
+          Assert (false, ExcNotImplemented());
+       }
+
+                                      // preset solution vector
+      solution(dof_number) = dof->second;
+    }
+
+  matrix.compress ();
+  solution.compress ();
+  right_hand_side.compress ();
+}
+
+
 #endif
 
 

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.