From: Denis Davydov Date: Tue, 16 Jun 2015 21:18:22 +0000 (+0200) Subject: extend ConstraintMatrix to cover complex-valued cases X-Git-Tag: v8.3.0-rc1~104^2~6 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7eedf05738954fe94af74be84859501a6d3e97fe;p=dealii.git extend ConstraintMatrix to cover complex-valued cases --- diff --git a/include/deal.II/lac/constraint_matrix.h b/include/deal.II/lac/constraint_matrix.h index c42b67265b..d2549ae3b9 100644 --- a/include/deal.II/lac/constraint_matrix.h +++ b/include/deal.II/lac/constraint_matrix.h @@ -717,12 +717,12 @@ public: * site. There is no locking mechanism inside this method to prevent data * races. */ - template + template void - 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; + const FullMatrix &local_matrix) const; /** * Enter a single value into a result vector, obeying constraints. @@ -818,7 +818,7 @@ public: */ template void - distribute_local_to_global (const FullMatrix &local_matrix, + distribute_local_to_global (const FullMatrix &local_matrix, const std::vector &local_dof_indices, MatrixType &global_matrix) const; @@ -851,7 +851,7 @@ public: */ template void - distribute_local_to_global (const FullMatrix &local_matrix, + distribute_local_to_global (const FullMatrix &local_matrix, const std::vector &row_indices, const std::vector &col_indices, MatrixType &global_matrix) const; @@ -874,8 +874,8 @@ public: */ 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, @@ -1239,8 +1239,8 @@ private: */ 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, @@ -1253,8 +1253,8 @@ private: */ 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, @@ -1310,12 +1310,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; /** * The assignment operator is not implemented for performance reasons. You @@ -1639,13 +1640,13 @@ template inline void ConstraintMatrix:: -distribute_local_to_global (const FullMatrix &local_matrix, +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); + Vector dummy(0); distribute_local_to_global (local_matrix, dummy, local_dof_indices, global_matrix, dummy, false, dealii::internal::bool2type::value>()); @@ -1653,12 +1654,13 @@ distribute_local_to_global (const FullMatrix &local_matrix, + template inline 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, @@ -1673,6 +1675,7 @@ distribute_local_to_global (const FullMatrix &local_matrix, + template inline void diff --git a/include/deal.II/lac/constraint_matrix.templates.h b/include/deal.II/lac/constraint_matrix.templates.h index dad4ca4144..c0194a7cf4 100644 --- a/include/deal.II/lac/constraint_matrix.templates.h +++ b/include/deal.II/lac/constraint_matrix.templates.h @@ -601,13 +601,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()); @@ -649,9 +649,9 @@ distribute_local_to_global (const Vector &local_vector, continue; } - const double matrix_entry = local_matrix(j,i); + const LocalType matrix_entry = local_matrix(j,i); - if (matrix_entry == 0) + if (matrix_entry == LocalType()) continue; const ConstraintLine &position_j = @@ -1535,16 +1535,17 @@ namespace internals // resolves constraints of one column at the innermost loop. goes through // the origin of each global entry and finds out which data we need to // collect. + template 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) + LocalType 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. @@ -1567,8 +1568,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, @@ -1593,7 +1594,7 @@ namespace internals const size_type i, const size_type column_start, const size_type column_end, - const FullMatrix &local_matrix, + const FullMatrix &local_matrix, size_type *&col_ptr, number *&val_ptr) { @@ -1610,14 +1611,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); @@ -1631,12 +1632,12 @@ namespace internals { for (size_type j=column_start; j (col_val); *col_ptr++ = global_cols.global_row(j); @@ -1651,15 +1652,15 @@ 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) { (void)row; - if (value != 0.) + if (value != LocalType ()) { while (matrix_values->column() < column) ++matrix_values; @@ -1674,14 +1675,14 @@ 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, + const FullMatrix &local_matrix, SparseMatrix *sparse_matrix) { if (column_end == column_start) @@ -1709,12 +1710,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 @@ -1805,8 +1806,8 @@ namespace internals ++matrix_values; // jump over diagonal element for (size_type j=column_start; jbegin(row)->value() += col_val; else @@ -1917,7 +1918,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, @@ -1925,9 +1926,9 @@ add_this_index: { if (global_rows.n_constraints() > 0) { - double average_diagonal = 0; + typename MatrixType::value_type average_diagonal = typename MatrixType::value_type(); for (size_type i=0; i(local_matrix.m()); for (size_type i=0; i &local_dof_indices, // Resolve the constraints from the vector and apply inhomogeneities. +template< typename LocalType> inline -double +LocalType ConstraintMatrix:: 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 { 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. @@ -2150,7 +2152,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, @@ -2255,12 +2257,12 @@ ConstraintMatrix::distribute_local_to_global ( // hand side. if (use_vectors == true) { - const double val = resolve_vector_entry (i, global_rows, + const number val = resolve_vector_entry (i, global_rows, local_vector, local_dof_indices, local_matrix); - if (val != 0) + if (val != number ()) global_vector(row) += static_cast(val); } } @@ -2275,12 +2277,12 @@ ConstraintMatrix::distribute_local_to_global ( template void ConstraintMatrix::distribute_local_to_global ( - const FullMatrix &local_matrix, + const FullMatrix &local_matrix, const std::vector &row_indices, const std::vector &col_indices, MatrixType &global_matrix) const { - typedef double number; + typedef typename MatrixType::value_type number; AssertDimension (local_matrix.m(), row_indices.size()); AssertDimension (local_matrix.n(), col_indices.size()); @@ -2331,8 +2333,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, @@ -2429,12 +2431,12 @@ distribute_local_to_global (const FullMatrix &local_matrix, if (use_vectors == true) { - const double val = resolve_vector_entry (i, global_rows, + const number val = resolve_vector_entry (i, global_rows, local_vector, local_dof_indices, local_matrix); - if (val != 0) + if (val != number ()) global_vector(global_indices[i]) += static_cast(val); } diff --git a/source/lac/constraint_matrix.cc b/source/lac/constraint_matrix.cc index 479d9f1894..eb8d83675b 100644 --- a/source/lac/constraint_matrix.cc +++ b/source/lac/constraint_matrix.cc @@ -1247,8 +1247,8 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector); #define MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \ template void ConstraintMatrix:: \ - distribute_local_to_global (const FullMatrix &, \ - const Vector &, \ + distribute_local_to_global (const FullMatrix &, \ + const Vector &, \ const std::vector &, \ MatrixType &, \ VectorType &, \ @@ -1256,17 +1256,17 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector); internal::bool2type) const #define MATRIX_FUNCTIONS(MatrixType) \ template void ConstraintMatrix:: \ - distribute_local_to_global > (const FullMatrix &, \ - const Vector &, \ - const std::vector &, \ - MatrixType &, \ - Vector &, \ - bool , \ - internal::bool2type) const + distribute_local_to_global > (const FullMatrix &, \ + const Vector &, \ + const std::vector &, \ + MatrixType &, \ + Vector &, \ + bool , \ + internal::bool2type) const #define BLOCK_MATRIX_VECTOR_FUNCTIONS(MatrixType, VectorType) \ template void ConstraintMatrix:: \ - distribute_local_to_global (const FullMatrix &, \ - const Vector &, \ + distribute_local_to_global (const FullMatrix &, \ + const Vector &, \ const std::vector &, \ MatrixType &, \ VectorType &, \ @@ -1274,32 +1274,29 @@ PARALLEL_VECTOR_FUNCTIONS(TrilinosWrappers::MPI::BlockVector); internal::bool2type) const #define BLOCK_MATRIX_FUNCTIONS(MatrixType) \ template void ConstraintMatrix:: \ - distribute_local_to_global > (const FullMatrix &, \ - const Vector &, \ - const std::vector &, \ - MatrixType &, \ - Vector &, \ - bool , \ - internal::bool2type) const + 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); MATRIX_FUNCTIONS(FullMatrix); -MATRIX_VECTOR_FUNCTIONS(SparseMatrix, Vector); +MATRIX_FUNCTIONS(FullMatrix >); BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix); BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrix); BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix, BlockVector); BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix, BlockVector); -BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrix, BlockVector); MATRIX_FUNCTIONS(SparseMatrixEZ); MATRIX_FUNCTIONS(SparseMatrixEZ); MATRIX_FUNCTIONS(ChunkSparseMatrix); MATRIX_FUNCTIONS(ChunkSparseMatrix); -MATRIX_VECTOR_FUNCTIONS(SparseMatrixEZ, Vector); -MATRIX_VECTOR_FUNCTIONS(ChunkSparseMatrix, Vector); // BLOCK_MATRIX_FUNCTIONS(BlockSparseMatrixEZ); // BLOCK_MATRIX_VECTOR_FUNCTIONS(BlockSparseMatrixEZ, Vector); @@ -1365,7 +1362,7 @@ BLOCK_SPARSITY_FUNCTIONS(TrilinosWrappers::BlockSparsityPattern); #define ONLY_MATRIX_FUNCTIONS(MatrixType) \ template void ConstraintMatrix::distribute_local_to_global (\ - const FullMatrix &, \ + const FullMatrix &, \ const std::vector &, \ const std::vector &, \ MatrixType &) const