#206 Problem using solver with hierarchical vector and matrix
Metadata
| Property | Value |
|---|---|
| Reported by | Patrick Leidenberger (mailings.pl@gmail.com) |
| Reported at | Nov 23, 2006 12:05 |
| Type | Bug Report |
| Version | Git (pre2.4) [autotools] |
| Operating System | Unspecified / All |
| Closed by | Markus Blatt (markus@dr-blatt.de) |
| Closed at | Nov 24, 2006 08:48 |
| Closed in version | Unknown |
| Resolution | Not a bug |
| Comment | If we need to have a direct solver for BCRSMatrix, please open a Feature Request. |
Description
Hi all,
I want to use a istl solver for a hierarchical matrix and vector structure, but I get an error from dune when I call the solver (dbjac(H ,v,d,w);)
I attached the modified example.cc from dune-istl/istl/tutorial and the error message.
Is this an error form dune or is my code wrong?
Best regards
Patrick
Error message:
if g++ -DHAVE_CONFIG_H -I. -I. -I../.. -I/home2/leidenberger/afs/private/dune/dune-common -I../.. -g -O2 -MT example.o -MD -MP -MF ".deps/example.Tpo" -c -o example.o example.cc;
then mv -f ".deps/example.Tpo" ".deps/example.Po"; else rm -f ".deps/example.Tpo"; exit 1; fi
../../dune/istl/gsetc.hh: In static member function static void Dune::algmeta_itsteps<0>::dbjac(const M&, X&, const Y&, const K&) [with M = Dune::BCRSMatrix<Dune::FieldMatrix<test_Iter::ValueType, 1, 1>, Dune::ISTLAllocator>, X = Dune::BlockVector<Dune::FieldVector<test_Iter::ValueType, 1>, Dune::ISTLAllocator>, Y = test_Iter()::EmpVectorTopoElemOrder, K = double]': ../../dune/istl/gsetc.hh:472: instantiated from static void Dune::algmeta_itsteps::dbjac(const M&, X&, const Y&, const K&) [with M = Dune::BCRSMatrix<test_Iter()::EmpMatrixTopoElemOrder, Dune::ISTLAllocator>, X = Dune::BlockVector<test_Iter()::EmpVectorTopoElemOrder, Dune::ISTLAllocator>, Y = Dune::BlockVector<test_Iter()::EmpVectorTopoElemOrder, Dune::ISTLAllocator>, K = double, int I = 1]'
../../dune/istl/gsetc.hh:545: instantiated from `void Dune::dbjac(const M&, X&, const Y&, const K&) [with M = test_Iter()::EmpMatrixTopoElem, X = test_Iter()::EmpVectorTopoElem, Y = test_Iter()::EmpVectorTopoElem, K = double]'
example.cc:119: instantiated from here
../../dune/istl/gsetc.hh:498: error: 'const class Dune::BCRSMatrix<Dune::FieldMatrix<test_Iter::ValueType, 1, 1>, Dune::ISTLAllocator>' has no member named 'solve'
make: *** [example.o] Error 1
example.cc:
// start with including some headers #include "config.h"
#include // for input/output to shell #include // for input/output to files #include // STL vector class #include
#include<math.h> // Yes, we do some math here #include<stdio.h> // There is nothing better than sprintf #include<sys/times.h> // for timing measurements
#include<dune/istl/istlexception.hh> #include<dune/istl/basearray.hh> #include<dune/common/fvector.hh> #include<dune/common/fmatrix.hh> #include<dune/istl/bvector.hh> #include<dune/istl/vbvector.hh> #include<dune/istl/bcrsmatrix.hh> #include<dune/istl/io.hh> #include<dune/istl/gsetc.hh> #include<dune/istl/ilu.hh> #include<dune/istl/operators.hh> #include<dune/istl/solvers.hh> #include<dune/istl/preconditioners.hh> #include<dune/istl/scalarproducts.hh>
void test_Iter () {
const int GG=2, HH =3,; typedef double ValueType; typedef Dune::BCRSMatrix<Dune::FieldMatrix<ValueType,1,1> > EmpMatrixTopoElemOrder; typedef Dune::BCRSMatrix EmpMatrixTopoElem; // typedef Dune::BCRSMatrix EmpMatrixTopo;
typedef Dune::BlockVector<Dune::FieldVector<ValueType,1> > EmpVectorTopoElemOrder; typedef Dune::BlockVector EmpVectorTopoElem; // typedef Dune::VariableBlockVector EmpVectorTopo;
// Create matrix. Dune::FieldMatrix<ValueType,1,1> F; F = 42.0;
EmpMatrixTopoElemOrder G(GG,GG,GG*GG,EmpMatrixTopoElemOrder::row_wise); for (EmpMatrixTopoElemOrder::CreateIterator i=G.createbegin(); i!=G.createend(); ++i) for (int j=0; j<GG; ++j) i.insert(j);
for (EmpMatrixTopoElemOrder::RowIterator i=G.begin(); i!=G.end(); ++i) for (EmpMatrixTopoElemOrder::ColIterator j=(*i).begin(); j!=(*i).end(); ++j) *j = F;
// std::cout << "============================================." << std::endl; // printmatrix(std::cout,G,"G","row",10,2); // std::cout << "=============================================."<< std::endl;
EmpMatrixTopoElem H(HH,HH,HH,EmpMatrixTopoElem::row_wise); for (EmpMatrixTopoElem::CreateIterator i=H.createbegin(); i!=H.createend(); ++i) for (int j=0; j<HH; ++j) { if (i.index() == j) { i.insert(j); } }
for (EmpMatrixTopoElem::RowIterator i=H.begin(); i!=H.end(); ++i) for (EmpMatrixTopoElem::ColIterator j=(*i).begin(); j!=(*i).end(); ++j) *j = G;
std::cout << "============================================." << std::endl; printmatrix(std::cout,H,"H","row",10,2); std::cout << "=============================================."<< std::endl;
EmpVectorTopoElem x,b,d,v; x.resize(HH);
EmpVectorTopoElem::iterator iter = x.begin();
for(iter; iter != x.end(); ++iter)
{
std::cout << "Hallo" << std::endl;
iter->resize(GG);
}
x = 0;
x[0][0] = 1.0;
x[1][0] = 21.0;
x[2][0] = 42.0;
b = x;
d = x;
v = x; // Vector for update.
printvector(std::cout,x,"exact solution","entry",10,10,2); b=0; H.umv(x,b); // set right hand side x=0; // initial guess
// solve in defect formulation std::cout.setf(std::ios_base::scientific, std::ios_base::floatfield); std::cout.precision(8);
d=b; H.mmv(x,d); // compute defect
double w=1.0; // damping factor
EmpMatrixTopoElem ILU(H);
// bilu0_decomposition(ILU); // printmatrix(std::cout,ILU,"ilu decomposition","row",12,4); for (int k=1; k<=20; k++) { v=0; // bilu_backsolve(ILU,v,d); // dbgs(H,v,d,w); // compute update dbjac(H ,v,d,w); // compute update // bsorf(H,v,d,w); // compute update // bsorb(H,v,d,w); // compute update x += v; // update solution H.mmv(v,d); // update defect // bltsolve(H,v,d,w); // compute update // x += v; // update solution // H.mmv(v,d); // update defect // butsolve(H,v,d,w); // compute update // x += v; // update solution // H.mmv(v,d); // update defect std::cout << k << " " << d.two_norm() << std::endl; if (d.two_norm()<1E-4) break; } }
int main (int argc , char ** argv) { try { test_Iter(); } catch (Dune::ISTLError& error) { std::cout << error << std::endl; } catch (Dune::Exception& error) { std::cout << error << std::endl; } catch (const std::bad_alloc& e) { std::cout << "memory exhausted" << std::endl; } catch (...) { std::cout << "unknown exception caught" << std::endl; }
return 0; }