C++ RCP::Compute方法代码示例

LinearOp DiagonalPreconditionerFactory::buildPreconditionerOperator(LinearOp & lo,PreconditionerState & state) const
{
  if(diagonalType_==BlkDiag) {
     // Sanity check the state
     DiagonalPrecondState & MyState = Teuchos::dyn_cast<DiagonalPrecondState>(state);
   
     // Get the underlying Epetra_CrsMatrix, if we have one
     Teuchos::RCP<const Epetra_Operator> eo=Thyra::get_Epetra_Operator(*lo);
     TEUCHOS_ASSERT(eo!=Teuchos::null);
     Teuchos::RCP<const Epetra_CrsMatrix> MAT = Teuchos::rcp_dynamic_cast<const Epetra_CrsMatrix>(eo);
     TEUCHOS_ASSERT(MAT!=Teuchos::null);
   
     // Create a new EpetraExt_PointToBlockDiagPermute for the state object, if we don't have one
     Teuchos::RCP<EpetraExt_PointToBlockDiagPermute> BDP;
     if(MyState.BDP_==Teuchos::null) {
       BDP = Teuchos::rcp(new EpetraExt_PointToBlockDiagPermute(*MAT));
       BDP->SetParameters(List_);
       BDP->Compute();
       MyState.BDP_ = BDP;
     }

     RCP<Epetra_FECrsMatrix> Hcrs=rcp(MyState.BDP_->CreateFECrsMatrix());
     return Thyra::epetraLinearOp(Hcrs);
   
     // Build the LinearOp object  (NTS: swapping the range and domain)
     // LinearOp MyOp = Teuchos::rcp(new DiagonalPreconditionerOp(MyState.BDP_,lo->domain(),lo->range()));
  }

  return getInvDiagonalOp(lo,diagonalType_);
}

