]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix SparseILU by a patch by Oliver Rheinbach <or@winfos.com>. Add a new test to check...
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 8 Feb 2001 14:35:57 +0000 (14:35 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Thu, 8 Feb 2001 14:35:57 +0000 (14:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@3872 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/news/2000/3.0.0-vs-3.1.0.html
deal.II/lac/include/lac/sparse_ilu.h
deal.II/lac/include/lac/sparse_ilu.templates.h
tests/lac/Makefile.in
tests/lac/sparse_ilu.cc [new file with mode: 0644]
tests/lac/sparse_ilu.checked [new file with mode: 0644]

index f21b2852ebf726681fa0eef8d0613d7b743065d5..1e48070b648cf2ec74ba698fb7f66cb2d0b4dfc5 100644 (file)
@@ -279,6 +279,12 @@ documentation, etc</a>.
 <h3>lac</h3>
 
 <ol>
+  <li> <p>
+       Fixed: <code class="class">SparseILU</code> had two bugs which
+       are now fixed.
+       <br>
+       (<a href="mailto:or@winfos.com">Oliver Rheinbach</a>, 2001/02/02)
+       </p>
 
   <li> <p>
        Improved: block classes have a variable number of blocks now,
index 15f1d3990e0768fccb2a547916415dc92821c4d5..b568c769da0a8ef5acff3dca1256d25282f1766d 100644 (file)
@@ -164,16 +164,20 @@ class SparseILU : protected SparseMatrix<number>
     void reinit (const SparsityPattern &sparsity);
 
                                     /**
-                                     * Perform the incomplete LU factorization
-                                     * of the given matrix.
+                                     * Perform the incomplete LU
+                                     * factorization of the given
+                                     * matrix.
                                      *
-                                     * Note that the sparsity structures of
-                                     * the decomposition and the matrix passed
-                                     * to this function need not be equal,
-                                     * but that the pattern used by this
-                                     * matrix needs to contain all elements
-                                     * used by the matrix to be decomposed.
-                                     * Fill-in is thus allowed.
+                                     * Note that the sparsity
+                                     * structures of the
+                                     * decomposition and the matrix
+                                     * passed to this function need
+                                     * not be equal, but that the
+                                     * pattern used by this matrix
+                                     * needs to contain all elements
+                                     * used by the matrix to be
+                                     * decomposed.  Fill-in is thus
+                                     * allowed.
                                      */
     template <typename somenumber>
     void decompose (const SparseMatrix<somenumber> &matrix,
@@ -189,7 +193,9 @@ class SparseILU : protected SparseMatrix<number>
                              const Vector<somenumber> &src) const;
 
                                     /**
-                                     * Same as @p{apply_decomposition}, format for LAC.
+                                     * Same as
+                                     * @p{apply_decomposition},
+                                     * format for LAC.
                                      */
     template <typename somenumber>
     void vmult (Vector<somenumber>       &dst,
index 753f3b40934aca9c0b233ec906182a90574bf680..c77bd18ac5751d617591d79d6d5f79162396a565 100644 (file)
@@ -176,7 +176,12 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                           // columns are sorted within each
                                           // row correctly, but excluding
                                           // the main diagonal entry
-         bool left_of_diagonal = true;
+         const int global_index_ki = sparsity(*col_ptr,row);
+
+         if (global_index_ki != -1)
+           diag_element(row) -= global_entry(global_index_ik) *
+                                global_entry(global_index_ki);
+
          for (const unsigned int * j = col_ptr+1;
               j<&column_numbers[rowstart_indices[row+1]];
               ++j)
@@ -191,26 +196,6 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
                                               // row linearly. I just didn't
                                               // have the time to figure out
                                               // the details.
-             
-                                              // check whether we have just
-                                              // traversed the diagonal
-                                              //
-                                              // note that j should never point
-                                              // to the diagonal itself!
-             if (left_of_diagonal && (*j > row))
-               {
-                 Assert (*j != row, ExcInternalError());
-                 
-                 left_of_diagonal = false;
-                                                  // a[i,i] -= a[i,k]*a[k,i]
-                 const int global_index_ki = sparsity(*col_ptr,row);
-
-                 if (global_index_ki != -1)
-                   diag_element(row) -= global_entry(global_index_ik) *
-                                        global_entry(global_index_ki);
-                 
-               };
-             
                      const int global_index_ij = j - &column_numbers[0],
                        global_index_kj = sparsity(*col_ptr,*j);
              if ((global_index_ij != -1) &&
@@ -221,6 +206,12 @@ void SparseILU<number>::decompose (const SparseMatrix<somenumber> &matrix,
        };
     };
 
+                                  // Here the very last diagonal
+                                  // element still has to be inverted
+                                  // because the for-loop doesn't do
+                                  // it...
+ diag_element(m()-1)=1./diag_element(m()-1);
+
 /*
   OLD CODE, rather crude first implementation with an algorithm taken
   from 'W. Hackbusch, G. Wittum: Incomplete Decompositions (ILU)-
index 96d4a33308b1dd903ac786c8cb05bfc52278829d..234fee5215521cb45f06c478f5d57f9dba533980 100644 (file)
@@ -1,6 +1,6 @@
 ############################################################
 # $Id$
-# Copyright (C) 2000 by the deal.II authors
+# Copyright (C) 2000, 2001 by the deal.II authors
 ############################################################
 
 ############################################################
@@ -141,6 +141,22 @@ solver.testcase: $(solver-o-files) $(libraries)
 ############################################################
 
 
+sparse_ilu-cc-files = sparse_ilu.cc testmatrix.cc
+
+ifeq ($(debug-mode),on)
+sparse_ilu-o-files = $(sparse_ilu-cc-files:.cc=.go)
+else
+sparse_ilu-o-files = $(sparse_ilu-cc-files:.cc=.o)
+endif
+
+sparse_ilu.testcase: $(sparse_ilu-o-files) $(libraries)
+       @echo =====linking======= $<
+       @$(CXX) $(LDFLAGS) -g  -o $@ $^ $(LIBS)
+
+
+############################################################
+
+
 full_matrix-cc-files = full_matrix.cc
 
 ifeq ($(debug-mode),on)
diff --git a/tests/lac/sparse_ilu.cc b/tests/lac/sparse_ilu.cc
new file mode 100644 (file)
index 0000000..d9cf0a6
--- /dev/null
@@ -0,0 +1,121 @@
+//----------------------------  solver.cc  ---------------------------
+//    $Id$
+//    Version: $Name$
+//
+//    Copyright (C) 2001 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.
+//
+//----------------------------  solver.cc  ---------------------------
+
+
+// make sure that the SparseILU applied with infinite fill-in
+// generates the exact inverse matrix
+
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <cstdlib>
+#include "testmatrix.h"
+#include <base/logstream.h>
+#include <lac/sparse_matrix.h>
+#include <lac/sparse_ilu.h>
+#include <lac/vector.h>
+
+
+
+int main()
+{
+  ofstream logfile("sparse_ilu.output");
+  logfile.setf(ios::fixed);
+  logfile.precision(3);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  
+
+  for (unsigned int size=4; size <= 16; size *= 2)
+    {
+      unsigned int dim = (size-1)*(size-1);
+
+      deallog << "Size " << size << " Unknowns " << dim << endl;
+      
+                                      // Make matrix
+      FDMatrix testproblem(size, size);
+      SparsityPattern structure(dim, dim, 5);
+      testproblem.build_structure(structure);
+      structure.compress();
+      SparseMatrix<double>  A(structure);
+      testproblem.laplacian(A);
+
+      
+      for (unsigned int test=0; test<2; ++test)
+       {
+         deallog << "Test " << test << endl;
+         
+                                          // generate sparse ILU.
+                                          //
+                                          // for test 1, test with
+                                          // full pattern.  for test
+                                          // 2, test with same
+                                          // pattern as A
+         SparsityPattern ilu_pattern (dim, dim,
+                                      (test==0 ? dim : 5));
+         switch (test)
+           {
+             case 0:
+                   for (unsigned int i=0; i<dim; ++i)
+                     for (unsigned int j=0; j<dim; ++j)
+                       ilu_pattern.add(i,j);
+                   break;
+
+             case 1:
+                   for (unsigned int i=0; i<dim; ++i)
+                     for (unsigned int j=0; j<dim; ++j)
+                       if (structure(i,j) != SparsityPattern::invalid_entry)
+                         ilu_pattern.add(i,j);
+                   break;
+
+             default:
+                   Assert (false, ExcNotImplemented());
+           };
+         ilu_pattern.compress();
+         SparseILU<double> ilu (ilu_pattern);
+         ilu.decompose (A);
+         
+                                          // now for three test vectors v
+                                          // determine norm of
+                                          // (I-BA)v, where B is ILU of A.
+                                          // since matrix is symmetric,
+                                          // likewise test for right
+                                          // preconditioner
+         Vector<double> v(dim);
+         Vector<double> tmp1(dim), tmp2(dim);
+         for (unsigned int i=0; i<3; ++i)
+           {
+             for (unsigned int j=0; j<dim; ++j)
+               v(j) = 1. * rand()/RAND_MAX;
+             
+             A.vmult (tmp1, v);
+             ilu.vmult (tmp2, tmp1);
+             tmp2 -= v;
+             const double left_residual = tmp2.l2_norm();
+             
+             ilu.vmult (tmp1, v);
+             A.vmult (tmp2, tmp1);
+             tmp2 -= v;
+             const double right_residual = tmp2.l2_norm();
+             
+             
+             deallog << "Residual with test vector " << i << ":  "
+                     << " left=" << left_residual
+                     << ", right=" << right_residual
+                     << endl;
+           };
+       };
+      
+    };
+};
+
diff --git a/tests/lac/sparse_ilu.checked b/tests/lac/sparse_ilu.checked
new file mode 100644 (file)
index 0000000..420bfc0
--- /dev/null
@@ -0,0 +1,28 @@
+JobId Thu Feb  8 15:21:09 2001
+DEAL::Size 4 Unknowns 9
+DEAL::Test 0
+DEAL::Residual with test vector 0:   left=0.000, right=0.000
+DEAL::Residual with test vector 1:   left=0.000, right=0.000
+DEAL::Residual with test vector 2:   left=0.000, right=0.000
+DEAL::Test 1
+DEAL::Residual with test vector 0:   left=0.425, right=0.447
+DEAL::Residual with test vector 1:   left=0.266, right=0.267
+DEAL::Residual with test vector 2:   left=0.297, right=0.377
+DEAL::Size 8 Unknowns 49
+DEAL::Test 0
+DEAL::Residual with test vector 0:   left=0.000, right=0.000
+DEAL::Residual with test vector 1:   left=0.000, right=0.000
+DEAL::Residual with test vector 2:   left=0.000, right=0.000
+DEAL::Test 1
+DEAL::Residual with test vector 0:   left=1.895, right=1.900
+DEAL::Residual with test vector 1:   left=1.705, right=1.711
+DEAL::Residual with test vector 2:   left=2.061, right=2.057
+DEAL::Size 16 Unknowns 225
+DEAL::Test 0
+DEAL::Residual with test vector 0:   left=0.000, right=0.000
+DEAL::Residual with test vector 1:   left=0.000, right=0.000
+DEAL::Residual with test vector 2:   left=0.000, right=0.000
+DEAL::Test 1
+DEAL::Residual with test vector 0:   left=5.900, right=5.911
+DEAL::Residual with test vector 1:   left=5.659, right=5.662
+DEAL::Residual with test vector 2:   left=5.914, right=5.921

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.