<h3>Specific improvements</h3>
<ol>
+<li> Improved: The SolverCG implementation now uses only three auxiliary
+vectors, down from previously four. Also, there are some shortcuts in case
+PreconditionIdentity is used that improve the solver's performance.
+<br>
+(Martin Kronbichler, 2013/04/11)
+</li>
+
<li> Fixed: The results section of step-23 did not show the movie in release 7.3
due to a poor HTML markup. This is now fixed.
<br>
DEAL_II_NAMESPACE_OPEN
+// forward declaration
+class PreconditionIdentity;
+
+
/*!@addtogroup Solvers */
/*@{*/
VECTOR *Vr;
VECTOR *Vp;
VECTOR *Vz;
- VECTOR *VAp;
/**
* Within the iteration loop, the
this->memory.free(Vr);
this->memory.free(Vp);
this->memory.free(Vz);
- this->memory.free(VAp);
deallog.pop();
}
deallog.push("cg");
// Memory allocation
- Vr = this->memory.alloc();
- Vp = this->memory.alloc();
- Vz = this->memory.alloc();
- VAp = this->memory.alloc();
+ Vr = this->memory.alloc();
+ Vz = this->memory.alloc();
+ Vp = this->memory.alloc();
// Should we build the matrix for
// eigenvalue computations?
bool do_eigenvalues = additional_data.compute_condition_number
try
{
// define some aliases for simpler access
- VECTOR &g = *Vr;
- VECTOR &h = *Vp;
- VECTOR &d = *Vz;
- VECTOR &Ad = *VAp;
+ VECTOR &g = *Vr;
+ VECTOR &d = *Vz;
+ VECTOR &h = *Vp;
// resize the vectors, but do not set
// the values since they'd be overwritten
// soon anyway.
g.reinit(x, true);
- h.reinit(x, true);
d.reinit(x, true);
- Ad.reinit(x, true);
+ h.reinit(x, true);
// Implementation taken from the DEAL
// library
int it=0;
return;
}
- precondition.vmult(h,g);
+ if (types_are_equal<PRECONDITIONER,PreconditionIdentity>::value == false)
+ {
+ precondition.vmult(h,g);
- d.equ(-1.,h);
+ d.equ(-1.,h);
- gh = g*h;
+ gh = g*h;
+ }
+ else
+ {
+ d.equ(-1.,g);
+ gh = res*res;
+ }
while (conv == SolverControl::iterate)
{
it++;
- A.vmult(Ad,d);
+ A.vmult(h,d);
- alpha = d*Ad;
+ alpha = d*h;
Assert(alpha != 0., ExcDivideByZero());
alpha = gh/alpha;
- g.add(alpha,Ad);
- x.add(alpha,d );
+ g.add(alpha,h);
+ x.add(alpha,d);
res = g.l2_norm();
print_vectors(it, x, g, d);
if (conv != SolverControl::iterate)
break;
- precondition.vmult(h,g);
+ if (types_are_equal<PRECONDITIONER,PreconditionIdentity>::value
+ == false)
+ {
+ precondition.vmult(h,g);
- beta = gh;
- Assert(beta != 0., ExcDivideByZero());
- gh = g*h;
- beta = gh/beta;
+ beta = gh;
+ Assert(beta != 0., ExcDivideByZero());
+ gh = g*h;
+ beta = gh/beta;
+ d.sadd(beta,-1.,h);
+ }
+ else
+ {
+ beta = gh;
+ gh = res*res;
+ beta = gh/beta;
+ d.sadd(beta,-1.,g);
+ }
if (additional_data.log_coefficients)
deallog << "alpha-beta:" << alpha << '\t' << beta << std::endl;
deallog << "Condition number estimate: " <<
T.eigenvalue(T.n()-1)/T.eigenvalue(0) << std::endl;
}
-
- d.sadd(beta,-1.,h);
}
}
catch (...)