From: wolf Date: Thu, 18 May 2000 13:02:53 +0000 (+0000) Subject: Add some more functions for block sparse matrices. Not all yet tested. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=fbee65b9cae4a1b522427433421d7fb2e59a418d;p=dealii-svn.git Add some more functions for block sparse matrices. Not all yet tested. git-svn-id: https://svn.dealii.org/trunk@2889 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/dofs/dof_constraints.h b/deal.II/deal.II/include/dofs/dof_constraints.h index 1c6b88f5d1..0dd6b3c9b4 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -234,6 +234,14 @@ class ConstraintMatrix : public Subscriptor */ void condense (SparsityPattern &sparsity) const; + /** + * Same function as above, but + * condenses square block sparsity + * patterns. + */ + template + void condense (BlockSparsityPattern &sparsity) const; + /** * Condense a given matrix. The associated * matrix struct should be condensed and @@ -259,49 +267,94 @@ class ConstraintMatrix : public Subscriptor template void condense (SparseMatrix &matrix) const; + /** + * Same function as above, but + * condenses square block sparse + * matrices. + */ + template + void condense (BlockSparseMatrix &sparsity) const; + /** * Condense the given vector #uncondensed# * into #condensed#. It is the user's * responsibility to guarantee that all * entries of #condensed# be zero! + * + * The #VectorType# may be a + * #Vector#, + * #Vector#, + * #BlockVector<...>#, or any + * other type having the same + * interface. */ - template - void condense (const Vector &uncondensed, - Vector &condensed) const; + template + void condense (const VectorType &uncondensed, + VectorType &condensed) const; /** - * Condense the given vector in-place. + * Condense the given vector + * in-place. The #VectorType# may + * be a #Vector#, + * #Vector#, + * #BlockVector<...>#, or any + * other type having the same + * interface. */ - template - void condense (Vector &vec) const; + template + void condense (VectorType &vec) const; /** - * Re-distribute the elements of the vector - * #condensed# to #uncondensed#. It is the user's - * responsibility to guarantee that all - * entries of #uncondensed# be zero! + * Re-distribute the elements of + * the vector #condensed# to + * #uncondensed#. It is the + * user's responsibility to + * guarantee that all entries of + * #uncondensed# be zero! + * + * This function undoes the + * action of #condense# somehow, + * but it should be noted that it + * is not the inverse of + * #condense#. * - * This function undoes the action of - * #condense# somehow, but it should be noted - * that it is not the inverse of #condense#- + * The #VectorType# may be a + * #Vector#, + * #Vector#, + * #BlockVector<...>#, or any + * other type having the same + * interface. */ - template - void distribute (const Vector &condensed, - Vector &uncondensed) const; + template + void distribute (const VectorType &condensed, + VectorType &uncondensed) const; /** - * Re-distribute the elements of the vector - * in-place. + * Re-distribute the elements of + * the vector in-place. The + * #VectorType# may be a + * #Vector#, + * #Vector#, + * #BlockVector<...>#, or any + * other type having the same + * interface. */ - template - void distribute (Vector &vec) const; + template + void distribute (VectorType &vec) const; /** - * Delete hanging nodes in a vector. - * Sets all hanging node values to zero. + * Delete hanging nodes in a + * vector. Sets all hanging node + * values to zero. The + * #VectorType# may be a + * #Vector#, + * #Vector#, + * #BlockVector<...>#, or any + * other type having the same + * interface. */ - template - void set_zero (Vector &vec) const; + template + void set_zero (VectorType &vec) const; /** diff --git a/deal.II/deal.II/include/dofs/dof_constraints.templates.h b/deal.II/deal.II/include/dofs/dof_constraints.templates.h index 31cc4664dc..03c0de63d5 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.templates.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h @@ -15,7 +15,6 @@ #include -#include template @@ -261,10 +260,160 @@ ConstraintMatrix::condense (SparseMatrix &uncondensed) const -template +template void -ConstraintMatrix::condense (const Vector &uncondensed, - Vector &condensed) const +ConstraintMatrix::condense (BlockSparseMatrix &uncondensed) const +{ + const BlockSparsityPattern & + sparsity = uncondensed.get_sparsity_pattern (); + + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed()); + Assert (sparsity.n_rows() == sparsity.n_cols(), + ExcMatrixNotSquare()); + Assert (sparsity.get_column_indices() == sparsity.get_row_indices(), + ExcMatrixNotSquare()); + + const BlockIndices & + index_mapping = sparsity.get_column_indices(); + + // store for each index whether it + // must be distributed or not. If entry + // is -1, no distribution is necessary. + // otherwise, the number states which + // line in the constraint matrix handles + // this index + vector distribute (sparsity.n_rows(), -1); + + for (unsigned int c=0; c(c); + + const unsigned int n_rows = sparsity.n_rows(); + for (unsigned int row=0; row + block_index = index_mapping.global_to_local(row); + + // regular line. loop over + // all columns and see + // whether this column must + // be distributed + for (unsigned int block_col=0; block_col +void +ConstraintMatrix::condense (const VectorType &uncondensed, + VectorType &condensed) const { Assert (sorted == true, ExcMatrixNotClosed()); Assert (condensed.size()+n_constraints() == uncondensed.size(), @@ -335,9 +484,9 @@ ConstraintMatrix::condense (const Vector &uncondensed, -template +template void -ConstraintMatrix::condense (Vector &vec) const +ConstraintMatrix::condense (VectorType &vec) const { Assert (sorted == true, ExcMatrixNotClosed()); @@ -366,9 +515,9 @@ ConstraintMatrix::condense (Vector &vec) const -template +template void -ConstraintMatrix::set_zero (Vector &vec) const +ConstraintMatrix::set_zero (VectorType &vec) const { Assert (sorted == true, ExcMatrixNotClosed()); @@ -392,10 +541,10 @@ ConstraintMatrix::set_zero (Vector &vec) const -template +template void -ConstraintMatrix::distribute (const Vector &condensed, - Vector &uncondensed) const +ConstraintMatrix::distribute (const VectorType &condensed, + VectorType &uncondensed) const { Assert (sorted == true, ExcMatrixNotClosed()); Assert (condensed.size()+n_constraints() == uncondensed.size(), @@ -465,9 +614,9 @@ ConstraintMatrix::distribute (const Vector &condensed, -template +template void -ConstraintMatrix::distribute (Vector &vec) const +ConstraintMatrix::distribute (VectorType &vec) const { Assert (sorted == true, ExcMatrixNotClosed()); diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index 1df1b5b143..c6d6aa4ff5 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -13,8 +13,13 @@ #include +#include + #include #include +#include +#include +#include #include #include #include @@ -28,11 +33,15 @@ ConstraintMatrix::ConstraintLine::operator < (const ConstraintLine &a) const }; + ConstraintMatrix::ConstraintMatrix () : - lines(), sorted(false) {}; + lines(), + sorted(false) +{}; -void ConstraintMatrix::add_line (const unsigned int line) { +void ConstraintMatrix::add_line (const unsigned int line) +{ Assert (sorted==false, ExcMatrixIsClosed()); // check whether line already exists; @@ -48,9 +57,11 @@ void ConstraintMatrix::add_line (const unsigned int line) { }; + void ConstraintMatrix::add_entry (const unsigned int line, const unsigned int column, - const double value) { + const double value) +{ Assert (sorted==false, ExcMatrixIsClosed()); vector::iterator line_ptr; @@ -89,7 +100,9 @@ void ConstraintMatrix::add_entry (const unsigned int line, }; -void ConstraintMatrix::close () { + +void ConstraintMatrix::close () +{ Assert (sorted==false, ExcMatrixIsClosed()); // sort the entries in the different lines @@ -148,14 +161,18 @@ void ConstraintMatrix::close () { }; -void ConstraintMatrix::clear () { + +void ConstraintMatrix::clear () +{ lines = vector(); sorted = false; }; + void ConstraintMatrix::condense (const SparsityPattern &uncondensed, - SparsityPattern &condensed) const { + SparsityPattern &condensed) const +{ Assert (sorted == true, ExcMatrixNotClosed()); Assert (uncondensed.is_compressed() == true, ExcMatrixNotClosed()); Assert (uncondensed.n_rows() == uncondensed.n_cols(), @@ -262,7 +279,9 @@ void ConstraintMatrix::condense (const SparsityPattern &uncondensed, }; -void ConstraintMatrix::condense (SparsityPattern &sparsity) const { + +void ConstraintMatrix::condense (SparsityPattern &sparsity) const +{ Assert (sorted == true, ExcMatrixNotClosed()); Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed()); Assert (sparsity.n_rows() == sparsity.n_cols(), @@ -276,6 +295,91 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const { // this index vector distribute(sparsity.n_rows(), -1); + for (unsigned int c=0; c(c); + + unsigned int n_rows = sparsity.n_rows(); + for (unsigned int row=0; row +void ConstraintMatrix::condense (BlockSparsityPattern &sparsity) const +{ + Assert (false, ExcInternalError()); +/* + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (sparsity.is_compressed() == false, ExcMatrixIsClosed()); + Assert (sparsity.n_rows() == sparsity.n_cols(), + ExcMatrixNotSquare()); + Assert (sparsity.get_column_indices() == sparsity.get_row_indices(), + ExcMatrixNotSquare()); + + // store for each index whether it + // must be distributed or not. If entry + // is -1, no distribution is necessary. + // otherwise, the number states which + // line in the constraint matrix handles + // this index + vector distribute(sparsity.n_rows(), -1); + for (unsigned int c=0; c(c); @@ -337,15 +441,18 @@ void ConstraintMatrix::condense (SparsityPattern &sparsity) const { }; sparsity.compress(); +*/ }; + unsigned int ConstraintMatrix::n_constraints () const { return lines.size(); }; + bool ConstraintMatrix::is_constrained (const unsigned int index) const { if (sorted == true) @@ -397,32 +504,82 @@ void ConstraintMatrix::print (ostream &out) const { AssertThrow (out, ExcIO()); }; -#include -#define number double -template void ConstraintMatrix::condense(const SparseMatrix &uncondensed, - SparseMatrix &condensed) const; -template void ConstraintMatrix::condense(SparseMatrix &uncondensed) const; -template void ConstraintMatrix::condense(const Vector &uncondensed, - Vector &condensed) const; -template void ConstraintMatrix::condense(Vector &vec) const; -template void ConstraintMatrix::set_zero(Vector &vec) const; -template void ConstraintMatrix::distribute(const Vector &condensed, - Vector &uncondensed) const; -template void ConstraintMatrix::distribute(Vector &vec) const; - -#undef number -#define number float - -template void ConstraintMatrix::condense(const SparseMatrix &uncondensed, - SparseMatrix &condensed) const; -template void ConstraintMatrix::condense(SparseMatrix &uncondensed) const; -template void ConstraintMatrix::condense(const Vector &uncondensed, - Vector &condensed) const; -template void ConstraintMatrix::condense(Vector &vec) const; -template void ConstraintMatrix::set_zero(Vector &vec) const; -template void ConstraintMatrix::distribute(const Vector &condensed, - Vector &uncondensed) const; -template void ConstraintMatrix::distribute(Vector &vec) const; - -#undef number + + + +// explicit instantiations +// +// define a list of functions for vectors and matrices, respectively, +// where the vector/matrix can be replaced using a preprocessor +// variable VectorType/MatrixType. note that we cannot do so by using +// a preprocessor function with one arg, since +// #vector_functions(BlockVector<2,double>)# is not recognized as one +// arg, and putting parentheses around the arg yields incorrect +// syntax... + +#define vector_functions \ + template void ConstraintMatrix::condense(const VectorType &uncondensed,\ + VectorType &condensed) const;\ + template void ConstraintMatrix::condense(VectorType &vec) const;\ + template void ConstraintMatrix::set_zero(VectorType &vec) const;\ + template void ConstraintMatrix::distribute(const VectorType &condensed,\ + VectorType &uncondensed) const;\ + template void ConstraintMatrix::distribute(VectorType &vec) const; + + + +#define matrix_functions_1 \ + template void ConstraintMatrix::condense(const MatrixType &uncondensed,\ + MatrixType &condensed) const; +#define matrix_functions_2 \ + template void ConstraintMatrix::condense(MatrixType &uncondensed) const; + + + +#define VectorType Vector +vector_functions; +#undef VectorType + +#define VectorType Vector +vector_functions; +#undef VectorType + +#define VectorType BlockVector<2,double> +vector_functions; +#undef VectorType + +#define VectorType BlockVector<3,double> +vector_functions; +#undef VectorType + + + +#define MatrixType SparseMatrix +matrix_functions_1; +matrix_functions_2; +#undef MatrixType + +#define MatrixType SparseMatrix +matrix_functions_1; +matrix_functions_2; +#undef MatrixType + +// block sparse matrices are only implemented for one of the two matrix functions +#define MatrixType BlockSparseMatrix +matrix_functions_2; +#undef MatrixType + +#define MatrixType BlockSparseMatrix +matrix_functions_2; +#undef MatrixType + + +#define MatrixType BlockSparsityPattern<2,2> +matrix_functions_2; +#undef MatrixType + +#define MatrixType BlockSparsityPattern<3,3> +matrix_functions_2; +#undef MatrixType +