7#include <dune/common/ftraits.hh>
23 template <
class Operator,
class Precoditioner,
class DiscreteFunction>
25 Precoditioner* preconditioner,
26 std::vector< DiscreteFunction >& tempMem,
28 const DiscreteFunction& b,
30 const int maxIterations,
31 const int toleranceCriteria,
32 std::ostream* os =
nullptr )
34 typedef typename DiscreteFunction::RangeFieldType RangeFieldType;
35 typedef typename Dune::FieldTraits< RangeFieldType >::real_type RealType;
37 const RealType tolerance = (epsilon * epsilon) * b.normSquaredDofs( );
39 assert( preconditioner ? tempMem.size() == 5 : tempMem.size() == 3 );
41 DiscreteFunction& h = tempMem[ 0 ];
47 DiscreteFunction& r = tempMem[ 1 ];
52 DiscreteFunction& p = tempMem[ 2 ];
56 DiscreteFunction& s = ( preconditioner ) ? tempMem[ 3 ] : p;
59 DiscreteFunction& q = ( preconditioner ) ? tempMem[ 4 ] : p;
62 (*preconditioner)( p, q );
66 RangeFieldType prevResiduum = 0;
67 RangeFieldType residuum = p.scalarProductDofs( q );
70 for( iterations = 0; (
std::real(residuum) > tolerance) && (iterations < maxIterations); ++iterations )
74 assert( residuum/prevResiduum == residuum/prevResiduum );
75 const RangeFieldType beta= residuum / prevResiduum;
89 RangeFieldType qdoth = q.scalarProductDofs( h );
90 const RangeFieldType alpha = residuum / qdoth;
91 assert( !std::isnan( alpha ) );
97 (*preconditioner)( p, s);
99 prevResiduum = residuum;
100 residuum = p.scalarProductDofs( s );
106 prevResiduum = residuum;
107 residuum = r.normSquaredDofs( );
112 (*os) <<
"Fem::CG it: " << iterations <<
" : sqr(residuum) " << residuum << std::endl;
116 return (iterations < maxIterations) ? iterations : -iterations;
double real(const complex< Dune::Fem::Double > &x)
Definition: double.hh:983
Definition: bindguard.hh:11
int cg(Operator &op, Precoditioner *preconditioner, std::vector< DiscreteFunction > &tempMem, DiscreteFunction &x, const DiscreteFunction &b, const double epsilon, const int maxIterations, const int toleranceCriteria, std::ostream *os=nullptr)
Definition: cg.hh:24
abstract operator
Definition: operator.hh:34
static const int absolute
Definition: cg.hh:17
static const int relative
Definition: cg.hh:18
static const int residualReduction
Definition: cg.hh:19