From 8568c5a0dbb7300a9e2a38adb226983400e7f1f8 Mon Sep 17 00:00:00 2001
From: Peter Bastian <peter@dune-project.org>
Date: Tue, 5 Oct 2004 11:40:35 +0000
Subject: [PATCH] testing of generic solver
[[Imported from SVN: r33]]
---
istl/tutorial/example.cc | 94 +++++++++++++++++++++++++++++++++++++---
1 file changed, 89 insertions(+), 5 deletions(-)
diff --git a/istl/tutorial/example.cc b/istl/tutorial/example.cc
index 59efab25f..60fdeb47d 100644
--- a/istl/tutorial/example.cc
+++ b/istl/tutorial/example.cc
@@ -17,6 +17,7 @@
#include "fmatrix.hh"
#include "bcrsmatrix.hh"
#include "io.hh"
+#include "gsetc.hh"
// a simple stop watch
class Timer
@@ -40,6 +41,12 @@ public:
cend = times(&buf);
return ((double)(cend-cstart))/100.0;
}
+
+ double gettime ()
+ {
+ return ((double)(cend-cstart))/100.0;
+ }
+
private:
clock_t cstart,cend;
};
@@ -214,11 +221,11 @@ void test_FieldMatrix ()
A /= 3.14;
A.umv(z,b);
- A.umvt(b,z);
- A.umvh(b,z);
+ A.umtv(b,z);
+ A.umhv(b,z);
A.usmv(-1.0,z,b);
- A.usmvt(-1.0,b,z);
- A.usmvh(-1.0,b,z);
+ A.usmtv(-1.0,b,z);
+ A.usmhv(-1.0,b,z);
std::cout << A.frobenius_norm() << " " << A.frobenius_norm2() << std::endl;
std::cout << A.infinity_norm() << " " << A.infinity_norm_real() << std::endl;
@@ -312,6 +319,78 @@ void test_IO ()
printmatrix(std::cout,B,"a block compressed sparse matrix","row",9,1);
}
+void test_Iter ()
+{
+ // block types
+ Timer t;
+ const int BlockSize=1;
+ typedef Dune::FieldVector<double,BlockSize> VB;
+ typedef Dune::FieldMatrix<double,BlockSize,BlockSize> MB;
+
+ // a fake discretization
+ t.start();
+ const int N=1000000, K1=1, K2=1000;
+ Dune::BCRSMatrix<MB> A(N,N,5*N,Dune::BCRSMatrix<MB>::row_wise);
+ for (Dune::BCRSMatrix<MB>::CreateIterator i=A.createbegin(); i!=A.createend(); ++i)
+ {
+ i.insert(i.index());
+ if (i.index()-K1>=0) i.insert(i.index()-K1);
+ if (i.index()-K2>=0) i.insert(i.index()-K2);
+ if (i.index()+K2< N) i.insert(i.index()+K2);
+ if (i.index()+K1< N) i.insert(i.index()+K1);
+ }
+ for (Dune::BCRSMatrix<MB>::RowIterator i=A.begin(); i!=A.end(); ++i)
+ for (Dune::BCRSMatrix<MB>::ColIterator j=(*i).begin(); j!=(*i).end(); ++j)
+ if (i.index()==j.index())
+ (*j) = 4;
+ else
+ (*j) = -1;
+ t.stop();
+ std::cout << "time for build=" << t.gettime() << " seconds." << std::endl;
+ // printmatrix(std::cout,A,"system matrix","row",10,2);
+
+ // set up system
+ Dune::BlockVector<VB> x(N),b(N),d(N);
+ b=0; b[0] = 1; b[N-1] = 1;
+ // printvector(std::cout,b,"rhs","entry",1,10,2);
+
+ // same in defect formulation
+ x = 0; // initial guess
+ std::cout.setf(std::ios_base::scientific, std::ios_base::floatfield);
+ std::cout.precision(8);
+ t.start();
+ d=b; A.mmv(x,d); // compute defect
+ std::cout << 0 << " " << d.two_norm() << std::endl;
+ Dune::BlockVector<VB> v(x); // memory for update
+ v = 1.23E-4;
+ double w=1.0; // damping factor
+ for (int k=1; k<=50; k++)
+ {
+ bgs_update(A,v,d); // compute update
+ x.axpy(w,v); // update solution
+ A.usmv(-w,v,d); // update defect
+ std::cout << k << " " << d.two_norm() << std::endl;
+ }
+ t.stop();
+ std::cout << "time for solve=" << t.gettime() << " seconds." << std::endl;
+ return;
+
+ // printvector(std::cout,x,"solution","entry",1,10,2);
+
+ // a simple iteration
+ x = 0; // initial guess
+ d=b; A.mmv(x,d);
+ std::cout.setf(std::ios_base::scientific, std::ios_base::floatfield);
+ std::cout.precision(8);
+ std::cout << 0 << " " << d.two_norm() << std::endl;
+ for (int k=1; k<=50; k++)
+ {
+ bsor(A,x,b,1.0);
+ d=b; A.mmv(x,d);
+ std::cout << k << " " << d.two_norm() << std::endl;
+ }
+ // printvector(std::cout,x,"solution","entry",1,10,2);
+}
int main (int argc , char ** argv)
{
@@ -321,7 +400,12 @@ int main (int argc , char ** argv)
// test_VariableBlockVector();
// test_FieldMatrix();
// test_BCRSMatrix();
- test_IO();
+ // test_IO();
+ test_Iter();
+ }
+ catch (Dune::ISTLError& error)
+ {
+ std::cout << error << std::endl;
}
catch (Dune::Exception& error)
{
--
GitLab