--- /dev/null
+//---------------------------- sparse_ilu.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2001, 2002, 2003, 2004, 2005, 2008, 2013 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.
+//
+//---------------------------- sparse_ilu.cc ---------------------------
+
+
+// make sure that the SparseMIC applied with infinite fill-in
+// generates the exact inverse matrix
+
+#include "../tests.h"
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <iomanip>
+#include <cstdlib>
+#include "testmatrix.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/sparse_mic.h>
+#include <deal.II/lac/vector.h>
+
+//TODO:[WB] find test that is less sensitive to floating point accuracy
+
+int main()
+{
+ std::ofstream logfile("sparse_mic/output");
+ deallog << std::fixed;
+ deallog << std::setprecision(3);
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+
+ for (unsigned int size=4; size <= 16; size *= 2)
+ {
+ unsigned int dim = (size-1)*(size-1);
+
+ deallog << "Size " << size << " Unknowns " << dim << std::endl;
+
+ // Make matrix
+ FDMatrix testproblem(size, size);
+ SparsityPattern structure(dim, dim, dim);
+ // testproblem.five_point_structure(structure);
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ structure.add(i,j);
+ structure.compress();
+ SparseMatrix<double> A(structure);
+ testproblem.five_point(A);
+
+
+ for (unsigned int test=0; test<2; ++test)
+ {
+ deallog << "Test " << test << std::endl;
+
+ // generate sparse ILU.
+ //
+ // for test 1, test with
+ // full pattern. for test
+ // 2, test with same
+ // pattern as A
+ SparsityPattern mic_pattern (dim, dim, dim);
+ switch (test)
+ {
+ case 0:
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ mic_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)
+ mic_pattern.add(i,j);
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ };
+ mic_pattern.compress();
+ SparseMIC<double> mic (mic_pattern);
+ mic.decompose (A);
+
+ // now for three test vectors v
+ // determine norm of
+ // (I-BA)v, where B is MIC 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. * std::rand()/RAND_MAX;
+
+ A.vmult (tmp1, v);
+ mic.vmult (tmp2, tmp1);
+ tmp2 -= v;
+ const double left_residual = tmp2.l2_norm() / v.l2_norm();
+
+ mic.vmult (tmp1, v);
+ A.vmult (tmp2, tmp1);
+ tmp2 -= v;
+ const double right_residual = tmp2.l2_norm() / v.l2_norm();
+
+
+ deallog << "Relative residual with test vector " << i << ": "
+ << " left=" << left_residual
+ << ", right=" << right_residual
+ << std::endl;
+ };
+ };
+
+ };
+}
--- /dev/null
+
+DEAL::Size 4 Unknowns 9
+DEAL::Test 0
+DEAL::Relative residual with test vector 0: left=0.071, right=0.170
+DEAL::Relative residual with test vector 1: left=0.043, right=0.097
+DEAL::Relative residual with test vector 2: left=0.073, right=0.146
+DEAL::Test 1
+DEAL::Relative residual with test vector 0: left=0.093, right=0.210
+DEAL::Relative residual with test vector 1: left=0.029, right=0.159
+DEAL::Relative residual with test vector 2: left=0.110, right=0.192
+DEAL::Size 8 Unknowns 49
+DEAL::Test 0
+DEAL::Relative residual with test vector 0: left=0.154, right=0.659
+DEAL::Relative residual with test vector 1: left=0.204, right=0.561
+DEAL::Relative residual with test vector 2: left=0.166, right=0.651
+DEAL::Test 1
+DEAL::Relative residual with test vector 0: left=0.230, right=0.649
+DEAL::Relative residual with test vector 1: left=0.148, right=0.629
+DEAL::Relative residual with test vector 2: left=0.190, right=0.702
+DEAL::Size 16 Unknowns 225
+DEAL::Test 0
+DEAL::Relative residual with test vector 0: left=0.305, right=1.589
+DEAL::Relative residual with test vector 1: left=0.427, right=1.498
+DEAL::Relative residual with test vector 2: left=0.333, right=1.572
+DEAL::Test 1
+DEAL::Relative residual with test vector 0: left=0.278, right=1.428
+DEAL::Relative residual with test vector 1: left=0.336, right=1.512
+DEAL::Relative residual with test vector 2: left=0.390, right=1.476
--- /dev/null
+
+DEAL::Size 4 Unknowns 9
+DEAL::Test 0
+DEAL::Residual with test vector 0: left=0, right=0
+DEAL::Residual with test vector 1: left=0, right=0
+DEAL::Residual with test vector 2: left=0, right=0
+DEAL::Test 1
+DEAL::Residual with test vector 0: left=0.295, right=0.333
+DEAL::Residual with test vector 1: left=0.383, right=0.385
+DEAL::Residual with test vector 2: left=0.299, right=0.285
+DEAL::Size 8 Unknowns 49
+DEAL::Test 0
+DEAL::Residual with test vector 0: left=0, right=0
+DEAL::Residual with test vector 1: left=0, right=0
+DEAL::Residual with test vector 2: left=0, right=0
+DEAL::Test 1
+DEAL::Residual with test vector 0: left=2.243, right=2.263
+DEAL::Residual with test vector 1: left=1.707, right=1.707
+DEAL::Residual with test vector 2: left=1.768, right=1.788
+DEAL::Size 16 Unknowns 225
+DEAL::Test 0
+DEAL::Residual with test vector 0: left=0, right=0
+DEAL::Residual with test vector 1: left=0, right=0
+DEAL::Residual with test vector 2: left=0, right=0
+DEAL::Test 1
+DEAL::Residual with test vector 0: left=5.672, right=5.683
+DEAL::Residual with test vector 1: left=5.931, right=5.944
+DEAL::Residual with test vector 2: left=5.988, right=5.996