[Dune] Problem with CG Solver

[Markus Blatt]

Hi,

On Tue, Nov 22, 2011 at 04:50:49PM +0100, S. Swayamjyoti wrote:
> Please find below my solver, could some one please point out why the
> residuals are not going to zero? Is it because of bad conditioning of the
> matrix?
> 
> *template<class GV>
> void P1Elements<GV>::solve_u_i()
> {
>    // make linear operator from A
>     Dune::MatrixAdapter<Matrix,ScalarField,ScalarField>
> op(ThomasFermiMatrix);
> 
>     // initialize preconditioner
>     Dune::SeqILUn<Matrix,ScalarField,ScalarField> ilu2(ThomasFermiMatrix,
> 1, 0.92)
> 
>     // the inverse operator
>     Dune::CGSolver<ScalarField> cg(op, ilu2, 1e-15, 5000, 2);
>     Dune::InverseOperatorResult r1;
> 
>     // initialue u to some arbitrary value to avoid u being the exact
>     // solution
>     u1.resize(b1.N(), false);
>     u1 = 0.0;
> 
>         b1 *= -1.0;
>     // finally solve the system
>     cg.apply(u1, b1, r1);
> } *
> 


that is hard to say. Might be the condition, but I would rather guess
that your b1 is not in the range of the matrix (i.e. there is no x,
such that Ax=b).

You could e.g. check whether this is a problem of the matrix condition
or cg by using a random left hand side x and then first compute the
right hand side b=A*x before solving A*u1=b.

Cheers,

Markus
-- 
Do you need more support with DUNE or HPC in general? 

Dr. Markus Blatt - HPC-Simulation-Software & Services http://www.dr-blatt.de
Rappoltsweilerstr. 5, 68229 Mannheim, Germany
Tel.: +49 (0) 160 97590858  Fax: +49 (0)322 1108991658