[Trilinos-Users] using teko

[Hoang Giang Bui]

Hi

This seems working. The preconditioner by Teko could compute the Schur
complement and the Aztec solver works. However, there is one issue which I
did not fully understand:

I created Teko linear operator for each of my blocks by:

Teko::LinearOp block_a  = Thyra::epetraLinearOp(Teuchos::rcp(&(*a)));
Teko::LinearOp block_b1 = Thyra::epetraLinearOp(Teuchos::rcp(&(*b1)));
Teko::LinearOp block_b2 = Thyra::epetraLinearOp(Teuchos::rcp(&(*b2)));
Teko::LinearOp block_c  = Thyra::epetraLinearOp(Teuchos::rcp(&(*c)));

Where a, b1, b2, c are of boost::shared_ptr< Epetra_FECrsMatrix >

creating the block:

Teuchos::RCP<Teko::Epetra::EpetraOperatorWrapper> P = Teuchos::rcp(new
Teko::Epetra::EpetraOperatorWrapper(Teko::block2x2(block_a,block_b1,block_b2,block_c)));

creating the factory:
Teuchos::RCP<Teko::BlockPreconditionerFactory> precFact = Teuchos::rcp(new
MultiphaseBlock2x2PreconditionerFactory(inverseA, inverseS));

build the preconditioner:

Teko::Epetra::EpetraBlockPreconditioner prec(precFact);
prec.buildPreconditioner(P);


pass to aztec and solve:
....
           aztec_solver.SetPrecOperator(&prec);
           aztec_solver.Iterate(mmax_iter, mtol);
....

However after solving, I tried to probe the value of matrix 'a' by
std::cout << (*a)[0][0] leads to a segmentation fault error. This prevents
me to propagate to next analysis step. I'm worry that the declaration
Teko::LinearOp block_a  = Thyra::epetraLinearOp(Teuchos::rcp(&(*a)) has
some implications. Can you advise me if there is anything wrong with the
scheme above ?

BR
Giang Bui








On Thu, Jan 3, 2013 at 11:52 PM, Cyr, Eric C <eccyr at sandia.gov> wrote:

>  Giang Bui,
>
>  Thanks for the question. I'll try to explain what is going on.
>
>     const Teko::LinearOp invA = Teko::buildInverse(*invAfact_, A);
>    const Teko::LinearOp invS = Teko::buildInverse(*invSfact_, S);
>
>  These operations create some approximation inverse of 'A' and 'S'. None
> of the standard inverse factories create a concrete Epetra_CrsMatrix (they
> are possibly Epetra_Operator objects). Now when you call
>
>     const Teko::LinearOp tmpP = Teko::explicitMultiply(invS,
> Teko::explicitMultiply(B2, invA));
>
>  What this says is that you want to form the product of invS*B2*invA
> explicitly. That means
> it will try to do first multiply B2*invA and store it as a CRS matrix.
> Generally this operation is not possible
> since B2 is probably a sparse operator and invA is likely an approximation
> of a dense
> inverse (using Aztec for instance). Teko will not allow this operation
> because it would be
> very, very expensive. What Teko does allow is "implicit" multiply
> operations. That is you
> can compute the action of invS*B2*invA by chaining together their
> application. To do that
> the pertinent line is
>
>     const Teko::LinearOp tmpP = Teko::multiply(invS, Teko::multiply(B2,
> invA));
>
>  So what is explicitMultiply and explicitAdd for then? Well sometimes you
> know you have
> a concrete sparse matrix and you want to explicitly form an operator so
> you can apply
> multigrid or ILU to it as an approximate inverse. For example, in the
> SIMPLE preconditioner
> for Navier-Stokes you need to form a Schur complement that looks like
>
>     S=C-B*H*Bt
>
>  Where C, B, and Bt are sparse CRS matrices and H is a diagonal matrix.
> Thus the Teko code
> for combining those is (not the most efficient version mind you)
>
>     const Teko::LinearOp S =
> Teko::explicitAdd(C,Teko::scale(-1.0,Teko::explicitMulitiply(B,H,Bt)));
>
>  This gives a concrete sparse CRS matrix that you can then apply ML,
> IFPACK or Amesos to. The
> efficient construction of these approximations is what makes block
> preconditioning challenging,
> but this is an issue that has to be resolved and understood at the
> mathematical level.
>
>  Please let me know if some I said is not clear. I appreciate your
> interest in Teko and would be eager
> to hear about your project and its goals, so please feel free to contact
> me.
>
>  Eric
>
>   From: "Heroux, Michael A" <maherou at sandia.gov>
> Date: Thu, 3 Jan 2013 15:30:26 -0700
> To: Eric Cyr <eccyr at sandia.gov>, Hoang Giang Bui <hgbk2008 at gmail.com>
> Subject: Re: [Trilinos-Users] using teko
>
>   Giang Bui,
>
>  I have copied Eric Cyr on this message.  He is the principal author of
> Teko.
>
>  Eric, Do you have any advice?
>
>  Thanks.
>
>  Mike
>
>   From: Hoang Giang Bui <hgbk2008 at gmail.com>
> Date: Tuesday, December 25, 2012 7:42 AM
> To: "trilinos-users at software.sandia.gov" <
> trilinos-users at software.sandia.gov>
> Subject: [EXTERNAL] [Trilinos-Users] using teko
>
>      Hi
>
>  When using Teko to build the preconditioner operator for Aztec solver I
> faced this error:
>
>  I tried to create the LinearOp of invS * B2 * invA. The invS is the
> approximation of S^{-1} and invA is approximation of A^{-1}. I created invS
> and invA by this code
>
> const Teko::LinearOp invA = Teko::buildInverse(*invAfact_, A);
> const Teko::LinearOp invS = Teko::buildInverse(*invSfact_, S);
>
>  Then
>
> const Teko::LinearOp tmpP = Teko::explicitMultiply(invS,
> Teko::explicitMultiply(B2, invA));
>
>  At this point, the program crashed and thow the error:
>
> Throw number = 1
>
> Throw test that evaluated to true: true
>
> dyn_cast<Thyra::EpetraLinearOp>(Thyra::LinearOpBase<double>) : Error, the
> object with the concrete type 'Thyra::DefaultInverseLinearOp<double>'
> (passed in through the interface type 'Thyra::LinearOpBase<double>')  does
> not support the interface 'Thyra::EpetraLinearOp' and the dynamic cast
> failed!
>
>  As I assumed, the type of invA and invS is
> Thyra::DefaultInverseLinearOp<double> and not inherited from
> Thyra::EpetraLinearOp, hence the dynamic_cast is failed. I also think that
> the InverseLinearOp does not try to create Crs_Matrix so it is not
> compatible with EpetraLinearOp.
>
>  Please help me on this case. It's important for me to create a block
> preconditioner for my problem.
>
> BR
> Giang Bui
>
>


-- 
With Best Regards !
Giang Bui
To learn and to excel
-------------- next part --------------
An HTML attachment was scrubbed...
URL: https://software.sandia.gov/pipermail/trilinos-users/attachments/20130104/9881718c/attachment.html