From fc241636081fa3323f60008d2ccea1e86955322d Mon Sep 17 00:00:00 2001
From: Peter Bastian <peter@dune-project.org>
Date: Fri, 22 Oct 2004 17:18:51 +0000
Subject: [PATCH] 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 | 54 ++++++++++++++++++++++++----------------
istl/tutorial/example.cc | 14 +++++++----
2 files changed, 41 insertions(+), 27 deletions(-)
diff --git a/istl/ilu.hh b/istl/ilu.hh
index 7c65a46d4..3cbab2e7e 100644
--- a/istl/ilu.hh
+++ b/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())
diff --git a/istl/tutorial/example.cc b/istl/tutorial/example.cc
index 65a4ef3e0..ac76300a0 100644
--- a/istl/tutorial/example.cc
+++ b/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;
--
GitLab