[Trilinos-Users] [EXTERNAL] Re: using teko

Cyr, Eric C eccyr at sandia.gov
Fri Jan 4 10:20:57 MST 2013


Giang Bui

I think that the rcp's associated with the teko linear operator are going
out of scope. Thus when the 'a' (boost::shared_ptr) is accessed the underlying
epetra operator has been deleted by the Teuchos::RCP. This is easily
resolved by ensuring the RCP does not delete the operator when it goes
out of scope and gets deconstructed. You have a few options in this regard.

I like this one (though others don't seem to use it as much as me)

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

This builds an RCP that understands that it _does_not_ own the underlying pointer
and thus should not try to delete it when the reference count goes to zero. It
does however imply that the user is responsible for make sure that the reference
stays in scope (is not deallocated) for as long as the RCP is required. An option with
a different statement but with the same meaning is to do

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

This says the rcp does not own the memory. A final option, which I have never
seen used but think it might be fun to try is

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));

This uses the boost::shared_ptr directly, not sure if the pointer is shared
between the RCP and boost::shared_ptr in some intelligent way or not.

Anyway, I see that you are looking at Multi Phase, I would be interested
to see what you are doing, if you wanted to discuss it or had a
preprint/publication to send on. Again, any questions with teko please
doesn't hesitate to contact me.

Best,
Eric Cyr

From: Hoang Giang Bui <hgbk2008 at gmail.com<mailto:hgbk2008 at gmail.com>>
Date: Fri, 4 Jan 2013 12:01:51 +0100
To: Eric Cyr <eccyr at sandia.gov<mailto:eccyr at sandia.gov>>, "trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>" <trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>>
Subject: [EXTERNAL] Re: [Trilinos-Users] using teko


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<mailto: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<mailto:maherou at sandia.gov>>
Date: Thu, 3 Jan 2013 15:30:26 -0700
To: Eric Cyr <eccyr at sandia.gov<mailto:eccyr at sandia.gov>>, Hoang Giang Bui <hgbk2008 at gmail.com<mailto: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<mailto:hgbk2008 at gmail.com>>
Date: Tuesday, December 25, 2012 7:42 AM
To: "trilinos-users at software.sandia.gov<mailto:trilinos-users at software.sandia.gov>" <trilinos-users at software.sandia.gov<mailto: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/d01ece0d/attachment.html 


More information about the Trilinos-Users mailing list