]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Allow some block matrix-usual vector operations.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Sep 2003 16:09:15 +0000 (16:09 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Tue, 9 Sep 2003 16:09:15 +0000 (16:09 +0000)
git-svn-id: https://svn.dealii.org/trunk@7968 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/c-4-0.html
tests/lac/sparse_matrices.cc
tests/lac/testmatrix.cc

index 8e3c3cba99f8ca263f824093890086141d5a4cd7..6686c6773e57f62494f7948725702e36d1d373f4 100644 (file)
@@ -142,6 +142,16 @@ contributor's names are abbreviated by WB (Wolfgang Bangerth), GK
 <h3>lac</h3>
 
 <ol>
+  <li> <p> New: Some of the member matrix-vector functions of the
+       <code class="class">BlockSparseMatrix</code> class that
+       previously could only be used with arguments of type <code
+       class="class">BlockVector</code> can now also be used with the
+       usual <code class="class">Vector</code> class provided the
+       block matrix has only one block row or column.
+       <br>
+       (WB 2002/09/09)
+       </p>
+
   <li> <p> Fixed: <code class="class">FullMatrix</code>::<code
        class="member">copy_from</code> didn't compile when copying
        from a sparse matrix. This is now fixed.
index 8f8ffc379e1753603a53b00e16af6b149e7b1199..1dafad0e11fe7e05f92afa9d40032efad3a4f3c8 100644 (file)
@@ -17,6 +17,7 @@
 #include "testmatrix.h"
 #include <base/logstream.h>
 #include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
 #include <lac/sparse_matrix_ez.h>
 #include <lac/vector.h>
 #include <lac/solver_richardson.h>
@@ -93,6 +94,49 @@ check_vmult_quadratic(std::vector<double>& residuals,
 }
 
 
+
+void
+check_vmult_quadratic(std::vector<double>& residuals,
+                     const BlockSparseMatrix<double> & A,
+                     const char* prefix)
+{
+  deallog.push(prefix);
+  
+  Vector<double> u(A.n());
+  Vector<double> f(A.m());
+  GrowingVectorMemory<> mem;
+
+  SolverControl control(10, 1.e-13, false);
+  SolverRichardson<> rich(control, mem, .01);
+  SolverRichardson<> prich(control, mem, 1.);
+  PreconditionIdentity identity;
+  PreconditionJacobi<BlockSparseMatrix<double> > jacobi;
+  jacobi.initialize(A, .5);
+
+  PreconditionBlock<BlockSparseMatrix<double>, float>::AdditionalData
+    data((unsigned int) std::sqrt(A.n()+.3), 1.2);
+  
+  PreconditionBlockJacobi<BlockSparseMatrix<double>, float> block_jacobi;
+  block_jacobi.initialize(A, data);
+  
+  u = 0.;
+  f = 1.;
+
+  PREC_CHECK(rich, solve, identity);
+  PREC_CHECK(prich, solve, jacobi);
+  u = 0.;
+  PREC_CHECK(prich, solve, block_jacobi);
+  
+  u = 0.;
+  deallog << "Transpose" << std::endl;
+  PREC_CHECK(rich, Tsolve, identity);
+  PREC_CHECK(prich, Tsolve, jacobi);
+  u = 0.;
+  PREC_CHECK(prich, Tsolve, block_jacobi);
+  deallog.pop();
+}
+
+
 template <class MATRIX>
 void
 check_iterator (const MATRIX& A)
@@ -238,43 +282,71 @@ int main()
 
   std::vector<double> A_res;
   std::vector<double> E_res;
-  
-  deallog << "Structure" << std::endl;
+
+                                  // usual sparse matrix
   SparsityPattern structure(dim, dim, 5);
-  testproblem.five_point_structure(structure);
-  structure.compress();
-  SparseMatrix<double>  A(structure);
-  deallog << "Assemble" << std::endl;
-  testproblem.five_point(A, true);
-  check_vmult_quadratic(A_res, A, "5-SparseMatrix<double>");
+  SparseMatrix<double>  A;
+  {
+    deallog << "Structure" << std::endl;
+    testproblem.five_point_structure(structure);
+    structure.compress();
+    A.reinit (structure);
+    deallog << "Assemble" << std::endl;
+    testproblem.five_point(A, true);
+    check_vmult_quadratic(A_res, A, "5-SparseMatrix<double>");
+#ifdef DEBUG
+    check_iterator(A);
+#endif
+  }
+
+                                  // block sparse matrix with only
+                                  // one block
+  {
+    deallog << "Structure" << std::endl;
+    BlockSparsityPattern block_structure(1,1);
+    block_structure.block(0,0).reinit (dim, dim, 5);
+    block_structure.collect_sizes ();
+    testproblem.five_point_structure(block_structure.block(0,0));
+    block_structure.compress();
+    BlockSparseMatrix<double>  block_A(block_structure);
+    deallog << "Assemble" << std::endl;
+    testproblem.five_point(block_A, true);
+    std::vector<double> block_A_res;
+    check_vmult_quadratic(block_A_res, block_A, "5-BlockSparseMatrix<double>");
 #ifdef DEBUG
-  check_iterator(A);
+    check_iterator(block_A);
 #endif
+  }
 
+                                  // ez sparse matrix
   SparseMatrixEZ<double> E(dim,dim,row_length,2);
-  deallog << "Assemble" << std::endl;
-  testproblem.five_point(E, true);
-  check_vmult_quadratic(E_res, E, "5-SparseMatrixEZ<double>");
+  {
+    deallog << "Assemble" << std::endl;
+    testproblem.five_point(E, true);
+    check_vmult_quadratic(E_res, E, "5-SparseMatrixEZ<double>");
 #ifdef DEBUG
-  check_iterator(E);
+    check_iterator(E);
 #endif
-  E.print_statistics(deallog, true);
-  E.add_scaled(-1., A);
-  if (E.l2_norm() < 1.e-14)
-    deallog << "Matrices are equal" << std::endl;
-  else
-    deallog << "Matrices differ!!" << std::endl;
+    E.print_statistics(deallog, true);
+    E.add_scaled(-1., A);
+    if (E.l2_norm() < 1.e-14)
+      deallog << "Matrices are equal" << std::endl;
+    else
+      deallog << "Matrices differ!!" << std::endl;
+  }
+  
+  {
+    structure.reinit(dim, dim, 5);
+    deallog << "Structure" << std::endl;
+    structure.reinit(dim, dim, 9);
+    testproblem.nine_point_structure(structure);
+    structure.compress();
+    A.reinit (structure);
+    deallog << "Assemble" << std::endl;
+    testproblem.nine_point(A);
+    check_vmult_quadratic(A_res, A, "9-SparseMatrix<double>");
+  }
   
-  A.clear();
-  deallog << "Structure" << std::endl;
-  structure.reinit(dim, dim, 9);
-  testproblem.nine_point_structure(structure);
-  structure.compress();
-  A.reinit(structure);
-  deallog << "Assemble" << std::endl;
-  testproblem.nine_point(A);
-  check_vmult_quadratic(A_res, A, "9-SparseMatrix<double>");
-
   E.clear();
   E.reinit(dim,dim,row_length,2);
   deallog << "Assemble" << std::endl;
index ad5ff5f51caa13f8bde21b16d760f94bdb6dc21b..d22d4f25f6f1703146cb5944ce8b74d26a5da82f 100644 (file)
@@ -14,6 +14,7 @@
 
 #include "testmatrix.h"
 #include <lac/sparse_matrix.h>
+#include <lac/block_sparse_matrix.h>
 #include <lac/sparse_matrix_ez.h>
 #include <lac/vector.h>
 
@@ -249,6 +250,8 @@ FDMatrix::gnuplot_print(std::ostream& s, const Vector<number>& V) const
 
 template void FDMatrix::five_point(SparseMatrix<long double>&, bool) const;
 template void FDMatrix::five_point(SparseMatrix<double>&, bool) const;
+template void FDMatrix::five_point(BlockSparseMatrix<long double>&, bool) const;
+template void FDMatrix::five_point(BlockSparseMatrix<double>&, bool) const;
 template void FDMatrix::five_point(SparseMatrixEZ<double>&, bool) const;
 template void FDMatrix::five_point(SparseMatrix<float>&, bool) const;
 template void FDMatrix::nine_point(SparseMatrix<long double>&, bool) const;

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.