ILU(n) works (as far as I could test it). Unfortunately only

n=0 and n=1 are practically useful. n=2 is already a direct
solver. (May I have still something wrong ...)

[[Imported from SVN: r978]]

istl/ilu.hh

@@ -160,7 +160,7 @@ namespace Dune {
     createiterator ci=ILU.createbegin();
     for (crowiterator i=A.begin(); i!=endi; ++i)
     {
-      std::cout << "symbolic factorization, row " << i.index() << std::endl;
+      //		std::cout << "in row " << i.index() << std::endl;
       map rowpattern;           // maps column index to generation
 
       // initialize pattern with row of A
@@ -170,22 +170,27 @@ namespace Dune {
       // eliminate entries in row which are to the left of the diagonal
       for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
       {
-        std::cout << "eliminating " << i.index() << "," << (*ik).first
-                  << std::endl;
-        coliterator endk = ILU[(*ik).first].end();                 // end of row k
-        coliterator kj = ILU[(*ik).first].find((*ik).first);                 // diagonal in k
-        for (++kj; kj!=endk; ++kj)                 // row k eliminates in row i
+        if ((*ik).second<n)
         {
-          int generation = static_cast<int>(firstmatrixelement(*kj));
-          mapiterator ij = rowpattern.find(kj.index());
-          if (ij==rowpattern.end())
-          {                         // new entry
+          //                            std::cout << "  eliminating " << i.index() << "," << (*ik).first
+          //                                              << " level " << (*ik).second << std::endl;
+
+          coliterator endk = ILU[(*ik).first].end();                       // end of row k
+          coliterator kj = ILU[(*ik).first].find((*ik).first);                       // diagonal in k
+          for (++kj; kj!=endk; ++kj)                       // row k eliminates in row i
+          {
+            int generation = (int) firstmatrixelement(*kj);
             if (generation<n)
-              rowpattern[kj.index()] = generation+1;
-          }
-          else
-          {                         // entry exists already
-            (*ij).second = std::min((*ij).second,generation+1);
+            {
+              mapiterator ij = rowpattern.find(kj.index());
+              if (ij==rowpattern.end())
+              {
+                //std::cout << "    new entry " << i.index() << "," << kj.index()
+                //                                                << " level " << (*ik).second+1 << std::endl;
+
+                rowpattern[kj.index()] = generation+1;
+              }
+            }
           }
         }
       }
@@ -195,17 +200,19 @@ namespace Dune {
         ci.insert((*ik).first);
       ++ci;           // now row i exist
 
-      // set generation index on new row
-      coliterator endik=ILU[i.index()].end();
-      for (coliterator ik=ILU[i.index()].begin(); ik!=endik; ++ik)
-        firstmatrixelement(*ik) = static_cast<K>(rowpattern[ik.index()]);
+      // write generation index into entries
+      coliterator endILUij = ILU[i.index()].end();;
+      for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
+        firstmatrixelement(*ILUij) = (K) rowpattern[ILUij.index()];
     }
 
-    // initialize new matrix with A
+    //	printmatrix(std::cout,ILU,"ilu pattern","row",10,2);
+
+    // copy entries of A
     for (crowiterator i=A.begin(); i!=endi; ++i)
     {
-      coliterator ILUij = ILU[i.index()].begin();
-      coliterator endILUij;
+      coliterator ILUij;
+      coliterator endILUij = ILU[i.index()].end();;
       for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
         (*ILUij) = 0;           // clear row
       ccoliterator Aij = (*i).begin();
@@ -214,7 +221,10 @@ namespace Dune {
       while (Aij!=endAij && ILUij!=endILUij)
       {
         if (Aij.index()==ILUij.index())
+        {
           *ILUij = *Aij;
+          ++Aij; ++ILUij;
+        }
         else
         {
           if (Aij.index()<ILUij.index())

istl/tutorial/example.cc

@@ -435,7 +435,7 @@ void test_Iter ()
 void test_Interface ()
 {
   // define Types
-  const int BlockSize = 1;
+  const int BlockSize = 6;
   typedef Dune::FieldVector<double,BlockSize> VB;
   typedef Dune::FieldMatrix<double,BlockSize,BlockSize> MB;
   typedef Dune::BlockVector<VB> Vector;
@@ -452,7 +452,7 @@ void test_Interface ()
     E[i][i] = -1;
 
   // make a block compressed row matrix with five point stencil
-  const int BW2=64, BW1=1, N=BW2*BW2;
+  const int BW2=100, BW1=1, N=BW2*BW2;
   Matrix A(N,N,5*N,Dune::BCRSMatrix<MB>::row_wise);
   for (Matrix::CreateIterator i=A.createbegin(); i!=A.createend(); ++i)
   {
@@ -469,6 +469,8 @@ void test_Interface ()
       else
         (*j) = E;
 
+  //  printmatrix(std::cout,A,"system matrix","row",10,2);
+
   // set up system
   Vector x(N),b(N);
   x=0; x[0]=1; x[N-1]=2; // prescribe known solution
@@ -477,10 +479,12 @@ void test_Interface ()
 
   // set up the high-level solver objects
   Dune::MatrixAdapter<Matrix,Vector,Vector> op(A);        // make linear operator from A
-  Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,1,1.78);     // SSOR preconditioner
+  Dune::SeqSSOR<Matrix,Vector,Vector> ssor(A,1,1);     // SSOR preconditioner
   Dune::SeqILU0<Matrix,Vector,Vector> ilu0(A);            // preconditioner object
-  Dune::CGSolver<Vector> cg(op,ilu0,1E-8,8000,2);         // an inverse operator
-  Dune::BiCGSTABSolver<Vector> bcgs(op,ilu0,1E-8,8000,2); // an inverse operator
+  Dune::SeqILUn<Matrix,Vector,Vector> ilun(A,1);          // preconditioner object
+  Dune::LoopSolver<Vector> loop(op,ilun,1E-8,8000,2);     // an inverse operator
+  Dune::CGSolver<Vector> cg(op,ilun,1E-8,8000,2);         // an inverse operator
+  Dune::BiCGSTABSolver<Vector> bcgs(op,ilun,1E-8,8000,2); // an inverse operator
 
   // call the solver
   Dune::InverseOperatorResult r;