From ce547d0f3f7c542d2cf23c899967c57dd4b733a0 Mon Sep 17 00:00:00 2001 From: young Date: Thu, 13 Feb 2014 13:02:23 +0000 Subject: [PATCH] Add constraint matrix patch written by Denis Davydov. git-svn-id: https://svn.dealii.org/branches/branch_petscscalar_complex@32479 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/cmake/config/template-arguments.in | 11 ++ .../include/deal.II/lac/constraint_matrix.h | 80 +++++++-- .../deal.II/lac/constraint_matrix.templates.h | 163 +++++++++--------- deal.II/source/lac/constraint_matrix.cc | 92 +++++++++- deal.II/source/lac/constraint_matrix.inst.in | 13 +- 5 files changed, 252 insertions(+), 107 deletions(-) diff --git a/deal.II/cmake/config/template-arguments.in b/deal.II/cmake/config/template-arguments.in index c95c0f21d6..74f9648d6d 100644 --- a/deal.II/cmake/config/template-arguments.in +++ b/deal.II/cmake/config/template-arguments.in @@ -79,6 +79,17 @@ EXTERNAL_PARALLEL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@ } + +REAL_EXTERNAL_SEQUENTIAL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_VECTOR@; + @DEAL_II_EXPAND_TRILINOS_BLOCKVECTOR@; + } + +COMPLEX_EXTERNAL_PARALLEL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@; + @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@; + @DEAL_II_EXPAND_PETSC_MPI_VECTOR@; + @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@ + } + VECTORS_WITH_MATRIX := { Vector; Vector ; Vector; diff --git a/deal.II/include/deal.II/lac/constraint_matrix.h b/deal.II/include/deal.II/lac/constraint_matrix.h index 75b91077b7..5dcfa922af 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.h @@ -812,12 +812,12 @@ public: * caller's site. There is no locking mechanism inside this method to * prevent data races. */ - template + template void distribute_local_to_global (const Vector &local_vector, const std::vector &local_dof_indices, VectorType &global_vector, - const FullMatrix &local_matrix) const; + const FullMatrix &local_matrix) const; /** * Enter a single value into a result vector, obeying constraints. @@ -910,6 +910,14 @@ public: distribute_local_to_global (const FullMatrix &local_matrix, const std::vector &local_dof_indices, MatrixType &global_matrix) const; + /** + * Same as above for complex. + */ + template + void + distribute_local_to_global (const FullMatrix > &local_matrix, + const std::vector &local_dof_indices, + MatrixType &global_matrix) const; /** * Does the same as the function above but can treat non quadratic matrices. @@ -945,6 +953,18 @@ public: VectorType &global_vector, bool use_inhomogeneities_for_rhs = false) const; + /** + * Same as above for complex. + */ + template + void + distribute_local_to_global (const FullMatrix > &local_matrix, + const Vector &local_vector, + const std::vector &local_dof_indices, + MatrixType &global_matrix, + VectorType &global_vector, + bool use_inhomogeneities_for_rhs = false) const; + /** * Do a similar operation as the distribute_local_to_global() function that * distributes writing entries into a matrix for constrained degrees of @@ -1312,10 +1332,10 @@ private: * This function actually implements the local_to_global function for * standard (non-block) matrices. */ - template + template void - distribute_local_to_global (const FullMatrix &local_matrix, - const Vector &local_vector, + distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, const std::vector &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, @@ -1326,10 +1346,10 @@ private: * This function actually implements the local_to_global function for block * matrices. */ - template + template void - distribute_local_to_global (const FullMatrix &local_matrix, - const Vector &local_vector, + distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, const std::vector &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, @@ -1385,12 +1405,13 @@ private: /** * Internal helper function for distribute_local_to_global function. */ - double + template + LocalType resolve_vector_entry (const size_type i, const internals::GlobalRowsFromLocal &global_rows, - const Vector &local_vector, + const Vector &local_vector, const std::vector &local_dof_indices, - const FullMatrix &local_matrix) const; + const FullMatrix &local_matrix) const; }; @@ -1720,6 +1741,23 @@ distribute_local_to_global (const FullMatrix &local_matrix, dealii::internal::bool2type::value>()); } +// Complex case (maybe can be merged with the function above) +template +inline +void +ConstraintMatrix:: +distribute_local_to_global (const FullMatrix > &local_matrix, + const std::vector &local_dof_indices, + MatrixType &global_matrix) const +{ + // create a dummy and hand on to the function actually implementing this + // feature in the cm.templates.h file. + Vector dummy(0); + distribute_local_to_global (local_matrix, dummy, local_dof_indices, + global_matrix, dummy, false, + dealii::internal::bool2type::value>()); +} + template @@ -1741,6 +1779,26 @@ distribute_local_to_global (const FullMatrix &local_matrix, } +// Complex case (maybe can be merged with the function above) +template +inline +void +ConstraintMatrix:: + distribute_local_to_global (const FullMatrix > &local_matrix, + const Vector > &local_vector, + const std::vector &local_dof_indices, + MatrixType &global_matrix, + VectorType &global_vector, + bool use_inhomogeneities_for_rhs) const +{ + // enter the internal function with the respective block information set, + // the actual implementation follows in the cm.templates.h file. + distribute_local_to_global (local_matrix, local_vector, local_dof_indices, + global_matrix, global_vector, use_inhomogeneities_for_rhs, + dealii::internal::bool2type::value>()); +} + + template inline diff --git a/deal.II/include/deal.II/lac/constraint_matrix.templates.h b/deal.II/include/deal.II/lac/constraint_matrix.templates.h index e8464dafc4..77ea578540 100644 --- a/deal.II/include/deal.II/lac/constraint_matrix.templates.h +++ b/deal.II/include/deal.II/lac/constraint_matrix.templates.h @@ -846,13 +846,13 @@ ConstraintMatrix::set_zero (VectorType &vec) const -template +template void ConstraintMatrix:: -distribute_local_to_global (const Vector &local_vector, +distribute_local_to_global (const Vector &local_vector, const std::vector &local_dof_indices, - VectorType &global_vector, - const FullMatrix &local_matrix) const + VectorType &global_vector, + const FullMatrix &local_matrix) const { Assert (sorted == true, ExcMatrixNotClosed()); AssertDimension (local_vector.size(), local_dof_indices.size()); @@ -862,10 +862,7 @@ distribute_local_to_global (const Vector &local_vector, const size_type n_local_dofs = local_vector.size(); if (lines.empty()) { - /* @whattodo Note: the Vector and FullMatrix coming in should be - complex */ - /* global_vector.add(local_dof_indices, local_vector); */ - Assert ((false), ExcMessage ("This function is corrupt: @whattodo")); + global_vector.add(local_dof_indices, local_vector); } else for (size_type i=0; i &local_vector, global_vector(local_dof_indices[j]) -= val * local_matrix(j,i); else { - const double matrix_entry = local_matrix(j,i); - if (matrix_entry == 0) + const LocalType matrix_entry = local_matrix(j,i); + if (matrix_entry == LocalType ()) // == 0 continue; const ConstraintLine &position_j = @@ -1855,15 +1852,15 @@ namespace internals // the origin of each global entry and finds out which data we need to // collect. static inline - double resolve_matrix_entry (const GlobalRowsFromLocal &global_rows, - const GlobalRowsFromLocal &global_cols, - const size_type i, - const size_type j, - const size_type loc_row, - const FullMatrix &local_matrix) + double resolve_matrix_entry (const GlobalRowsFromLocal &global_rows, + const GlobalRowsFromLocal &global_cols, + const size_type i, + const size_type j, + const size_type loc_row, + const FullMatrix &local_matrix) { const size_type loc_col = global_cols.local_row(j); - double col_val; + LocalType col_val; // case 1: row has direct contribution in local matrix. decide whether col // has a direct contribution. if not, set the value to zero. @@ -1886,8 +1883,8 @@ namespace internals // the direct and indirect references in the given column. for (size_type q=0; q + template inline void - resolve_matrix_row (const GlobalRowsFromLocal &global_rows, - const GlobalRowsFromLocal &global_cols, - const size_type i, - const size_type column_start, - const size_type column_end, - const FullMatrix &local_matrix, - size_type *&col_ptr, - number *&val_ptr) + resolve_matrix_row (const GlobalRowsFromLocal &global_rows, + const GlobalRowsFromLocal &global_cols, + const size_type i, + const size_type column_start, + const size_type column_end, + const FullMatrix &local_matrix, + size_type *&col_ptr, + number *&val_ptr) { if (column_end == column_start) return; @@ -1929,14 +1926,14 @@ namespace internals global_cols.have_indirect_rows() == false) { AssertIndexRange(loc_row, local_matrix.m()); - const double *matrix_ptr = &local_matrix(loc_row, 0); + const LocalType *matrix_ptr = &local_matrix(loc_row, 0); for (size_type j=column_start; j (col_val); *col_ptr++ = global_cols.global_row(j); @@ -1950,12 +1947,12 @@ namespace internals { for (size_type j=column_start; j (col_val); *col_ptr++ = global_cols.global_row(j); @@ -1970,14 +1967,14 @@ namespace internals // SparseMatrix. namespace dealiiSparseMatrix { - template + template static inline - void add_value (const double value, + void add_value (const LocalType value, const size_type row, const size_type column, SparseMatrixIterator &matrix_values) { - if (value != 0.) + if (value != LocalType ()) { while (matrix_values->column() < column) ++matrix_values; @@ -1992,15 +1989,15 @@ namespace internals // similar as before, now with shortcut for deal.II sparse matrices. this // lets us avoid using extra arrays, and does all the operations just in // place, i.e., in the respective matrix row - template + template inline void - resolve_matrix_row (const GlobalRowsFromLocal &global_rows, - const size_type i, - const size_type column_start, - const size_type column_end, - const FullMatrix &local_matrix, - SparseMatrix *sparse_matrix) + resolve_matrix_row (const GlobalRowsFromLocal &global_rows, + const size_type i, + const size_type column_start, + const size_type column_end, + const FullMatrix &local_matrix, + SparseMatrix *sparse_matrix) { if (column_end == column_start) return; @@ -2027,12 +2024,12 @@ namespace internals if (global_rows.have_indirect_rows() == false) { AssertIndexRange (loc_row, local_matrix.m()); - const double *matrix_ptr = &local_matrix(loc_row, 0); + const LocalType *matrix_ptr = &local_matrix(loc_row, 0); for (size_type j=column_start; jbegin(row)->value() += matrix_ptr[loc_row]; for (size_type j=column_start; jbegin(row)->value() += col_val; else @@ -2123,7 +2120,7 @@ namespace internals ++matrix_values; // jump over diagonal element for (size_type j=column_start; jbegin(row)->value() += col_val; @@ -2235,7 +2232,7 @@ add_this_index: inline void set_matrix_diagonals (const internals::GlobalRowsFromLocal &global_rows, const std::vector &local_dof_indices, - const FullMatrix &local_matrix, + const FullMatrix &local_matrix, const ConstraintMatrix &constraints, MatrixType &global_matrix, VectorType &global_vector, @@ -2243,18 +2240,18 @@ add_this_index: { if (global_rows.n_constraints() > 0) { - double average_diagonal = 0; + LocalType average_diagonal = 0; for (size_type i=0; i(local_matrix.m()); for (size_type i=0; i &local_vector, + const Vector &local_vector, const std::vector &local_dof_indices, - const FullMatrix &local_matrix) const + const FullMatrix &local_matrix) const { const size_type loc_row = global_rows.local_row(i); const size_type n_inhomogeneous_rows = global_rows.n_inhomogeneities(); - double val = 0; + LocalType val = 0; // has a direct contribution from some local entry. If we have inhomogeneous // constraints, compute the contribution of the inhomogeneity in the current // row. @@ -2471,7 +2464,7 @@ resolve_vector_entry (const size_type i, for (size_type q=0; q void ConstraintMatrix::distribute_local_to_global ( - const FullMatrix &local_matrix, - const Vector &local_vector, + const FullMatrix &local_matrix, + const Vector &local_vector, const std::vector &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, @@ -2576,12 +2569,12 @@ ConstraintMatrix::distribute_local_to_global ( // hand side. if (use_vectors == true) { - const double val = resolve_vector_entry (i, global_rows, - local_vector, - local_dof_indices, - local_matrix); + const LocalType val = resolve_vector_entry (i, global_rows, + local_vector, + local_dof_indices, + local_matrix); - if (val != 0) + if (val != LocalType ()) global_vector(row) += static_cast(val); } } @@ -2656,8 +2649,8 @@ ConstraintMatrix::distribute_local_to_global ( template void ConstraintMatrix:: -distribute_local_to_global (const FullMatrix &local_matrix, - const Vector &local_vector, +distribute_local_to_global (const FullMatrix &local_matrix, + const Vector &local_vector, const std::vector &local_dof_indices, MatrixType &global_matrix, VectorType &global_vector, @@ -2754,12 +2747,12 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (use_vectors == true) { - const double val = resolve_vector_entry (i, global_rows, - local_vector, - local_dof_indices, - local_matrix); + const LocalType val = resolve_vector_entry (i, global_rows, + local_vector, + local_dof_indices, + local_matrix); - if (val != 0) + if (val != LocalType ()) global_vector(global_indices[i]) += static_cast(val); } diff --git a/deal.II/source/lac/constraint_matrix.cc b/deal.II/source/lac/constraint_matrix.cc index 6c8cdbc796..bb5041cb0b 100644 --- a/deal.II/source/lac/constraint_matrix.cc +++ b/deal.II/source/lac/constraint_matrix.cc @@ -1732,6 +1732,27 @@ ConstraintMatrix::resolve_indices (std::vector &indices template void ConstraintMatrix::distribute(const VectorType &condensed,\ VectorType &uncondensed) const +#define COMPLEX_VECTOR_FUNCTIONS(VectorType) \ + template void ConstraintMatrix::condense(const VectorType &uncondensed, \ + VectorType &condensed) const; \ + template void ConstraintMatrix::condense(VectorType &vec) const; \ + template void ConstraintMatrix::condense(const SparseMatrix &uncondensed, \ + const VectorType &uncondensed_vector, \ + SparseMatrix &condensed, \ + VectorType &condensed_vector) const; \ + template void ConstraintMatrix::condense(const SparseMatrix &uncondensed, \ + const VectorType &uncondensed_vector, \ + SparseMatrix &condensed, \ + VectorType &condensed_vector) const; \ + template void ConstraintMatrix:: \ + distribute_local_to_global (const Vector > &, \ + const std::vector &, \ + VectorType &, \ + const FullMatrix > &) const; \ + template void ConstraintMatrix::distribute(const VectorType &condensed, \ + VectorType &uncondensed) const + + #define PARALLEL_VECTOR_FUNCTIONS(VectorType) \ template void ConstraintMatrix:: \ distribute_local_to_global (const Vector &, \ @@ -1739,10 +1760,12 @@ ConstraintMatrix::resolve_indices (std::vector &indices VectorType &, \ const FullMatrix &) const - #ifdef DEAL_II_WITH_PETSC VECTOR_FUNCTIONS(PETScWrappers::MPI::Vector); VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockVector); +// if PETSC_COMPLEX @whattodo check this works +VECTOR_FUNCTIONS_COMPLEX(PETScWrappers::MPI::Vector); +VECTOR_FUNCTIONS_COMPLEX(PETScWrappers::MPI::BlockVector); #endif #ifdef DEAL_II_WITH_TRILINOS @@ -1787,6 +1810,46 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector); bool , \ internal::bool2type) const +#define COMPLEX_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \ + template void ConstraintMatrix:: \ + distribute_local_to_global (const FullMatrix > &, \ + const Vector > &, \ + const std::vector &, \ + MatrixType &, \ + VectorType &, \ + bool , \ + internal::bool2type) const + +#define COMPLEX_MATRIX_FUNCTIONS(MatrixType) \ + template void ConstraintMatrix:: \ + distribute_local_to_global > > (const FullMatrix > &, \ + const Vector > &, \ + const std::vector &, \ + MatrixType &, \ + Vector > &, \ + bool , \ + internal::bool2type) const + +#define COMPLEX_BLOCK_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \ + template void ConstraintMatrix:: \ + distribute_local_to_global (const FullMatrix > &, \ + const Vector > &, \ + const std::vector &, \ + MatrixType &, \ + VectorType &, \ + bool , \ + internal::bool2type) const + +#define COMPLEX_BLOCK_MATRIX_FUNCTIONS(MatrixType) \ + template void ConstraintMatrix:: \ + distribute_local_to_global > > (const FullMatrix > &, \ + const Vector > &, \ + const std::vector &, \ + MatrixType &, \ + Vector > &, \ + bool , \ + internal::bool2type) const + MATRIX_FUNCTIONS(SparseMatrix); MATRIX_FUNCTIONS(SparseMatrix); MATRIX_FUNCTIONS(FullMatrix); @@ -1810,14 +1873,25 @@ MATRIX_VECTOR_FUNCTIONS(ChunkSparseMatrix, Vector); // BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ, Vector); #ifdef DEAL_II_WITH_PETSC -MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix); -BLOCK_MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix); -MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix); -BLOCK_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix); -MATRIX_VECTOR_FUNCTIONS(PETScWrappers::SparseMatrix, PETScWrappers::Vector); -BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::BlockSparseMatrix, PETScWrappers::BlockVector); -MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector); -BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix ,PETScWrappers::MPI::BlockVector); +// #if PETSC_COMPLEX @whattodo: Implement this and check +COMPLEX_MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix); +COMPLEX_BLOCK_MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix); +COMPLEX_MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix); +COMPLEX_BLOCK_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix); +COMPLEX_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::SparseMatrix, PETScWrappers::Vector); +COMPLEX_BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::BlockSparseMatrix, PETScWrappers::BlockVector); +COMPLEX_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector); +COMPLEX_BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix ,PETScWrappers::MPI::BlockVector); +// #else +// MATRIX_FUNCTIONS(PETScWrappers::SparseMatrix); +// BLOCK_MATRIX_FUNCTIONS(PETScWrappers::BlockSparseMatrix); +// MATRIX_FUNCTIONS(PETScWrappers::MPI::SparseMatrix); +// BLOCK_MATRIX_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix); +// MATRIX_VECTOR_FUNCTIONS(PETScWrappers::SparseMatrix, PETScWrappers::Vector); +// BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::BlockSparseMatrix, PETScWrappers::BlockVector); +// MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::SparseMatrix, PETScWrappers::MPI::Vector); +// BLOCK_MATRIX_VECTOR_FUNCTIONS(PETScWrappers::MPI::BlockSparseMatrix ,PETScWrappers::MPI::BlockVector); +// #endif #endif #ifdef DEAL_II_WITH_TRILINOS diff --git a/deal.II/source/lac/constraint_matrix.inst.in b/deal.II/source/lac/constraint_matrix.inst.in index 6320aaa188..be0efea90d 100644 --- a/deal.II/source/lac/constraint_matrix.inst.in +++ b/deal.II/source/lac/constraint_matrix.inst.in @@ -40,12 +40,21 @@ for (V: EXTERNAL_SEQUENTIAL_VECTORS) { template void ConstraintMatrix::condense(const V&, V&) const; template void ConstraintMatrix::condense(V&vec) const; - template void ConstraintMatrix::distribute_local_to_global ( - const Vector&, const std::vector &, V&, const FullMatrix&) const; template void ConstraintMatrix::distribute(const V&, V&) const; template void ConstraintMatrix::set_zero(V&) const; } +for (V: EXTERNAL_SEQUENTIAL_VECTORS) + { + template void ConstraintMatrix::distribute_local_to_global ( + const Vector&, const std::vector &, V&, const FullMatrix&) const; + } + +for (V: COMPLEX_EXTERNAL_SEQUENTIAL_VECTORS) + { + template void ConstraintMatrix::distribute_local_to_global ( + const Vector >&, const std::vector &, V&, const FullMatrix >&) const; + } for (V: EXTERNAL_PARALLEL_VECTORS) { -- 2.39.5