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

[Hoang Giang Bui]

Beautiful!

Thanks for your help

BR..Bui


On Sat, Jan 5, 2013 at 11:48 PM, Bartlett, Roscoe A. <bartlettra at ornl.gov>wrote:

> You likely need to add the include:****
>
> ** **
>
> #include "Teuchos_RCPBoostSharedPtrConversions.hpp"****
>
> ** **
>
> -Ross****
>
> ** **
>
> ** **
>
> ** **
>
> *From:* Hoang Giang Bui [mailto:hgbk2008 at gmail.com]
> *Sent:* Friday, January 04, 2013 9:08 PM
> *To:* Bartlett, Roscoe A.; trilinos-users at software.sandia.gov
>
> *Subject:* Re: [Trilinos-Users] [EXTERNAL] Re: using teko****
>
> ** **
>
> Sounds good but it doesn't compile again. The configure output contains:
> ...
> -- Processing enabled TPL: Boost
> --   TPL_Boost_INCLUDE_DIRS='/opt/boost/include'
> ...
> -- Setting TeuchosCore_ENABLE_Boost=ON since TPL_ENABLE_Boost=ON
> ...****
>
> Is that Boost enabled as you mention?****
>
> BR
> Giang Bui****
>
> ** **
>
> On Fri, Jan 4, 2013 at 6:52 PM, Bartlett, Roscoe A. <bartlettra at ornl.gov>
> wrote:****
>
> The boost::shared_ptr is handled correctly when passed to the RCP. Just
> enable support for boost.
>
> -Ross
>
> Sent from my android phone. Excuse the terse message.****
>
>
>
>
> -----Original Message-----
> *From: *Cyr, Eric C [eccyr at sandia.gov]
> *Sent: *Friday, January 04, 2013 12:22 PM Eastern Standard Time
> *To: *Hoang Giang Bui; trilinos-users at software.sandia.gov
> *Subject: *Re: [Trilinos-Users] [EXTERNAL] Re: using teko****
>
> 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>
> *Date:* Fri, 4 Jan 2013 12:01:51 +0100
> *To:* Eric Cyr <eccyr at sandia.gov>, "trilinos-users at software.sandia.gov" <
> 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> 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****
>
>
>
>
> --
> With Best Regards !
> Giang Bui****
>
> To learn and to excel****
>



-- 
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/20130106/448f8b9c/attachment-0001.html