<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,
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,
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,
// 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)
// 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) &&
};
};
+ // 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)-
############################################################
# $Id$
-# Copyright (C) 2000 by the deal.II authors
+# Copyright (C) 2000, 2001 by the deal.II authors
############################################################
############################################################
############################################################
+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)
--- /dev/null
+//---------------------------- 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;
+ };
+ };
+
+ };
+};
+
--- /dev/null
+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