]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
SparseMatrixEZ used in tests
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Jul 2002 17:50:13 +0000 (17:50 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Mon, 15 Jul 2002 17:50:13 +0000 (17:50 +0000)
git-svn-id: https://svn.dealii.org/trunk@6246 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/lac/include/lac/sparse_matrix_ez.templates.h
tests/Makefile.rules
tests/lac/Makefile
tests/lac/solver.cc
tests/lac/sparse_ilu.cc
tests/lac/sparse_matrices.cc [new file with mode: 0644]
tests/lac/sparse_matrices.checked [new file with mode: 0644]
tests/lac/sparsity_pattern.cc
tests/lac/testmatrix.cc
tests/lac/testmatrix.h

index 81ca366be0d221f5e68a00c8474b8f43dea74732..60f3c0a988550316e126c649c2334bab39ff74d3 100644 (file)
@@ -8,7 +8,7 @@ SparseMatrixEZ<number>::Row::set(unsigned int column,
                                 const number& value)
 {
                                   // Store end of vector
-  const std::vector<Entry>::const_iterator end_col = end();
+  const typename std::vector<Entry>::iterator end_col = end();
 
                                   // Create Entry for inserting
   const Entry e(column, value);
@@ -16,7 +16,7 @@ SparseMatrixEZ<number>::Row::set(unsigned int column,
                                   // Find position for inserting
                                   // should return first Entry with
                                   // higher column number.
-  std::vector<Entry>::iterator col = lower_bound(begin(), end, e);
+  typename std::vector<Entry>::iterator col = lower_bound(begin(), end_col, e);
 
                                   // Insert Entry if column did not exist.
                                   // Edit existing entry otherwise.
@@ -24,7 +24,7 @@ SparseMatrixEZ<number>::Row::set(unsigned int column,
     values.push_back(e);
   else
     if (col->column == column)
-      col->values = value;
+      col->value = value;
     else
       values.insert(col, e);
 }
@@ -35,15 +35,15 @@ void
 SparseMatrixEZ<number>::Row::add(unsigned int column,
                                 const number& value)
 {
-  const std::vector<Entry>::const_iterator end_col = end();
+  const typename std::vector<Entry>::const_iterator end_col = end();
 
                                   // Create Entry for inserting
-  const Entry e(column, value);
+  const typename SparseMatrixEZ<number>::Entry e(column, value);
   
                                   // Find position for inserting
                                   // should return first Entry with
                                   // higher column number.
-  std::vector<Entry>::iterator col = lower_bound(begin(), end, e);
+  typename std::vector<Entry>::iterator col = lower_bound(begin(), end, e);
 
                                   // Insert Entry if column did not exist.
                                   // Edit existing entry otherwise.
@@ -67,7 +67,7 @@ SparseMatrixEZ<number>::SparseMatrixEZ()
 
 
 template <typename number>
-SparseMatrixEZ<number>::SparseMatrixEZ(const SparseMatrixEZ&)
+SparseMatrixEZ<number>::SparseMatrixEZ(const SparseMatrixEZ<number>&)
 {
   Assert(false, ExcNotImplemented());
 }
@@ -116,7 +116,7 @@ SparseMatrixEZ<number>::clear()
 
 template <typename number>
 bool
-SparseMatrixEZ<number>::empty()
+SparseMatrixEZ<number>::empty() const
 {
   return ((n_columns == 0) && (rows.size()==0));
 }
@@ -131,12 +131,12 @@ SparseMatrixEZ<number>::vmult (Vector<somenumber>& dst,
   Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
 
-  std::vector<Row>::const_iterator row = rows.begin();
-  const std::vector<Row>::const_iterator end_row = rows.end();
+  typename std::vector<Row>::const_iterator row = rows.begin();
+  const typename std::vector<Row>::const_iterator end_row = rows.end();
   for (unsigned int i=0; row != end_row; ++i, ++row)
     {
-      std::vector<Entry>::const_iterator col = row->begin();
-      const std::vector<Entry>::const_iterator end_col = row->end();
+      typename std::vector<Entry>::const_iterator col = row->begin();
+      const typename std::vector<Entry>::const_iterator end_col = row->end();
 
       double s = 0.;
       for (;col != end_col; ++col)
@@ -166,12 +166,12 @@ SparseMatrixEZ<number>::vmult_add (Vector<somenumber>& dst,
   Assert(m() == dst.size(), ExcDimensionMismatch(m(),dst.size()));
   Assert(n() == src.size(), ExcDimensionMismatch(n(),src.size()));
 
-  std::vector<Row>::const_iterator row = rows.begin();
-  const std::vector<Row>::const_iterator end_row = rows.end();
+  typename std::vector<Row>::const_iterator row = rows.begin();
+  const typename std::vector<Row>::const_iterator end_row = rows.end();
   for (unsigned int i=0; row != end_row; ++i, ++row)
     {
-      std::vector<Entry>::const_iterator col = row->begin();
-      const std::vector<Entry>::const_iterator end_col = row->end();
+      typename std::vector<Entry>::const_iterator col = row->begin();
+      const typename std::vector<Entry>::const_iterator end_col = row->end();
 
       double s = 0.;
       for (;col != end_col; ++col)
@@ -189,28 +189,26 @@ SparseMatrixEZ<number>::Tvmult_add (Vector<somenumber>& dst,
   Assert(n() == dst.size(), ExcDimensionMismatch(n(),dst.size()));
   Assert(m() == src.size(), ExcDimensionMismatch(m(),src.size()));
 
-  std::vector<Row>::const_iterator row = rows.begin();
-  const std::vector<Row>::const_iterator end_row = rows.end();
+  typename std::vector<Row>::const_iterator row = rows.begin();
+  const typename std::vector<Row>::const_iterator end_row = rows.end();
   for (unsigned int i=0; row != end_row; ++i, ++row)
     {
-      std::vector<Entry>::const_iterator col = row->begin();
-      const std::vector<Entry>::const_iterator end_col = row->end();
+      typename std::vector<Entry>::const_iterator col = row->begin();
+      const typename std::vector<Entry>::const_iterator end_col = row->end();
 
-      double s = 0.;
       for (;col != end_col; ++col)
        dst(col->column) += col->value * src(i);
     }
 }
 
 template <typename number>
-template <typename somenumber>
 unsigned int
 SparseMatrixEZ<number>::memory_consumption() const
 {
   unsigned int result = sizeof (*this)
                        + sizeof(Row) * sizeof (rows);
 
-  for (std::vector<Row>::const_iterator r = rows.begin();
+  for (typename std::vector<Row>::const_iterator r = rows.begin();
        r != rows.end(); ++r)
     result += r->size() * sizeof(Entry);
 
index b4726d7dde957b861f8c0cabde76674c66ea2b51..91927dd36acbe18c9cb49a75a99b5309531268ff 100644 (file)
@@ -27,6 +27,10 @@ endif
 %.go : %.cc Makefile
        @echo =====debug========= $<
        @$(CXX) $(flags) -c $< -o $@
+
+%.o : %.cc Makefile
+       @echo =====optimized===== $<
+       @$(CXX) $(CXXFLAGS.o) -c $< -o $@
 %.exe :
        @echo =====linking======= $@
        @$(CXX) $(LDFLAGS) -o $@ $^ $(LIBS)
index 19ab30819abd20aa97df23822e0eb568f208a57d..63f23b2ce6495c29c0cc81d4a50dbc236add10d4 100644 (file)
@@ -1,8 +1,6 @@
-# Generated automatically from Makefile.in by configure.
-# Generated automatically from Makefile.in by configure.
 ############################################################
 # $Id$
-# Copyright (C) 2000, 2001 by the deal.II authors
+# Copyright (C) 2000, 2001, 2002 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -15,6 +13,8 @@ debug-mode = on
 
 libraries = $(lib-lac.g)  \
             $(lib-base.g)
+libraries_o = $(lib-lac.o)  \
+            $(lib-base.o)
 
 default: run-tests 
 
@@ -26,6 +26,7 @@ block_vector_iterator.exe : block_vector_iterator.go                        $(li
 full_matrix.exe           : full_matrix.go                                  $(libraries)
 matrix_out.exe            : matrix_out.go                                   $(libraries)
 solver.exe                : solver.go                         testmatrix.go $(libraries)
+sparse_matrices.exe       : sparse_matrices.go                testmatrix.go $(libraries)
 sparsity_pattern.exe      : sparsity_pattern.go               testmatrix.go $(libraries)
 sparse_ilu.exe            : sparse_ilu.go                     testmatrix.go $(libraries)
 vector-vector.exe         : vector-vector.go                                $(libraries)
@@ -33,7 +34,8 @@ vector-vector.exe         : vector-vector.go                                $(li
 #mg.exe                   : mg.go                             testmatrix.go $(libraries)
 
 
-tests = block_matrices block_vector block_vector_iterator full_matrix \
+tests = sparse_matrices block_matrices block_vector block_vector_iterator \
+        full_matrix \
         matrix_out solver sparsity_pattern sparse_ilu vector-vector
 
 ############################################################
index f5520064639b9d46b7c805ed8fbe6526cb93af4d..5d495b85999e7f44012a62e3b329b013c252caf5 100644 (file)
@@ -72,10 +72,10 @@ int main()
                                       // Make matrix
       FDMatrix testproblem(size, size);
       SparsityPattern structure(dim, dim, 5);
-      testproblem.build_structure(structure);
+      testproblem.five_point_structure(structure);
       structure.compress();
       SparseMatrix<double>  A(structure);
-      testproblem.laplacian(A);
+      testproblem.five_point(A);
 
       PreconditionIdentity prec_no;
       PreconditionSOR<> prec_sor;
index 723ec0743bc89f0eda8b24d4ff64644ad2ffef4a..9a683024895a414fdf9e65ec822f69e06b6078fc 100644 (file)
@@ -46,10 +46,10 @@ int main()
                                       // Make matrix
       FDMatrix testproblem(size, size);
       SparsityPattern structure(dim, dim, 5);
-      testproblem.build_structure(structure);
+      testproblem.five_point_structure(structure);
       structure.compress();
       SparseMatrix<double>  A(structure);
-      testproblem.laplacian(A);
+      testproblem.five_point(A);
 
       
       for (unsigned int test=0; test<2; ++test)
diff --git a/tests/lac/sparse_matrices.cc b/tests/lac/sparse_matrices.cc
new file mode 100644 (file)
index 0000000..77ccb2a
--- /dev/null
@@ -0,0 +1,107 @@
+//----------------------------  testmatrix.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 1998, 1999, 2000, 2002 by the deal.II authors
+//
+//    This file is subject to QPL and may not be  distributed
+//    without copyright and license information. Please refer
+//    to the file deal.II/doc/license.html for the  text  and
+//    further information on this license.
+//
+//----------------------------  testmatrix.cc  ---------------------------
+
+
+#include "testmatrix.h"
+#include <base/logstream.h>
+#include <lac/sparse_matrix.h>
+#include <lac/sparse_matrix_ez.h>
+#include <lac/vector.h>
+#include <lac/solver_richardson.h>
+#include <lac/precondition.h>
+
+#include <lac/sparse_matrix_ez.templates.h>
+
+#include <fstream>
+
+template<class MATRIX>
+void
+check_vmult_quadratic(const MATRIX& A,
+                     const char* prefix)
+{
+  deallog.push(prefix);
+  
+  Vector<double> u(A.n());
+  Vector<double> f(A.m());
+  GrowingVectorMemory<> mem;
+
+  SolverControl control(50, 1.e-3, false);
+  SolverRichardson<> rich(control, mem, .01);
+  PreconditionIdentity prec;
+
+  u = 0.;
+  f = 1.;
+
+  try
+    {
+      rich.solve(A, u, f, prec);
+    }
+  catch (...)
+    {
+    }
+  deallog << "Transpose" << std::endl;
+  try
+    {
+      rich.Tsolve(A, u, f, prec);
+    }
+  catch (...)
+    {
+    }
+  deallog.pop();
+}
+
+int main()
+{
+  std::ofstream logfile("sparse_matrices.output");
+  logfile.setf(std::ios::fixed);
+  logfile.precision(2);
+  deallog.attach(logfile);
+
+                                  // Switch between regression test
+                                  // and benchmark
+#ifdef DEBUG  
+  deallog.depth_console(0);
+  deallog.log_execution_time(true);
+  deallog.log_time_differences(true);
+  const unsigned int size = 1000;
+#else
+  const unsigned int size = 100;
+#endif
+  
+  FDMatrix testproblem (size, size);
+  unsigned int dim = (size-1)*(size-1);
+  
+  SparsityPattern structure(dim, dim, 5);
+  testproblem.five_point_structure(structure);
+  structure.compress();
+  SparseMatrix<double>  A(structure);
+  testproblem.five_point(A);
+  check_vmult_quadratic(A, "5-SparseMatrix<double>");
+
+  SparseMatrixEZ<double> E(dim,dim);
+  testproblem.five_point(E);
+  check_vmult_quadratic(E, "5-SparseMatrixEZ<double>");
+
+  A.clear();
+  structure.reinit(dim, dim, 9);
+  testproblem.nine_point_structure(structure);
+  structure.compress();
+  A.reinit(structure);
+  check_vmult_quadratic(A, "9-SparseMatrix<double>");
+
+  E.clear();
+  E.reinit(dim,dim);
+  testproblem.nine_point(E);
+  check_vmult_quadratic(E, "9-SparseMatrixEZ<double>");
+  
+}
diff --git a/tests/lac/sparse_matrices.checked b/tests/lac/sparse_matrices.checked
new file mode 100644 (file)
index 0000000..85588b3
--- /dev/null
@@ -0,0 +1,29 @@
+
+9.84  :DEAL:5-SparseMatrix<double>:Richardson::Starting 
+52.82 :DEAL:5-SparseMatrix<double>:Richardson::Failure step 50 
+0.00  :DEAL:5-SparseMatrix<double>::Transpose
+1.56  :DEAL:5-SparseMatrix<double>:RichardsonT::Starting 
+80.04 :DEAL:5-SparseMatrix<double>:RichardsonT::Failure step 50 
+0.00  :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00  :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+16.61 :DEAL:5-SparseMatrixEZ<double>:Richardson::Starting 
+67.84 :DEAL:5-SparseMatrixEZ<double>:Richardson::Failure step 50 
+0.00  :DEAL:5-SparseMatrixEZ<double>::Transpose
+1.52  :DEAL:5-SparseMatrixEZ<double>:RichardsonT::Starting 
+78.19 :DEAL:5-SparseMatrixEZ<double>:RichardsonT::Failure step 50 
+0.00  :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00  :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+13.09 :DEAL:9-SparseMatrix<double>:Richardson::Starting 
+73.20 :DEAL:9-SparseMatrix<double>:Richardson::Failure step 50 
+0.00  :DEAL:9-SparseMatrix<double>::Transpose
+2.35  :DEAL:9-SparseMatrix<double>:RichardsonT::Starting 
+119.03:DEAL:9-SparseMatrix<double>:RichardsonT::Failure step 50 
+0.00  :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00  :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
+32.49 :DEAL:9-SparseMatrixEZ<double>:Richardson::Starting 
+94.45 :DEAL:9-SparseMatrixEZ<double>:Richardson::Failure step 50 
+0.00  :DEAL:9-SparseMatrixEZ<double>::Transpose
+2.14  :DEAL:9-SparseMatrixEZ<double>:RichardsonT::Starting 
+109.61:DEAL:9-SparseMatrixEZ<double>:RichardsonT::Failure step 50 
+0.00  :DEAL::GrowingVectorMemory:Overall allocated vectors: 4
+0.00  :DEAL::GrowingVectorMemory:Maximum allocated vectors: 2
index 0de153d6794c250a7ded7be2f9bf358ec5bf9ad7..bc688c6083c4156f30bf9cea99809473c8f73e44 100644 (file)
@@ -33,7 +33,7 @@ main ()
                                   // generate usual sparsity pattern
   const unsigned int N = 15;
   SparsityPattern sp1((N-1)*(N-1), (N-1)*(N-1), 5);
-  FDMatrix(N,N).build_structure (sp1);
+  FDMatrix(N,N).five_point_structure (sp1);
   sp1.compress ();
   deallog << sp1.n_rows() << " " << sp1.n_cols() << " "
          << sp1.bandwidth() << " " << sp1.n_nonzero_elements()
index 7220583398deb4a891558f528ab2e39e1e7e2e02..1f922f2e54e2f001a094a67fd812fc76dd57233f 100644 (file)
@@ -2,7 +2,7 @@
 //    $Id$
 //    Version: $Name$
 //
-//    Copyright (C) 1998, 1999, 2000 by the deal.II authors
+//    Copyright (C) 1998, 1999, 2000, 2002 by the deal.II authors
 //
 //    This file is subject to QPL and may not be  distributed
 //    without copyright and license information. Please refer
 
 #include "testmatrix.h"
 #include <lac/sparse_matrix.h>
+#include <lac/sparse_matrix_ez.h>
 #include <lac/vector.h>
 
+#include <lac/sparse_matrix_ez.templates.h>
+
 FDMatrix::FDMatrix(unsigned int nx, unsigned int ny)
                :
                nx(nx), ny(ny)
 {}
 
 void
-FDMatrix::build_structure(SparsityPattern& structure) const
+FDMatrix::five_point_structure(SparsityPattern& structure) const
 {
   for(unsigned int i=0;i<=ny-2;i++)
     {
@@ -55,9 +58,119 @@ FDMatrix::build_structure(SparsityPattern& structure) const
     }
 }
 
-template<typename number>
 void
-FDMatrix::laplacian(SparseMatrix<number>& A) const
+FDMatrix::nine_point_structure(SparsityPattern& structure) const
+{
+  for(unsigned int i=0;i<=ny-2;i++)
+    {
+      for(unsigned int j=0;j<=nx-2; j++)
+       {
+                                          // Number of the row to be entered
+         unsigned int row = j+(nx-1)*i;
+         structure.add(row, row);
+         if (j>0)
+           {
+             structure.add(row-1, row);
+             structure.add(row, row-1);
+             if (i>0)
+               {
+                 structure.add(row-1, row-(nx-1));
+                 structure.add(row-(nx-1), row-1);
+               }
+             if (i<ny-2)
+               {
+                 structure.add(row-1, row+(nx-1));
+                 structure.add(row+(nx-1), row-1);
+               }
+           }
+         if (j<nx-2)
+           {
+             structure.add(row+1, row);
+             structure.add(row, row+1);
+             if (i>0)
+               {
+                 structure.add(row+1, row-(nx-1));
+                 structure.add(row-(nx-1), row+1);
+               }
+             if (i<ny-2)
+               {
+                 structure.add(row+1, row+(nx-1));
+                 structure.add(row+(nx-1), row+1);
+               }
+           }
+         if (i>0)
+           {
+             structure.add(row-(nx-1), row);
+             structure.add(row, row-(nx-1));
+           }
+         if (i<ny-2)
+           {
+             structure.add(row+(nx-1), row);
+             structure.add(row, row+(nx-1));
+           }
+       }
+    }
+}
+
+template<typename MATRIX>
+void
+FDMatrix::nine_point(MATRIX& A) const
+{
+   for(unsigned int i=0;i<=ny-2;i++)
+    {
+      for(unsigned int j=0;j<=nx-2; j++)
+       {
+                                          // Number of the row to be entered
+         unsigned int row = j+(nx-1)*i;
+         
+         A.set(row, row, 20.);
+         if (j>0)
+           {
+             A.set(row-1, row, -4.);
+             A.set(row, row-1, -4.);
+             if (i>0)
+               {
+                 A.set(row-1, row-(nx-1), -1.);
+                 A.set(row-(nx-1), row-1, -1.);
+               }
+             if (i<ny-2)
+               {
+                 A.set(row-1, row+(nx-1), -1.);
+                 A.set(row+(nx-1), row-1, -1.);
+               }
+           }
+         if (j<nx-2)
+           {
+             A.set(row+1, row, -4.);
+             A.set(row, row+1, -4.);
+             if (i>0)
+               {
+                 A.set(row+1, row-(nx-1), -1.);
+                 A.set(row-(nx-1), row+1, -1.);
+               }
+             if (i<ny-2)
+               {
+                 A.set(row+1, row+(nx-1), -1.);
+                 A.set(row+(nx-1), row+1, -1.);
+               }
+           }
+         if (i>0)
+           {
+             A.set(row-(nx-1), row, -4.);
+             A.set(row, row-(nx-1), -4.);
+           }
+         if (i<ny-2)
+           {
+             A.set(row+(nx-1), row, -4.);
+             A.set(row, row+(nx-1), -4.);
+           }
+       }
+    } 
+}
+
+template<typename MATRIX>
+void
+FDMatrix::five_point(MATRIX& A) const
 {
    for(unsigned int i=0;i<=ny-2;i++)
     {
@@ -110,9 +223,14 @@ FDMatrix::gnuplot_print(std::ostream& s, const Vector<number>& V) const
 
 
 
-template void FDMatrix::laplacian(SparseMatrix<long double>&) const;
-template void FDMatrix::laplacian(SparseMatrix<double>&) const;
-template void FDMatrix::laplacian(SparseMatrix<float>&) const;
+template void FDMatrix::five_point(SparseMatrix<long double>&) const;
+template void FDMatrix::five_point(SparseMatrix<double>&) const;
+template void FDMatrix::five_point(SparseMatrixEZ<double>&) const;
+template void FDMatrix::five_point(SparseMatrix<float>&) const;
+template void FDMatrix::nine_point(SparseMatrix<long double>&) const;
+template void FDMatrix::nine_point(SparseMatrix<double>&) const;
+template void FDMatrix::nine_point(SparseMatrixEZ<double>&) const;
+template void FDMatrix::nine_point(SparseMatrix<float>&) const;
 template void FDMatrix::gnuplot_print(std::ostream& s, const Vector<long double>& V) const;
 template void FDMatrix::gnuplot_print(std::ostream& s, const Vector<double>& V) const;
 template void FDMatrix::gnuplot_print(std::ostream& s, const Vector<float>& V) const;
index 78952e3b3c690fb44b64a52ef194857cc95b51fc..d671d9c4dd0705d45dbba1481ae27ee98eef466b 100644 (file)
@@ -19,13 +19,24 @@ class FDMatrix
                                     /**
                                      * Generate the matrix structure.
                                      */
-    void build_structure(SparsityPattern& structure) const;
+    void five_point_structure(SparsityPattern& structure) const;
+    
+                                    /**
+                                     * Generate the matrix structure.
+                                     */
+    void nine_point_structure(SparsityPattern& structure) const;
     
                                     /**
                                      * Fill the matrix with values.
                                      */
-    template <typename number>
-    void laplacian(SparseMatrix<number>&) const;
+    template <typename MATRIX>
+    void five_point(MATRIX&) const;
+
+                                    /**
+                                     * Fill the matrix with values.
+                                     */
+    template <typename MATRIX>
+    void nine_point(MATRIX&) const;
 
     template <typename number>
     void gnuplot_print(std::ostream&, const Vector<number>&) 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.