From: wolf Date: Mon, 2 May 2005 22:04:24 +0000 (+0000) Subject: Implement apply_boundary_values also for parallel petsc matrices. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=6dccd67766eb0b0ce4f48ece0d616b68ddf2685d;p=dealii-svn.git Implement apply_boundary_values also for parallel petsc matrices. git-svn-id: https://svn.dealii.org/trunk@10596 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/matrices.h b/deal.II/deal.II/include/numerics/matrices.h index 34498542ae..4c0f8fe612 100644 --- a/deal.II/deal.II/include/numerics/matrices.h +++ b/deal.II/deal.II/include/numerics/matrices.h @@ -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 &boundary_values, + PETScWrappers::MPI::SparseMatrix &matrix, + PETScWrappers::MPI::Vector &solution, + PETScWrappers::MPI::Vector &right_hand_side, + const bool eliminate_columns = true); #endif /** diff --git a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc index d5227f3728..69463f6190 100644 --- a/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc +++ b/deal.II/deal.II/source/numerics/matrices.all_dimensions.cc @@ -20,7 +20,9 @@ #include #ifdef DEAL_II_USE_PETSC +# include # include +# include # include #endif @@ -556,7 +558,7 @@ apply_boundary_values (const std::map &boundary_values, { first_nonzero_diagonal_entry = matrix.diag_element(i); break; - }; + } std::map::const_iterator @@ -592,19 +594,9 @@ apply_boundary_values (const std::map &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 &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 &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::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