C++ Thyra::epetraLinearOp方法代码示例

void tLSCIntegrationTest::loadStableSystem()
{
   Epetra_CrsMatrix *F=0, *B=0, *Bt=0,*Qu=0;

   F=0; B=0; Bt=0; Qu=0;

   // read in stable discretization
   TEST_FOR_EXCEPT(EpetraExt::MatrixMarketFileToCrsMatrix("./data/lsc_F_2.mm",*velMap_,*velMap_,*velMap_,F));
   TEST_FOR_EXCEPT(EpetraExt::MatrixMarketFileToCrsMatrix("./data/lsc_B_2.mm",*prsMap_,*prsMap_,*velMap_,B));
   TEST_FOR_EXCEPT(EpetraExt::MatrixMarketFileToCrsMatrix("./data/lsc_Bt_2.mm",*velMap_,*velMap_,*prsMap_,Bt));
   TEST_FOR_EXCEPT(EpetraExt::MatrixMarketFileToCrsMatrix("./data/lsc_Qu_2.mm",*velMap_,*velMap_,*velMap_,Qu));

   // set stable matrix pointers
   sF_  = rcp(F); sB_  = rcp(B); sBt_ = rcp(Bt); sQu_ = rcp(Qu);

   Teko::LinearOp C;
   Teko::LinearOp tA_ = Teko::block2x2<double>(epetraLinearOp(sF_),epetraLinearOp(sBt_),epetraLinearOp(sB_),C,"A");
   sA_ = rcp(new Teko::Epetra::EpetraOperatorWrapper(tA_));

   // build an exporter to work around issue with MMFileToVector
   Epetra_Export exporter(*fullMap_,sA_->OperatorRangeMap());

   // read in RHS vector
   {
      Epetra_Vector *vfull=0, *temp=0;
 
      // read in rhs file 
      TEST_FOR_EXCEPT(EpetraExt::MatrixMarketFileToVector("./data/lsc_rhs.mm",*fullMap_,vfull));

      // MMFileToVector is immplemented incompletely...thats why an exporter is used
      temp = new Epetra_Vector(sA_->OperatorRangeMap());
      temp->Export(*vfull,exporter,Insert);
      rhs_ = rcp(temp);

      delete vfull;
   }

   // read in solution vector
   {
      Epetra_Vector *vfull=0, *temp=0;
 
      // read in exact solution file 
      TEST_FOR_EXCEPT(EpetraExt::MatrixMarketFileToVector("./data/lsc_exact_2.mm",*fullMap_,vfull));

      // MMFileToVector is immplemented incompletely...thats why an exporter is used
      temp = new Epetra_Vector(sA_->OperatorRangeMap());
      temp->Export(*vfull,exporter,Insert);
      sExact_ = rcp(temp);

      delete vfull;
   }
}

bool tLSCIntegrationTest::test_withmassStable(int verbosity,std::ostream & os)
{
   Teuchos::ParameterList paramList;
   solveList(paramList,8);

   RCP<Teko::InverseFactory> invFact = Teko::invFactoryFromParamList(paramList,"ML");
   TEUCHOS_ASSERT(invFact!=Teuchos::null);

   bool status = false;
   bool allPassed = true;

   // load everything
   loadStableSystem();

   // if you get here you automatically pass the first test
   if(verbosity>=10 ) {
      os << std::endl << "   tLSCIntegrationTest::test_withmassStable: loading system ... " 
         << toString(true) << std::endl;
   }

   LinearOp Qu = epetraLinearOp(sQu_);
   const RCP<Teko::NS::LSCStrategy> strategy = rcp(new Teko::NS::InvLSCStrategy(invFact,Qu));
   const RCP<Teko::BlockPreconditionerFactory> precFact = rcp(new Teko::NS::LSCPreconditionerFactory(strategy));
   const RCP<Teko::Epetra::EpetraBlockPreconditioner> prec = rcp(new Teko::Epetra::EpetraBlockPreconditioner(precFact));
   prec->buildPreconditioner(sA_);

   // B. Build solver and solve system
   Epetra_Vector x(sA_->OperatorDomainMap());
   x.Scale(0.0);

   // build Epetra problem
   Epetra_LinearProblem problem(&*sA_,&x,&*rhs_); // this doesn't take const arguments!

   AztecOO solver(problem);
   solver.SetAztecOption(AZ_solver,AZ_gmres);
   solver.SetAztecOption(AZ_precond,AZ_none);
   solver.SetAztecOption(AZ_kspace,50);
   solver.SetAztecOption(AZ_output,AZ_none);
   solver.SetPrecOperator(&*prec);

   solver.Iterate(1000,1e-8);

   // check iteration count
   status = (solver.NumIters()<=16);
   if(not status || verbosity>=10 ) { 
      os << std::endl << "   tLSCIntegrationTest::test_withmassStable " << toString(status) 
                      << ": # of iterations = " << solver.NumIters() << " (should be 16)" << std::endl;
   }
   allPassed &= status;
 
   // check exact answer (versus IFISS solution)
   x.Update(-1.0,*sExact_,1.0); // x = x - x*
   double errnorm,exactnorm,relerr;
   x.Norm2(&errnorm);
   sExact_->Norm2(&exactnorm);
   status = ((relerr = errnorm/exactnorm) <= tolerance_);
   if(not status || verbosity>=10 ) { 
      os << std::endl << "   tLSCIntegrationTest::test_withmassStable " << toString(status) 
                      << ": error in solution = " << std::scientific << relerr << " <= " << tolerance_ << std::endl;
   }
   allPassed &= status;

   return allPassed;
}