From: guido Date: Mon, 26 Apr 1999 20:30:17 +0000 (+0000) Subject: number templates in ConstraintMatrix X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d8414bfd4127bf423ee94982dd1c88401ce42749;p=dealii-svn.git number templates in ConstraintMatrix git-svn-id: https://svn.dealii.org/trunk@1211 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 1a86af56d6..92a14f7bc8 100644 --- a/deal.II/deal.II/include/dofs/dof_constraints.h +++ b/deal.II/deal.II/include/dofs/dof_constraints.h @@ -196,8 +196,9 @@ class ConstraintMatrix { * The constraint matrix object must be * closed to call this function. */ - void condense (const SparseMatrix &uncondensed, - SparseMatrix &condensed) const; + template + void condense (const SparseMatrix &uncondensed, + SparseMatrix &condensed) const; /** * This function does much the same as @@ -206,7 +207,8 @@ class ConstraintMatrix { * documentation of this class for more * detailed information. */ - void condense (SparseMatrix &matrix) const; + template + void condense (SparseMatrix &matrix) const; /** * Condense the given vector #uncondensed# @@ -214,13 +216,15 @@ class ConstraintMatrix { * responsibility to guarantee that all * entries of #condensed# be zero! */ - void condense (const Vector &uncondensed, - Vector &condensed) const; + template + void condense (const Vector &uncondensed, + Vector &condensed) const; /** * Condense the given vector in-place. */ - void condense (Vector &vec) const; + template + void condense (Vector &vec) const; /** * Re-distribute the elements of the vector @@ -232,14 +236,16 @@ class ConstraintMatrix { * #condense# somehow, but it should be noted * that it is not the inverse of #condense#- */ - void distribute (const Vector &condensed, - Vector &uncondensed) const; + template + void distribute (const Vector &condensed, + Vector &uncondensed) const; /** * Re-distribute the elements of the vector * in-place. */ - void distribute (Vector &vec) const; + template + void distribute (Vector &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 new file mode 100644 index 0000000000..4e964898b8 --- /dev/null +++ b/deal.II/deal.II/include/dofs/dof_constraints.templates.h @@ -0,0 +1,442 @@ +template +void +ConstraintMatrix::condense (const SparseMatrix &uncondensed, + SparseMatrix &condensed) const +{ + const SparseMatrixStruct &uncondensed_struct = uncondensed.get_sparsity_pattern (); + + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed()); + Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed()); + Assert (uncondensed_struct.n_rows() == uncondensed_struct.n_cols(), + ExcMatrixNotSquare()); + Assert (condensed.n() == condensed.m(), + ExcMatrixNotSquare()); + Assert (condensed.n()+n_constraints() == uncondensed.n(), + ExcWrongDimension()); + + // store for each line of the matrix + // its new line number + // after compression. If the shift is + // -1, this line will be condensed away + vector new_line; + + new_line.reserve (uncondensed_struct.n_rows()); + + vector::const_iterator next_constraint = lines.begin(); + unsigned int shift = 0; + unsigned int n_rows = uncondensed_struct.n_rows(); + + if (next_constraint == lines.end()) + // if no constraint is to be handled + for (unsigned int row=0; row!=n_rows; ++row) + new_line.push_back (row); + else + for (unsigned int row=0; row!=n_rows; ++row) + if (row == next_constraint->line) + { + // this line is constrained + new_line.push_back (-1); + // note that #lines# is ordered + ++shift; + ++next_constraint; + if (next_constraint == lines.end()) + // nothing more to do; finish rest + // of loop + { + for (unsigned int i=row+1; i::const_iterator c = lines.begin(); + while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) + ++c; + + for (unsigned int q=0; q!=c->entries.size(); ++q) + // distribute to rows with + // appropriate weight + condensed.add (new_line[row], new_line[c->entries[q].first], + uncondensed.global_entry(j) * c->entries[q].second); + } + else + // line must be distributed + { + for (unsigned int j=uncondensed_struct.get_rowstart_indices()[row]; + jentries.size(); ++q) + condensed.add (new_line[next_constraint->entries[q].first], + new_line[uncondensed_struct.get_column_numbers()[j]], + uncondensed.global_entry(j) * + next_constraint->entries[q].second); + + else + // not only this line but + // also this col is constrained + { + // let c point to the constraint + // of this column + vector::const_iterator c = lines.begin(); + while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) + ++c; + + for (unsigned int p=0; p!=c->entries.size(); ++p) + for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) + condensed.add (new_line[next_constraint->entries[q].first], + new_line[c->entries[p].first], + uncondensed.global_entry(j) * + next_constraint->entries[q].second * + c->entries[p].second); + }; + + ++next_constraint; + }; +}; + + + +template +void +ConstraintMatrix::condense (SparseMatrix &uncondensed) const +{ + const SparseMatrixStruct &sparsity = uncondensed.get_sparsity_pattern (); + + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed()); + Assert (sparsity.n_rows() == sparsity.n_cols(), + 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); + + int n_rows = sparsity.n_rows(); + for (int row=0; row +void +ConstraintMatrix::condense (const Vector &uncondensed, + Vector &condensed) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (condensed.size()+n_constraints() == uncondensed.size(), + ExcWrongDimension()); + + // store for each line of the vector + // its new line number + // after compression. If the shift is + // -1, this line will be condensed away + vector new_line; + + new_line.reserve (uncondensed.size()); + + vector::const_iterator next_constraint = lines.begin(); + unsigned int shift = 0; + unsigned int n_rows = uncondensed.size(); + + if (next_constraint == lines.end()) + // if no constraint is to be handled + for (unsigned int row=0; row!=n_rows; ++row) + new_line.push_back (row); + else + for (unsigned int row=0; row!=n_rows; ++row) + if (row == next_constraint->line) + { + // this line is constrained + new_line.push_back (-1); + // note that #lines# is ordered + ++shift; + ++next_constraint; + if (next_constraint == lines.end()) + // nothing more to do; finish rest + // of loop + { + for (unsigned int i=row+1; ientries.size(); ++q) + condensed(new_line[next_constraint->entries[q].first]) + += + uncondensed(row) * next_constraint->entries[q].second; + + ++next_constraint; + }; +}; + + +template +void +ConstraintMatrix::condense (Vector &vec) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + + if (lines.size() == 0) + // nothing to do + return; + + vector::const_iterator next_constraint = lines.begin(); + for (unsigned int row=0; rowline) + // line must be distributed + { + for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) + vec(next_constraint->entries[q].first) + += + vec(row) * next_constraint->entries[q].second; + // set entry to zero + vec(row) = 0.; + + ++next_constraint; + if (next_constraint == lines.end()) + // nothing more to do + break; + }; +}; + + + +template +void +ConstraintMatrix::distribute (const Vector &condensed, + Vector &uncondensed) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + Assert (condensed.size()+n_constraints() == uncondensed.size(), + ExcWrongDimension()); + + // store for each line of the new vector + // its old line number before + // distribution. If the shift is + // -1, this line was condensed away + vector old_line; + + old_line.reserve (uncondensed.size()); + + vector::const_iterator next_constraint = lines.begin(); + unsigned int shift = 0; + unsigned int n_rows = uncondensed.size(); + + if (next_constraint == lines.end()) + // if no constraint is to be handled + for (unsigned int row=0; row!=n_rows; ++row) + old_line.push_back (row); + else + for (unsigned int row=0; row!=n_rows; ++row) + if (row == next_constraint->line) + { + // this line is constrained + old_line.push_back (-1); + // note that #lines# is ordered + ++shift; + ++next_constraint; + if (next_constraint == lines.end()) + // nothing more to do; finish rest + // of loop + { + for (unsigned int i=row+1; ientries.size(); ++i) + uncondensed(line) += (condensed(old_line[next_constraint->entries[i].first]) * + next_constraint->entries[i].second); + ++next_constraint; + }; +}; + + +template +void +ConstraintMatrix::distribute (Vector &vec) const +{ + Assert (sorted == true, ExcMatrixNotClosed()); + + vector::const_iterator next_constraint = lines.begin(); + for (; next_constraint != lines.end(); ++next_constraint) + { + // make entry in line next_constraint.line + // first set it to zero + vec(next_constraint->line) = 0.; + // then add the different contributions + for (unsigned int i=0; ientries.size(); ++i) + vec(next_constraint->line) += (vec(next_constraint->entries[i].first) * + next_constraint->entries[i].second); + }; +}; + + diff --git a/deal.II/deal.II/source/dofs/dof_constraints.cc b/deal.II/deal.II/source/dofs/dof_constraints.cc index cbbeed1db4..1af3e7944b 100644 --- a/deal.II/deal.II/source/dofs/dof_constraints.cc +++ b/deal.II/deal.II/source/dofs/dof_constraints.cc @@ -307,435 +307,6 @@ void ConstraintMatrix::condense (SparseMatrixStruct &sparsity) const { -void ConstraintMatrix::condense (const SparseMatrix &uncondensed, - SparseMatrix &condensed) const { - const SparseMatrixStruct &uncondensed_struct = uncondensed.get_sparsity_pattern (); - - Assert (sorted == true, ExcMatrixNotClosed()); - Assert (uncondensed_struct.is_compressed() == true, ExcMatrixNotClosed()); - Assert (condensed.get_sparsity_pattern().is_compressed() == true, ExcMatrixNotClosed()); - Assert (uncondensed_struct.n_rows() == uncondensed_struct.n_cols(), - ExcMatrixNotSquare()); - Assert (condensed.n() == condensed.m(), - ExcMatrixNotSquare()); - Assert (condensed.n()+n_constraints() == uncondensed.n(), - ExcWrongDimension()); - - // store for each line of the matrix - // its new line number - // after compression. If the shift is - // -1, this line will be condensed away - vector new_line; - - new_line.reserve (uncondensed_struct.n_rows()); - - vector::const_iterator next_constraint = lines.begin(); - unsigned int shift = 0; - unsigned int n_rows = uncondensed_struct.n_rows(); - - if (next_constraint == lines.end()) - // if no constraint is to be handled - for (unsigned int row=0; row!=n_rows; ++row) - new_line.push_back (row); - else - for (unsigned int row=0; row!=n_rows; ++row) - if (row == next_constraint->line) - { - // this line is constrained - new_line.push_back (-1); - // note that #lines# is ordered - ++shift; - ++next_constraint; - if (next_constraint == lines.end()) - // nothing more to do; finish rest - // of loop - { - for (unsigned int i=row+1; i::const_iterator c = lines.begin(); - while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) - ++c; - - for (unsigned int q=0; q!=c->entries.size(); ++q) - // distribute to rows with - // appropriate weight - condensed.add (new_line[row], new_line[c->entries[q].first], - uncondensed.global_entry(j) * c->entries[q].second); - } - else - // line must be distributed - { - for (unsigned int j=uncondensed_struct.get_rowstart_indices()[row]; - jentries.size(); ++q) - condensed.add (new_line[next_constraint->entries[q].first], - new_line[uncondensed_struct.get_column_numbers()[j]], - uncondensed.global_entry(j) * - next_constraint->entries[q].second); - - else - // not only this line but - // also this col is constrained - { - // let c point to the constraint - // of this column - vector::const_iterator c = lines.begin(); - while (c->line != (unsigned int)uncondensed_struct.get_column_numbers()[j]) - ++c; - - for (unsigned int p=0; p!=c->entries.size(); ++p) - for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) - condensed.add (new_line[next_constraint->entries[q].first], - new_line[c->entries[p].first], - uncondensed.global_entry(j) * - next_constraint->entries[q].second * - c->entries[p].second); - }; - - ++next_constraint; - }; -}; - - - -void ConstraintMatrix::condense (SparseMatrix &uncondensed) const { - const SparseMatrixStruct &sparsity = uncondensed.get_sparsity_pattern (); - - Assert (sorted == true, ExcMatrixNotClosed()); - Assert (sparsity.is_compressed() == true, ExcMatrixNotClosed()); - Assert (sparsity.n_rows() == sparsity.n_cols(), - 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); - - int n_rows = sparsity.n_rows(); - for (int row=0; row &uncondensed, - Vector &condensed) const { - Assert (sorted == true, ExcMatrixNotClosed()); - Assert (condensed.size()+n_constraints() == uncondensed.size(), - ExcWrongDimension()); - - // store for each line of the vector - // its new line number - // after compression. If the shift is - // -1, this line will be condensed away - vector new_line; - - new_line.reserve (uncondensed.size()); - - vector::const_iterator next_constraint = lines.begin(); - unsigned int shift = 0; - unsigned int n_rows = uncondensed.size(); - - if (next_constraint == lines.end()) - // if no constraint is to be handled - for (unsigned int row=0; row!=n_rows; ++row) - new_line.push_back (row); - else - for (unsigned int row=0; row!=n_rows; ++row) - if (row == next_constraint->line) - { - // this line is constrained - new_line.push_back (-1); - // note that #lines# is ordered - ++shift; - ++next_constraint; - if (next_constraint == lines.end()) - // nothing more to do; finish rest - // of loop - { - for (unsigned int i=row+1; ientries.size(); ++q) - condensed(new_line[next_constraint->entries[q].first]) - += - uncondensed(row) * next_constraint->entries[q].second; - - ++next_constraint; - }; -}; - - - -void ConstraintMatrix::condense (Vector &vec) const { - Assert (sorted == true, ExcMatrixNotClosed()); - - if (lines.size() == 0) - // nothing to do - return; - - vector::const_iterator next_constraint = lines.begin(); - for (unsigned int row=0; rowline) - // line must be distributed - { - for (unsigned int q=0; q!=next_constraint->entries.size(); ++q) - vec(next_constraint->entries[q].first) - += - vec(row) * next_constraint->entries[q].second; - // set entry to zero - vec(row) = 0.; - - ++next_constraint; - if (next_constraint == lines.end()) - // nothing more to do - break; - }; -}; - - - - -void ConstraintMatrix::distribute (const Vector &condensed, - Vector &uncondensed) const { - Assert (sorted == true, ExcMatrixNotClosed()); - Assert (condensed.size()+n_constraints() == uncondensed.size(), - ExcWrongDimension()); - - // store for each line of the new vector - // its old line number before - // distribution. If the shift is - // -1, this line was condensed away - vector old_line; - - old_line.reserve (uncondensed.size()); - - vector::const_iterator next_constraint = lines.begin(); - unsigned int shift = 0; - unsigned int n_rows = uncondensed.size(); - - if (next_constraint == lines.end()) - // if no constraint is to be handled - for (unsigned int row=0; row!=n_rows; ++row) - old_line.push_back (row); - else - for (unsigned int row=0; row!=n_rows; ++row) - if (row == next_constraint->line) - { - // this line is constrained - old_line.push_back (-1); - // note that #lines# is ordered - ++shift; - ++next_constraint; - if (next_constraint == lines.end()) - // nothing more to do; finish rest - // of loop - { - for (unsigned int i=row+1; ientries.size(); ++i) - uncondensed(line) += (condensed(old_line[next_constraint->entries[i].first]) * - next_constraint->entries[i].second); - ++next_constraint; - }; -}; - - - -void ConstraintMatrix::distribute (Vector &vec) const { - Assert (sorted == true, ExcMatrixNotClosed()); - - vector::const_iterator next_constraint = lines.begin(); - for (; next_constraint != lines.end(); ++next_constraint) - { - // make entry in line next_constraint.line - // first set it to zero - vec(next_constraint->line) = 0.; - // then add the different contributions - for (unsigned int i=0; ientries.size(); ++i) - vec(next_constraint->line) += (vec(next_constraint->entries[i].first) * - next_constraint->entries[i].second); - }; -}; - - - - unsigned int ConstraintMatrix::n_constraints () const { return lines.size(); }; @@ -751,3 +322,29 @@ 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::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::distribute(const Vector &condensed, + Vector &uncondensed) const; +template void ConstraintMatrix::distribute(Vector &vec) const;