441 const Teuchos::RCP<
const LinearOpSourceBase<double> > &fwdOpSrc,
442 const Teuchos::RCP<
const PreconditionerBase<double> > &prec,
443 const Teuchos::RCP<
const LinearOpSourceBase<double> > &approxFwdOpSrc,
444 const bool reusePrec,
445 LinearOpWithSolveBase<double> *Op
451 using Teuchos::rcp_dynamic_cast;
452 using Teuchos::rcp_const_cast;
453 using Teuchos::set_extra_data;
454 using Teuchos::get_optional_extra_data;
455 using Teuchos::get_optional_nonconst_extra_data;
456 using Teuchos::outArg;
457 typedef EpetraExt::ProductOperator PO;
459 const Teuchos::RCP<Teuchos::FancyOStream> out = this->getOStream();
460 const Teuchos::EVerbosityLevel verbLevel = this->getVerbLevel();
461 Teuchos::OSTab tab(out);
462 if(out.get() &&
static_cast<int>(verbLevel) >
static_cast<int>(Teuchos::VERB_LOW))
463 *out <<
"\nEntering Thyra::AztecOOLinearOpWithSolveFactory::initializeOp_impl(...) ...\n";
465 typedef Teuchos::VerboseObjectTempState<PreconditionerFactoryBase<double> > VOTSPF;
466 VOTSPF precFactoryOutputTempState(
precFactory_,out,verbLevel);
469 TEUCHOS_TEST_FOR_EXCEPT(Op==NULL);
470 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc.get()==NULL);
471 TEUCHOS_TEST_FOR_EXCEPT(fwdOpSrc->getOp().get()==NULL);
479 Teuchos::RCP<const LinearOpBase<double> >
480 tmpFwdOp = fwdOpSrc->getOp(),
481 tmpApproxFwdOp = ( approxFwdOpSrc.get() ? approxFwdOpSrc->getOp() : Teuchos::null );
482 Teuchos::RCP<const LinearOpBase<double> > fwdOp;
483 Teuchos::RCP<const LinearOpBase<double> > approxFwdOp;
484 if (
dynamic_cast<const EpetraLinearOpBase*
>(tmpFwdOp.get())!=0 )
487 approxFwdOp = tmpApproxFwdOp;
491 fwdOp = makeEpetraWrapper(tmpFwdOp);
495 dynamic_cast<const EpetraLinearOpBase*
>(&*tmpApproxFwdOp.get())
498 approxFwdOp = makeEpetraWrapper(tmpApproxFwdOp);
506 *aztecOp = &Teuchos::dyn_cast<AztecOOLinearOpWithSolve>(*Op);
511 Teuchos::RCP<const Epetra_Operator> epetra_epetraFwdOp;
512 EOpTransp epetra_epetraFwdOpTransp;
513 EApplyEpetraOpAs epetra_epetraFwdOpApplyAs;
514 EAdjointEpetraOp epetra_epetraFwdOpAdjointSupport;
515 double epetra_epetraFwdOpScalar;
516 epetraFwdOpViewExtractor_->getEpetraOpView(
518 outArg(epetra_epetraFwdOp), outArg(epetra_epetraFwdOpTransp),
519 outArg(epetra_epetraFwdOpApplyAs), outArg(epetra_epetraFwdOpAdjointSupport),
520 outArg(epetra_epetraFwdOpScalar)
522 TEUCHOS_TEST_FOR_EXCEPTION(
523 epetra_epetraFwdOp.get()==NULL, std::logic_error
524 ,
"Error, The input fwdOp object must be fully initialized "
525 "before calling this function!"
531 Teuchos::RCP<PreconditionerBase<double> > myPrec;
532 Teuchos::RCP<const PreconditionerBase<double> > precUsed;
541 ( !aztecOp->isExternalPrec()
542 ? Teuchos::rcp_const_cast<PreconditionerBase<double> >(
543 aztecOp->extract_prec())
560 RCP<const LinearOpBase<double> > rightPrecOp;
561 if (precUsed.get()) {
562 RCP<const LinearOpBase<double> > unspecified = precUsed->getUnspecifiedPrecOp();
563 RCP<const LinearOpBase<double> > left = precUsed->getLeftPrecOp();
564 RCP<const LinearOpBase<double> > right = precUsed->getRightPrecOp();
565 TEUCHOS_TEST_FOR_EXCEPTION(
566 !( left.get() || right.get() || unspecified.get() ), std::logic_error
567 ,
"Error, at least one preconditoner linear operator objects must be set!"
569 if(unspecified.get()) {
570 rightPrecOp = unspecified;
574 TEUCHOS_TEST_FOR_EXCEPTION(
575 left.get(),std::logic_error
576 ,
"Error, we can not currently handle a left"
577 " preconditioner with the AztecOO/Thyra adapters!"
582 double wrappedPrecOpScalar = 0.0;
583 EOpTransp wrappedPrecOpTransp = NOTRANS;
584 RCP<const LinearOpBase<double> > wrappedPrecOp = null;
585 RCP<const EpetraLinearOpBase> epetraPrecOp;
586 Teuchos::RCP<const Epetra_Operator> epetra_epetraPrecOp;
587 EOpTransp epetra_epetraPrecOpTransp;
588 EApplyEpetraOpAs epetra_epetraPrecOpApplyAs;
589 EAdjointEpetraOp epetra_epetraPrecOpAdjointSupport;
590 EOpTransp overall_epetra_epetraPrecOpTransp=NOTRANS;
591 if(rightPrecOp.get()) {
592 RCP<const LinearOpBase<double> > tmpWrappedPrecOp;
594 rightPrecOp,&wrappedPrecOpScalar,&wrappedPrecOpTransp,&tmpWrappedPrecOp);
595 if(
dynamic_cast<const EpetraLinearOpBase*
>(&*tmpWrappedPrecOp) ) {
596 wrappedPrecOp = tmpWrappedPrecOp;
599 wrappedPrecOp = makeEpetraWrapper(tmpWrappedPrecOp);
601 epetraPrecOp = rcp_dynamic_cast<const EpetraLinearOpBase>(
603 epetraPrecOp->getEpetraOpView(
604 outArg(epetra_epetraPrecOp), outArg(epetra_epetraPrecOpTransp),
605 outArg(epetra_epetraPrecOpApplyAs), outArg(epetra_epetraPrecOpAdjointSupport));
606 TEUCHOS_TEST_FOR_EXCEPTION(
607 epetra_epetraPrecOp.get()==NULL,std::logic_error
608 ,
"Error, The input prec object and its embedded preconditioner"
609 " operator must be fully initialized before calling this function!"
618 overall_epetra_epetraPrecOpTransp
620 real_trans(wrappedPrecOpTransp),
621 real_trans(epetra_epetraPrecOpTransp)
629 if(approxFwdOp.get()) {
632 unwrap(approxFwdOp,&wrappedPrecOpScalar,&wrappedPrecOpTransp,&wrappedPrecOp);
633 epetraPrecOp = rcp_dynamic_cast<const EpetraLinearOpBase>(
635 epetraPrecOp->getEpetraOpView(
636 outArg(epetra_epetraPrecOp), outArg(epetra_epetraPrecOpTransp),
637 outArg(epetra_epetraPrecOpApplyAs), outArg(epetra_epetraPrecOpAdjointSupport)
639 TEUCHOS_TEST_FOR_EXCEPTION(
640 epetra_epetraPrecOp.get()==NULL,std::logic_error
641 ,
"Error, The input approxFwdOp object must be fully initialized"
642 " before calling this function!"
652 overall_epetra_epetraPrecOpTransp
654 real_trans(wrappedPrecOpTransp),
655 real_trans(epetra_epetraPrecOpTransp)
663 RCP<const Epetra_RowMatrix>
664 rowmatrix_epetraFwdOp = rcp_dynamic_cast<const Epetra_RowMatrix>(
666 rowmatrix_epetraPrecOp = rcp_dynamic_cast<const Epetra_RowMatrix>(
667 epetra_epetraPrecOp);
673 enum ELocalPrecType {
674 PT_NONE, PT_AZTEC_FROM_OP, PT_AZTEC_FROM_APPROX_FWD_MATRIX,
675 PT_FROM_PREC_OP, PT_UPPER_BOUND
677 ELocalPrecType localPrecType = PT_UPPER_BOUND;
678 if( precUsed.get()==NULL && approxFwdOp.get()==NULL && !
useAztecPrec_ ) {
680 localPrecType = PT_NONE;
682 else if( precUsed.get()==NULL && approxFwdOp.get()==NULL &&
useAztecPrec_ ) {
685 localPrecType = PT_AZTEC_FROM_OP;
690 localPrecType = PT_AZTEC_FROM_APPROX_FWD_MATRIX;
692 else if( precUsed.get() ) {
695 localPrecType = PT_FROM_PREC_OP;
697 TEUCHOS_TEST_FOR_EXCEPTION
698 (localPrecType == PT_UPPER_BOUND, std::logic_error,
699 "AztecOOLinearOpWithSolveFactory::initializeOp_impl(...): "
700 "localPrecType == PT_UPPER_BOUND. This means that previously, "
701 "this value might have been used uninitialized. "
702 "Please report this bug to the Stratimikos developers.");
708 RCP<AztecOO> aztecFwdSolver, aztecAdjSolver;
714 Teuchos::RCP<const LinearOpBase<double> > old_fwdOp;
715 Teuchos::RCP<const LinearOpSourceBase<double> > old_fwdOpSrc;
716 Teuchos::RCP<const PreconditionerBase<double> > old_prec;
717 bool old_isExternalPrec;
718 Teuchos::RCP<const LinearOpSourceBase<double> > old_approxFwdOpSrc;
719 Teuchos::RCP<AztecOO> old_aztecFwdSolver;
720 Teuchos::RCP<AztecOO> old_aztecAdjSolver;
721 double old_aztecSolverScalar;
722 aztecOp->uninitialize(
732 ,&old_aztecSolverScalar
734 if( old_aztecFwdSolver.get()==NULL ) {
742 aztecFwdSolver = old_aztecFwdSolver;
743 aztecAdjSolver = old_aztecAdjSolver;
744 startingOver =
false;
747 Ptr<bool> constructedAztecPreconditioner;
751 !is_null(constructedAztecPreconditioner = get_optional_nonconst_extra_data<bool>(
752 aztecFwdSolver,
"AOOLOWSF::constructedAztecPreconditoner") )
754 *constructedAztecPreconditioner
757 aztecFwdSolver->DestroyPreconditioner();
758 *constructedAztecPreconditioner =
false;
762 Ptr<bool> setPreconditionerOperator;
764 localPrecType != PT_FROM_PREC_OP
765 && !is_null( setPreconditionerOperator = get_optional_nonconst_extra_data<bool>(
766 aztecFwdSolver,
"AOOLOWSF::setPreconditonerOperator") )
767 && *setPreconditionerOperator
783 aztecFwdSolver = rcp(
new AztecOO());
784 aztecFwdSolver->SetAztecOption(AZ_diagnostics,AZ_none);
785 aztecFwdSolver->SetAztecOption(AZ_keep_info,1);
788 epetra_epetraFwdOpAdjointSupport==EPETRA_OP_ADJOINT_SUPPORTED
790 localPrecType!=PT_AZTEC_FROM_OP && localPrecType!=PT_AZTEC_FROM_APPROX_FWD_MATRIX
793 aztecAdjSolver = rcp(
new AztecOO());
794 aztecAdjSolver->SetAztecOption(AZ_diagnostics,AZ_none);
805 &
paramList_->sublist(ForwardSolve_name).sublist(AztecOO_Settings_name),
810 &
paramList_->sublist(AdjointSolve_name).sublist(AztecOO_Settings_name),
818 RCP<const Epetra_Operator>
819 aztec_epetra_epetraFwdOp,
820 aztec_epetra_epetraAdjOp;
822 RCP<const Epetra_Operator>
824 = { epetra_epetraFwdOp };
827 = { epetra_epetraFwdOpTransp==NOTRANS ? Teuchos::NO_TRANS : Teuchos::TRANS };
830 = { epetra_epetraFwdOpApplyAs==EPETRA_OP_APPLY_APPLY
831 ? PO::APPLY_MODE_APPLY
832 : PO::APPLY_MODE_APPLY_INVERSE };
834 epetraOpsTransp[0] == Teuchos::NO_TRANS
836 epetraOpsApplyMode[0] == PO::APPLY_MODE_APPLY
839 aztec_epetra_epetraFwdOp = epetra_epetraFwdOp;
843 aztec_epetra_epetraFwdOp = rcp(
844 new PO(1,epetraOps,epetraOpsTransp,epetraOpsApplyMode));
849 aztec_epetra_epetraFwdOp.get() != aztecFwdSolver->GetUserOperator()
854 aztecFwdSolver->SetUserOperator(
855 const_cast<Epetra_Operator*
>(&*aztec_epetra_epetraFwdOp));
857 aztec_epetra_epetraFwdOp, AOOLOWSF_aztec_epetra_epetraFwdOp_str,
858 Teuchos::inOutArg(aztecFwdSolver), Teuchos::POST_DESTROY,
false
862 if( aztecAdjSolver.get() ) {
863 epetraOpsTransp[0] = (
864 epetra_epetraFwdOpTransp==NOTRANS
869 epetraOpsTransp[0] == Teuchos::NO_TRANS
871 epetraOpsApplyMode[0] == PO::APPLY_MODE_APPLY
874 aztec_epetra_epetraAdjOp = epetra_epetraFwdOp;
877 aztec_epetra_epetraAdjOp = rcp(
878 new PO(1,epetraOps,epetraOpsTransp,epetraOpsApplyMode));
880 aztecAdjSolver->SetUserOperator(
881 const_cast<Epetra_Operator*
>(&*aztec_epetra_epetraAdjOp));
883 aztec_epetra_epetraAdjOp, AOOLOWSF_aztec_epetra_epetraAdjOp_str,
884 Teuchos::inOutArg(aztecAdjSolver), Teuchos::POST_DESTROY,
false
891 RCP<const Epetra_Operator>
892 aztec_fwd_epetra_epetraPrecOp,
893 aztec_adj_epetra_epetraPrecOp;
894 bool setAztecPreconditioner =
false;
895 switch(localPrecType) {
902 case PT_AZTEC_FROM_OP: {
907 if( startingOver || !reusePrec ) {
908 TEUCHOS_TEST_FOR_EXCEPTION(
909 rowmatrix_epetraFwdOp.get()==NULL, std::logic_error,
910 "AztecOOLinearOpWithSolveFactory::initializeOp_impl(...): "
911 "Error, There is no preconditioner given by client, but the client "
912 "passed in an Epetra_Operator for the forward operator of type \'"
913 <<typeName(*epetra_epetraFwdOp)<<
"\' that does not "
914 "support the Epetra_RowMatrix interface!"
916 TEUCHOS_TEST_FOR_EXCEPTION(
917 epetra_epetraFwdOpTransp!=NOTRANS, std::logic_error,
918 "AztecOOLinearOpWithSolveFactory::initializeOp_impl(...):"
919 " Error, There is no preconditioner given by client and the client "
920 "passed in an Epetra_RowMatrix for the forward operator but the "
921 "overall transpose is not NOTRANS and therefore we can can just "
922 "hand this over to aztec without making a copy which is not supported here!"
924 aztecFwdSolver->SetPrecMatrix(
925 const_cast<Epetra_RowMatrix*
>(&*rowmatrix_epetraFwdOp));
927 rowmatrix_epetraFwdOp, AOOLOWSF_rowmatrix_epetraFwdOp_str,
928 Teuchos::inOutArg(aztecFwdSolver), Teuchos::POST_DESTROY,
false
931 setAztecPreconditioner =
true;
934 case PT_AZTEC_FROM_APPROX_FWD_MATRIX: {
939 if( startingOver || !reusePrec ) {
940 TEUCHOS_TEST_FOR_EXCEPTION(
941 rowmatrix_epetraPrecOp.get()==NULL, std::logic_error
942 ,
"AztecOOLinearOpWithSolveFactor::initializeOp_impl(...): The client "
943 "passed in an Epetra_Operator for the preconditioner matrix of type \'"
944 <<typeName(*epetra_epetraPrecOp)<<
"\' that does not "
945 "support the Epetra_RowMatrix interface!"
947 TEUCHOS_TEST_FOR_EXCEPTION(
948 overall_epetra_epetraPrecOpTransp!=NOTRANS, std::logic_error
949 ,
"AztecOOLinearOpWithSolveFactor::initializeOp_impl(...): Error, The client "
950 "passed in an Epetra_RowMatrix for the preconditoner matrix but the overall "
951 "transpose is not NOTRANS and therefore we can can just "
952 "hand this over to aztec without making a copy which is not supported here!"
954 aztecFwdSolver->SetPrecMatrix(
955 const_cast<Epetra_RowMatrix*
>(&*rowmatrix_epetraPrecOp));
957 rowmatrix_epetraPrecOp, AOOLOWSF_rowmatrix_epetraPrecOp_str,
958 Teuchos::inOutArg(aztecFwdSolver), Teuchos::POST_DESTROY,
false
961 setAztecPreconditioner =
true;
964 case PT_FROM_PREC_OP: {
969 RCP<const Epetra_Operator>
971 = { epetra_epetraPrecOp };
974 = { overall_epetra_epetraPrecOpTransp==NOTRANS
980 theEpetraOpsApplyMode[]
981 = { epetra_epetraPrecOpApplyAs==EPETRA_OP_APPLY_APPLY
982 ? PO::APPLY_MODE_APPLY_INVERSE
983 : PO::APPLY_MODE_APPLY };
985 theEpetraOpsTransp[0] == Teuchos::NO_TRANS
987 epetra_epetraPrecOpApplyAs==EPETRA_OP_APPLY_APPLY_INVERSE
990 aztec_fwd_epetra_epetraPrecOp = epetra_epetraPrecOp;
993 aztec_fwd_epetra_epetraPrecOp = rcp(
new PO(1,theEpetraOps,theEpetraOpsTransp,theEpetraOpsApplyMode));
995 aztecFwdSolver->SetPrecOperator(
996 const_cast<Epetra_Operator*
>(&*aztec_fwd_epetra_epetraPrecOp));
998 aztec_fwd_epetra_epetraPrecOp, AOOLOWSF_aztec_fwd_epetra_epetraPrecOp_str,
999 Teuchos::inOutArg(aztecFwdSolver), Teuchos::POST_DESTROY,
false
1003 aztecAdjSolver.get()
1005 epetra_epetraPrecOpAdjointSupport == EPETRA_OP_ADJOINT_SUPPORTED
1008 theEpetraOpsTransp[0] = (
1009 overall_epetra_epetraPrecOpTransp==NOTRANS
1014 theEpetraOpsTransp[0] == Teuchos::NO_TRANS
1016 epetra_epetraPrecOpApplyAs==EPETRA_OP_APPLY_APPLY_INVERSE
1019 aztec_adj_epetra_epetraPrecOp = epetra_epetraPrecOp;
1022 aztec_adj_epetra_epetraPrecOp = rcp(
1023 new PO(1,theEpetraOps,theEpetraOpsTransp,theEpetraOpsApplyMode));
1025 aztecAdjSolver->SetPrecOperator(
1026 const_cast<Epetra_Operator*
>(&*aztec_adj_epetra_epetraPrecOp));
1028 aztec_adj_epetra_epetraPrecOp, AOOLOWSF_aztec_adj_epetra_epetraPrecOp_str,
1029 Teuchos::inOutArg(aztecAdjSolver), Teuchos::POST_DESTROY,
false
1031 set_extra_data<bool>(
1032 true, AOOLOWSF_setPrecondtionerOperator_str,
1033 Teuchos::inOutArg(aztecFwdSolver), Teuchos::POST_DESTROY,
false
1039 TEUCHOS_TEST_FOR_EXCEPT(
true);
1045 if(setAztecPreconditioner) {
1046 if( startingOver || !reusePrec ) {
1047 double condNumEst = -1.0;
1048 TEUCHOS_TEST_FOR_EXCEPT(0!=aztecFwdSolver->ConstructPreconditioner(condNumEst));
1050 set_extra_data<bool>(
1051 true, AOOLOWSF_constructedAztecPreconditoner_str,
1052 Teuchos::inOutArg(aztecFwdSolver), Teuchos::POST_DESTROY,
false
1063 if(aztecAdjSolver.get() && aztecAdjSolver->GetPrecOperator()) {
1064 aztecOp->initialize(
1065 fwdOp, fwdOpSrc,precUsed, prec.get()!=NULL, approxFwdOpSrc,
1066 aztecFwdSolver,
true, aztecAdjSolver,
true, epetra_epetraFwdOpScalar
1070 aztecOp->initialize(
1071 fwdOp, fwdOpSrc, precUsed, prec.get()!=NULL, approxFwdOpSrc,
1072 aztecFwdSolver,
true, null,
false, epetra_epetraFwdOpScalar
1080 aztecOp->setOStream(this->getOStream());
1081 if(!is_null(this->getOverridingOStream()))
1082 aztecOp->setOverridingOStream(this->getOverridingOStream());
1083 aztecOp->setVerbLevel(this->getVerbLevel());
1090 if(out.get() &&
static_cast<int>(verbLevel) >
static_cast<int>(Teuchos::VERB_LOW))
1091 *out <<
"\nLeaving Thyra::AztecOOLinearOpWithSolveFactory::initializeOp_impl(...) ...\n";