[Fwd: Re: [Dune] Performace of ISTL]

[Patrick Leidenberger]

Hello Markus

> 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> > >,
exact
>   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
At the moment I we have only 1st order base functions, so we get only
one degree of freedom per edge.
If this code runs we want to use higher order base functions, which will
lead to degrees of freedom per face and element, more than 1 degree of
freedom per entity.

The matrix is able to handle this different degrees of freedom, due to
the hierarchical matrix structure. But at the moment I do not use this
feature of my code.

To make a quick and dirty test with trilinos (I don't know if trilinos
supports hierarchical matrices like dune), I copy the matrix entires not
with general code, but I copy all entries!



> - 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.
That is not the point! I tell you the times of the first solving step.
Both vectors are initialized with zero! By the way, a few weeks ago I
tested the dune solver with an initial guess, but I do not save one
iteration step, there was no speed up. Does the initial guess in Dune works?

Is there a loss of speed if I use a hierarchical matrix with istl? Will
it be better to copy the hierarchical matrix into a non hierarchical matrix?

Perhaps I will test this next week.

> -
> 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.
Yes, but without any loss of information. Only [1][i][0] entries are set.
> 
>>    }
>>    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?
I make a lot of calculations with the dofp1_ in other places of the
code. Trilinos is only for testing in the code, generally I want use the
istl solver. And I want to create as less redundant code as possible.
> 
>> [ .... 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.
.. see above.

Best regards

Patrick