From 490e8a2b5c945d0334dc68b109be6d9b3046d0de Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Thu, 8 Feb 2001 14:35:57 +0000 Subject: [PATCH] Fix SparseILU by a patch by Oliver Rheinbach . Add a new test to check for exact inversion in case we allow the ILU to use all elements of the matrix. Also check that the residual is small in case we use a rather sparse pattern for the ILU (the same one in fact as for the original five-point matrix. 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 | 6 + deal.II/lac/include/lac/sparse_ilu.h | 26 ++-- .../lac/include/lac/sparse_ilu.templates.h | 33 ++--- tests/lac/Makefile.in | 18 ++- tests/lac/sparse_ilu.cc | 121 ++++++++++++++++++ tests/lac/sparse_ilu.checked | 28 ++++ 6 files changed, 200 insertions(+), 32 deletions(-) create mode 100644 tests/lac/sparse_ilu.cc create mode 100644 tests/lac/sparse_ilu.checked diff --git a/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html b/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html index f21b2852eb..1e48070b64 100644 --- a/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html +++ b/deal.II/doc/news/2000/3.0.0-vs-3.1.0.html @@ -279,6 +279,12 @@ documentation, etc.

lac

    +
  1. + Fixed: SparseILU had two bugs which + are now fixed. +
    + (Oliver Rheinbach, 2001/02/02) +

  2. Improved: block classes have a variable number of blocks now, diff --git a/deal.II/lac/include/lac/sparse_ilu.h b/deal.II/lac/include/lac/sparse_ilu.h index 15f1d3990e..b568c769da 100644 --- a/deal.II/lac/include/lac/sparse_ilu.h +++ b/deal.II/lac/include/lac/sparse_ilu.h @@ -164,16 +164,20 @@ class SparseILU : protected SparseMatrix 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 void decompose (const SparseMatrix &matrix, @@ -189,7 +193,9 @@ class SparseILU : protected SparseMatrix const Vector &src) const; /** - * Same as @p{apply_decomposition}, format for LAC. + * Same as + * @p{apply_decomposition}, + * format for LAC. */ template void vmult (Vector &dst, diff --git a/deal.II/lac/include/lac/sparse_ilu.templates.h b/deal.II/lac/include/lac/sparse_ilu.templates.h index 753f3b4093..c77bd18ac5 100644 --- a/deal.II/lac/include/lac/sparse_ilu.templates.h +++ b/deal.II/lac/include/lac/sparse_ilu.templates.h @@ -176,7 +176,12 @@ void SparseILU::decompose (const SparseMatrix &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::decompose (const SparseMatrix &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::decompose (const SparseMatrix &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)- diff --git a/tests/lac/Makefile.in b/tests/lac/Makefile.in index 96d4a33308..234fee5215 100644 --- a/tests/lac/Makefile.in +++ b/tests/lac/Makefile.in @@ -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 index 0000000000..d9cf0a67aa --- /dev/null +++ b/tests/lac/sparse_ilu.cc @@ -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 +#include +#include +#include +#include "testmatrix.h" +#include +#include +#include +#include + + + +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 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 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 v(dim); + Vector tmp1(dim), tmp2(dim); + for (unsigned int i=0; i<3; ++i) + { + for (unsigned int j=0; j