[Dune] Performace of ISTL

[Markus Blatt]

Hallo,

I do not totally understand your code, but it seems to me that you
solve different systems with trilinos and istl:

- You seem to have 3 levels of recursiveness in your ISTL vectors
  (like in BlockVector<BlockVector<FieldVector<double,1> > >, but when
  copying it to Trilinos you simply disregard the first level and just 
  copy vector[1]. Therefore it seems like you only solve a part of
  DUNE problem in Trilinos
- After the computation you copy your solution back to ISTL and use it
  as an initial guess. Therefore your initial defect is already much
  smaller than the one you had with Trilinos.
  You should really use the same initial guess for both solvers.
-
On Thu, Apr 19, 2007 at 04:51:52PM +0200, Patrick Leidenberger wrote:
> 
> [  .... Code ... ] 
> #ifdef TRILINOS_SOLVER
>  #ifdef HAVE_TRILINOS
>   
>    int vecSize = yi[1].size();
>    Epetra_Map map(vecSize, 0, Comm_);
> 
>    /** Create trilinos vector. */
>    Epetra_Vector yiTrilinos_(map);
>    Epetra_Vector xiTrilinos_(map);
> 
>    /** Copy from dune to trilinos vector. */
>    for (unsigned int vecIter = 0;  vecIter < vecSize; ++vecIter)
>    {
>      yiTrilinos_[vecIter] = yi[1][vecIter][0];
>      xiTrilinos_[vecIter] =  EmpClass::dof1_[1][vecIter][0];
                               ^^^^^^^^^^^^^^^^^^
Here you diregard the first level of the hierarchy.

>    }
>    Epetra_LinearProblem problem(&(*EmpClass::aijTrilinos_), 
> &xiTrilinos_, &yiTrilinos_);
>    AztecOO solver(problem);
>    int options[AZ_OPTIONS_SIZE];
>    double params[AZ_PARAMS_SIZE];
>    AZ_defaults(options, params);
> 
>    /** Solve the linear system. */
>    
> solver.Iterate(EmpClass::solverMaximumIteration_,EmpClass::solverReduction_);
> 
>    /** Copy the solution vector from trilinos to a dune vector. */
>    for (unsigned int vecIter = 0;  vecIter < vecSize; ++vecIter)
>    {
>      EmpClass::dofp1_[1][vecIter][0] = xiTrilinos_[vecIter];
>    }

why do you copy it back?

>
> [ .... Code ... ]
>  typedef Dune::MatrixAdapter<EmpMatrixTopo,EmpVectorTopo,EmpVectorTopo> 
> Operator;
>  Operator op(EmpClass::aij_);
>

Here you give dune all levels of the hierarchy and do not disregard
the first recusion level.

>
> [ ... Code ... ]
>

cheers,

Markus