//.........这里部分代码省略.........
    Teuchos::RCP<Epetra_Vector> solution_vec =
      interface->getSolution();
    Teuchos::RCP<Epetra_Vector> rhs_vec =
      Teuchos::rcp(new Epetra_Vector(*solution_vec));
    Teuchos::RCP<Epetra_Vector> lhs_vec =
      Teuchos::rcp(new Epetra_Vector(*solution_vec));
    Teuchos::RCP<Epetra_CrsMatrix> jacobian_matrix =
      interface->getJacobian();


    if (verbose)
      p.out() << "Evaluating F and J" << std::endl;
    solution_vec->PutScalar(1.0);
    interface->computeF(*solution_vec, *rhs_vec);
    rhs_vec->Scale(-1.0);
    interface->computeJacobian(*solution_vec, *jacobian_matrix);

    double norm =0.0;
    rhs_vec->Norm2(&norm);
    if (verbose)
      p.out() << "Step 0, ||F|| = " << norm << std::endl;


    if (verbose)
      p.out() << "Creating Ifpack preconditioner" << std::endl;

    Ifpack Factory;
    Teuchos::RCP<Ifpack_Preconditioner> PreconditionerPtr;
    PreconditionerPtr = Teuchos::rcp(Factory.Create("ILU",
          jacobian_matrix.get(),0));
    Teuchos::ParameterList teuchosParams;
    PreconditionerPtr->SetParameters(teuchosParams);
    PreconditionerPtr->Initialize();
    PreconditionerPtr->Compute();


    if (verbose)
      p.out() << "Creating Aztec Solver" << std::endl;

    Teuchos::RCP<AztecOO> aztecSolverPtr = Teuchos::rcp(new AztecOO());
    if (verbose)
      aztecSolverPtr->SetAztecOption(AZ_output, AZ_last);
    else
      aztecSolverPtr->SetAztecOption(AZ_output, AZ_none);

    // *******************************
    // Reuse Test
    // *******************************

    if (verbose) {
      p.out() << "**********************************************" << std::endl;
      p.out() << "Testing Newton solve with prec reuse" << std::endl;
      p.out() << "**********************************************" << std::endl;
    }

    int step_number = 0;
    int max_steps = 20;
    bool converged = false;
    int total_linear_iterations = 0;
    while (norm > 1.0e-8 && step_number < max_steps) {

      step_number++;

      if (verbose)
        p.out() << "Step " << step_number << ", ||F|| = " << norm << std::endl;

int main(int argc, char *argv[]) {
  using std::cout;
  using std::endl;
  int i;

#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  int MyPID = Comm.MyPID();

  // Number of dimension of the domain
  int space_dim = 2;

  // Size of each of the dimensions of the domain
  std::vector<double> brick_dim( space_dim );
  brick_dim[0] = 1.0;
  brick_dim[1] = 1.0;

  // Number of elements in each of the dimensions of the domain
  std::vector<int> elements( space_dim );
  elements[0] = 10;
  elements[1] = 10;

  // Create problem
  Teuchos::RCP<ModalProblem> testCase = Teuchos::rcp( new ModeLaplace2DQ2(Comm, brick_dim[0], elements[0], brick_dim[1], elements[1]) );

  // Get the stiffness and mass matrices
  Teuchos::RCP<Epetra_CrsMatrix> K = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getStiffness()), false );
  Teuchos::RCP<Epetra_CrsMatrix> M = Teuchos::rcp( const_cast<Epetra_CrsMatrix *>(testCase->getMass()), false );

  //
  // ************Construct preconditioner*************
  //
  Teuchos::ParameterList ifpackList;

  // allocates an IFPACK factory. No data is associated
  // to this object (only method Create()).
  Ifpack Factory;

  // create the preconditioner. For valid PrecType values,
  // please check the documentation
  std::string PrecType = "ICT"; // incomplete Cholesky
  int OverlapLevel = 0; // must be >= 0. If Comm.NumProc() == 1,
                        // it is ignored.

  Teuchos::RCP<Ifpack_Preconditioner> Prec = Teuchos::rcp( Factory.Create(PrecType, &*K, OverlapLevel) );
  assert(Prec != Teuchos::null);

  // specify parameters for ICT
  ifpackList.set("fact: drop tolerance", 1e-4);
  ifpackList.set("fact: ict level-of-fill", 0.);
  // the combine mode is on the following:
  // "Add", "Zero", "Insert", "InsertAdd", "Average", "AbsMax"
  // Their meaning is as defined in file Epetra_CombineMode.h
  ifpackList.set("schwarz: combine mode", "Add");
  // sets the parameters
  IFPACK_CHK_ERR(Prec->SetParameters(ifpackList));

  // initialize the preconditioner. At this point the matrix must
  // have been FillComplete()'d, but actual values are ignored.
  IFPACK_CHK_ERR(Prec->Initialize());

  // Builds the preconditioners, by looking for the values of
  // the matrix.
  IFPACK_CHK_ERR(Prec->Compute());

  //
  //*******************************************************/
  // Set up Belos Block CG operator for inner iteration
  //*******************************************************/
  //
  int blockSize = 3; // block size used by linear solver and eigensolver [ not required to be the same ]
  int maxits = K->NumGlobalRows(); // maximum number of iterations to run
  //
  // Create the Belos::LinearProblem
  //
  Teuchos::RCP<Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator> >
    My_LP = Teuchos::rcp( new Belos::LinearProblem<double,Epetra_MultiVector,Epetra_Operator>() );
  My_LP->setOperator( K );

  // Create the Belos preconditioned operator from the Ifpack preconditioner.
  // NOTE:  This is necessary because Belos expects an operator to apply the
  //        preconditioner with Apply() NOT ApplyInverse().
  Teuchos::RCP<Epetra_Operator> belosPrec = Teuchos::rcp( new Epetra_InvOperator( Prec.get() ) );
  My_LP->setLeftPrec( belosPrec );
  //
  // Create the ParameterList for the Belos Operator
  //
  Teuchos::RCP<Teuchos::ParameterList> My_List = Teuchos::rcp( new Teuchos::ParameterList() );
  My_List->set( "Solver", "BlockCG" );
  My_List->set( "Maximum Iterations", maxits );
  My_List->set( "Block Size", 1 );
  My_List->set( "Convergence Tolerance", 1e-12 );
  //
  // Create the Belos::EpetraOperator
//.........这里部分代码省略.........