C++ epetra::Vector类代码示例

NOX::Epetra::LinearSystemAmesos::
LinearSystemAmesos(
  Teuchos::ParameterList& printingParams, 
  Teuchos::ParameterList& linearSolverParams, 
  const Teuchos::RCP<NOX::Epetra::Interface::Required>& iReq, 
  const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac, 
  const Teuchos::RCP<Epetra_Operator>& J,
  const NOX::Epetra::Vector& cloneVector,
  const Teuchos::RCP<NOX::Epetra::Scaling> s):
  amesosProblem(Teuchos::null),
  amesosSolver(Teuchos::null),
  factory(),
  isValidFactorization(false),
  jacInterfacePtr(iJac),
  jacPtr(J),
  leftHandSide(Teuchos::rcp(new Epetra_Vector(cloneVector.getEpetraVector()))),
  rightHandSide(Teuchos::rcp(new Epetra_Vector(cloneVector.getEpetraVector()))),
  scaling(s),
  timer(cloneVector.getEpetraVector().Comm()),
  utils(printingParams)
{
  amesosProblem = Teuchos::rcp(new Epetra_LinearProblem(
				      dynamic_cast<Epetra_CrsMatrix *>(jacPtr.get()),
				      leftHandSide.get(),
				      rightHandSide.get()));

  Amesos_BaseSolver * tmp = factory.Create(linearSolverParams.get("Amesos Solver","Amesos_Klu"), 
      *amesosProblem);
  TEUCHOS_TEST_FOR_EXCEPTION ( tmp == 0, Teuchos::Exceptions::InvalidParameterValue, 
      "Invalid Amesos Solver: " << linearSolverParams.get<string>("Amesos Solver"));
  amesosSolver = Teuchos::rcp(tmp);

  amesosSolver->SetParameters(linearSolverParams);
}

//***********************************************************************
bool NOX::Epetra::LinearSystemStratimikos::
applyJacobian(const NOX::Epetra::Vector& input, 
	      NOX::Epetra::Vector& result) const
{
  jacPtr->SetUseTranspose(false);
  int status = jacPtr->Apply(input.getEpetraVector(), 
  		  	     result.getEpetraVector());
  return (status == 0);
}

//***********************************************************************
NOX::Epetra::LinearSystemStratimikos::
LinearSystemStratimikos(
 Teuchos::ParameterList& printParams, 
 Teuchos::ParameterList& stratSolverParams,  
 const Teuchos::RCP<NOX::Epetra::Interface::Jacobian>& iJac, 
 const Teuchos::RCP<Epetra_Operator>& jacobian,
 const NOX::Epetra::Vector& cloneVector,
 const Teuchos::RCP<NOX::Epetra::Scaling> s):
  utils(printParams),
  jacInterfacePtr(iJac),
  jacType(EpetraOperator),
  jacPtr(jacobian),
  precType(EpetraOperator),
  precMatrixSource(UseJacobian),
  scaling(s),
  conditionNumberEstimate(0.0),
  isPrecConstructed(false),
  precQueryCounter(0),
  maxAgeOfPrec(1),
  timer(cloneVector.getEpetraVector().Comm()),
  timeCreatePreconditioner(0.0),
  timeApplyJacbianInverse(0.0),
  getLinearSolveToleranceFromNox(false)
{
  jacType = getOperatorType(*jacPtr);
  precType = jacType;

  reset(stratSolverParams.sublist("NOX Stratimikos Options"));

  // Allocate solver
  initializeStratimikos(stratSolverParams.sublist("Stratimikos"));
  tmpVectorPtr = Teuchos::rcp(new NOX::Epetra::Vector(cloneVector));
}

//***********************************************************************
bool NOX::Epetra::LinearSystemStratimikos::
computeJacobian(const NOX::Epetra::Vector& x)
{
  bool success = jacInterfacePtr->computeJacobian(x.getEpetraVector(), 
						  *jacPtr);
  return success;
}

bool
LOCA::Epetra::TransposeLinearSystem::LeftPreconditioning::
applyJacobianTransposeInverse(Teuchos::ParameterList &params,
                  const NOX::Epetra::Vector &input,
                  NOX::Epetra::Vector &result)
{
  // Create preconditioned operator
  Teuchos::RCP<Epetra_Operator> left_prec_jac =
    Teuchos::rcp(new LOCA::Epetra::LeftPreconditionedOp(jac, prec));

  // Replace Jacobian operator with transposed, left preconditioned op
  linsys->setJacobianOperatorForSolve(left_prec_jac);

  // Create identity operator as a right preconditioner
  Teuchos::RCP<Epetra_Operator> identity_prec =
    Teuchos::rcp(new LOCA::Epetra::IdentityOp(
               Teuchos::rcp(&(jac->Comm()), false),
               Teuchos::rcp(&(jac->OperatorDomainMap()), false)));

  // Replace preconditioner with identity
  linsys->setPrecOperatorForSolve(identity_prec);

  // Precondition the RHS
  NOX::Epetra::Vector prec_input(input);
  prec->ApplyInverse(input.getEpetraVector(), prec_input.getEpetraVector());

  // Solve the system
  bool res = linsys->applyJacobianInverse(params, prec_input, result);

  return res;
}

bool NOX::Epetra::LinearSystemMPBD::
applyJacobianInverse(Teuchos::ParameterList &params, 
		     const NOX::Epetra::Vector &input, 
		     NOX::Epetra::Vector &result)
{
  TEUCHOS_FUNC_TIME_MONITOR("Total deterministic solve Time");

  // Extract blocks
  EpetraExt::BlockVector input_block(View, *base_map, 
				     input.getEpetraVector());
  EpetraExt::BlockVector result_block(View, *base_map, 
				      result.getEpetraVector());
  result_block.PutScalar(0.0);
   
  
  Teuchos::ParameterList& block_solver_params = 
    params.sublist("Deterministic Solver Parameters");
  
  // Solve block linear systems
  bool final_status = true;
  bool status;
  for (int i=0; i<num_mp_blocks; i++) {
    NOX::Epetra::Vector nox_input(input_block.GetBlock(i), 
				  NOX::Epetra::Vector::CreateView);
    NOX::Epetra::Vector nox_result(result_block.GetBlock(i), 
				   NOX::Epetra::Vector::CreateView);
    
    block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));

    if (precStrategy == STANDARD)
      block_solver->setPrecOperatorForSolve(precs[i]);
    else if (precStrategy == ON_THE_FLY) {
      block_solver->createPreconditioner(*(prec_x->GetBlock(i)), 
					 block_solver_params, false);
    }

    status = block_solver->applyJacobianInverse(block_solver_params, nox_input, 
						nox_result);
    final_status = final_status && status;
  }

  return final_status;
}

NOX::Epetra::Vector::Vector(const NOX::Epetra::Vector& source,
                NOX::CopyType type)
{
  vectorSpace = source.vectorSpace;

  switch (type) {

  case DeepCopy:        // default behavior

    epetraVec = Teuchos::rcp(new Epetra_Vector(source.getEpetraVector()));
    break;

  case ShapeCopy:

    epetraVec =
      Teuchos::rcp(new Epetra_Vector(source.getEpetraVector().Map()));
    break;

  }
}

bool NOX::Epetra::LinearSystemMPBD::
recomputePreconditioner(const NOX::Epetra::Vector& x, 
			Teuchos::ParameterList& p) const
{
  EpetraExt::BlockVector mp_x_block(View, *base_map, x.getEpetraVector());
  Teuchos::ParameterList& solverParams = 
    p.sublist("Deterministic Solver Parameters");
  bool total_success = true;

  if (precStrategy == STANDARD) {
    for (int i=0; i<num_mp_blocks; i++) {
      block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(i));
      if (precs[i] != Teuchos::null)
	block_solver->setPrecOperatorForSolve(precs[i]);
      bool success = 
	block_solver->recomputePreconditioner(*(mp_x_block.GetBlock(i)), 
					      solverParams);
      precs[i] = block_solver->getGeneratedPrecOperator();
      total_success = total_success && success;
    }
  }

  else if (precStrategy == MEAN) {
    block_solver->setJacobianOperatorForSolve(block_ops->getCoeffPtr(0));
    bool success = 
      block_solver->recomputePreconditioner(*(mp_x_block.GetBlock(0)), 
					    solverParams);
    total_success = total_success && success;
  }

  else if (precStrategy == ON_THE_FLY) {
    if (prec_x == Teuchos::null)
      prec_x = Teuchos::rcp(new EpetraExt::BlockVector(mp_x_block));
    else
      *prec_x = mp_x_block;
  }

  return total_success;
}

//***********************************************************************
bool NOX::Epetra::LinearSystemStratimikos::
createPreconditioner(const NOX::Epetra::Vector& x, Teuchos::ParameterList& p, 
		     bool recomputeGraph) const
{
  using Teuchos::RCP;
  using Teuchos::rcp;

  NOX_FUNC_TIME_MONITOR("NOX: Total Preconditioner Generation Time");

  double startTime = timer.WallTime();  

  if (utils.isPrintType(Utils::LinearSolverDetails))
    utils.out() << "\n       Creating a new preconditioner" << std::endl;;


//  Teuchos::RCP<Thyra::PreconditionerBase<double> > precObj;

  if (precMatrixSource == UseJacobian) {
    RCP<Thyra::PreconditionerFactoryBase<double> > precFactory =
      lowsFactory->getPreconditionerFactory();

    if (precFactory != Teuchos::null) {
      precObj = precFactory->createPrec();

      // Wrap Thyra objects around Epetra Op
      RCP<const Thyra::LinearOpBase<double> > linearOp =
        Thyra::epetraLinearOp(jacPtr);
  
      RCP<const Thyra::LinearOpSourceBase<double> > losb =
           rcp(new Thyra::DefaultLinearOpSource<double>(linearOp));

      // Computation of prec (e.g. ilu) happens here:
      precFactory->initializePrec(losb, precObj.get());

      // Get underlying Epetra operator
      Teuchos::RCP<Thyra::LinearOpBase<double> > pop;
      pop = precObj->getNonconstRightPrecOp();
      if (pop == Teuchos::null)
      	pop = precObj->getNonconstUnspecifiedPrecOp();
      solvePrecOpPtr = 
      	Teuchos::rcp_dynamic_cast<Thyra::EpetraLinearOp>(pop,true)->epetra_op();
    }
    else // no preconditioner
      precObj = Teuchos::null;

  }
  else if (precMatrixSource == SeparateMatrix) {

    // Compute matrix for use in preconditioning
    precInterfacePtr->computePreconditioner(x.getEpetraVector(), 
				      *precPtr, &p);

    RCP<Thyra::PreconditionerFactoryBase<double> > precFactory =
      lowsFactory->getPreconditionerFactory();

    if (precFactory != Teuchos::null) {
      precObj = precFactory->createPrec();

      // Send Prec Matrix to 
      RCP<const Thyra::LinearOpBase<double> > precOp =
        Thyra::epetraLinearOp(precPtr);

      RCP<const Thyra::LinearOpSourceBase<double> > losb =
           rcp(new Thyra::DefaultLinearOpSource<double>(precOp));

      // Computation of prec (e.g. ilu) happens here:
      precFactory->initializePrec(losb, precObj.get());

      // Get underlying Epetra operator
      Teuchos::RCP<Thyra::LinearOpBase<double> > pop;
      pop = precObj->getNonconstRightPrecOp();
      if (pop == Teuchos::null)
      	pop = precObj->getNonconstUnspecifiedPrecOp();
      solvePrecOpPtr = 
      	Teuchos::rcp_dynamic_cast<Thyra::EpetraLinearOp>(pop,true)->epetra_op();
    }
    else // no preconditioner
      precObj = Teuchos::null;

  }
  else if (precMatrixSource == UserDefined_) {

    precInterfacePtr->computePreconditioner(x.getEpetraVector(),
					    *precPtr, &p);

    // Wrap the preconditioner so that apply() calls ApplyInverse()
    RCP<const Thyra::LinearOpBase<double> > precOp =
      Thyra::epetraLinearOp(precPtr, 
			    Thyra::NOTRANS, 
			    Thyra::EPETRA_OP_APPLY_APPLY_INVERSE);

    RCP<Thyra::DefaultPreconditioner<double> > precObjDef =
       rcp(new Thyra::DefaultPreconditioner<double>);
    precObjDef->initializeRight(precOp);
    precObj = precObjDef;
    solvePrecOpPtr = precPtr;
  }

  isPrecConstructed = true; 
//.........这里部分代码省略.........

//***********************************************************************
bool NOX::Epetra::LinearSystemStratimikos::
applyJacobianInverse(Teuchos::ParameterList &p,
		     const NOX::Epetra::Vector& input, 
		     NOX::Epetra::Vector& result)
{
  using Teuchos::RCP;
  using Teuchos::rcp;

  NOX_FUNC_TIME_MONITOR("NOX: Total Linear Solve Time");

  double startTime = timer.WallTime();

  // Need non-const version of the input vector
  // Epetra_LinearProblem requires non-const versions so we can perform
  // scaling of the linear problem.
  NOX::Epetra::Vector& nonConstInput = const_cast<NOX::Epetra::Vector&>(input);
  
  // Zero out the delta X of the linear problem if requested by user.
  if (zeroInitialGuess) result.init(0.0);

  // Wrap Thyra objects around Epetra and NOX objects
  Teuchos::RCP<const Thyra::LinearOpBase<double> > linearOp =
    Thyra::epetraLinearOp(jacPtr);

  // Set the linear Op and  precomputed prec on this lows
  if (precObj == Teuchos::null) 
    Thyra::initializeOp(*lowsFactory, linearOp, lows.ptr());
  else 
    Thyra::initializePreconditionedOp<double>(
      *lowsFactory, linearOp, precObj, lows.ptr());

  Teuchos::RCP<Epetra_Vector> resultRCP =
    Teuchos::rcp(&result.getEpetraVector(), false);
  Teuchos::RCP<Epetra_Vector> inputRCP =
    Teuchos::rcp(&nonConstInput.getEpetraVector(), false);

  Teuchos::RCP<Thyra::VectorBase<double> >
    x = Thyra::create_Vector(resultRCP , linearOp->domain() );
  Teuchos::RCP<const Thyra::VectorBase<double> >
    b = Thyra::create_Vector(inputRCP, linearOp->range() );

  // Alter the convergence tolerance, if Inexact Newton
  Teuchos::RCP<Thyra::SolveCriteria<double> > solveCriteria;
  if (getLinearSolveToleranceFromNox) {
    Thyra::SolveMeasureType solveMeasure(
        Thyra::SOLVE_MEASURE_NORM_RESIDUAL,
        Thyra::SOLVE_MEASURE_NORM_INIT_RESIDUAL );
    solveCriteria = Teuchos::rcp(new Thyra::SolveCriteria<double>(
        solveMeasure, p.get<double>("Tolerance") ) );
  }

  // Solve the linear system for x
  Thyra::SolveStatus<double> status =
    lows->solve(Thyra::NOTRANS, *b, x.ptr(), solveCriteria.ptr());

  // MOVE TO FUNCTION: Update statistics: solves, iters, iters_total, achieved tol
  ++linearSolveCount;
  if (status.extraParameters != Teuchos::null) {
    if (status.extraParameters->isParameter("Belos/Iteration Count")) {
      linearSolveIters_last = status.extraParameters->get<int>("Belos/Iteration Count");
      linearSolveIters_total += linearSolveIters_last;
    }
    if (status.extraParameters->isParameter("Belos/Achieved Tolerance")) 
      linearSolveAchievedTol = status.extraParameters->get<double>("Belos/Achieved Tolerance");
    if (status.extraParameters->isParameter("AztecOO/Iteration Count")) {

      linearSolveIters_last = status.extraParameters->get<int>("AztecOO/Iteration Count");
      linearSolveIters_total += linearSolveIters_last;
    }
    if (status.extraParameters->isParameter("AztecOO/Achieved Tolerance")) 
      linearSolveAchievedTol = status.extraParameters->get<double>("AztecOO/Achieved Tolerance");
  }

  // Dump solution of linear system
#ifdef HAVE_NOX_EPETRAEXT
  if (p.get("Write Linear System", false)) {
    std::ostringstream iterationNumber;
    iterationNumber << linearSolveCount;
    
    std::string prefixName = p.get("Write Linear System File Prefix", 
				   "NOX_LinSys");
    std::string postfixName = iterationNumber.str();
    postfixName += ".mm";

    std::string lhsFileName = prefixName + "_LHS_" + postfixName;
    std::string rhsFileName = prefixName + "_RHS_" + postfixName;
    std::string jacFileName = prefixName + "_Jacobian_" + postfixName;
    EpetraExt::MultiVectorToMatrixMarketFile(lhsFileName.c_str(), 
					     result.getEpetraVector());
    EpetraExt::MultiVectorToMatrixMarketFile(rhsFileName.c_str(), 
					     input.getEpetraVector());

    Epetra_RowMatrix* printMatrix = NULL;
    printMatrix = dynamic_cast<Epetra_RowMatrix*>(jacPtr.get()); 
    if (printMatrix == NULL) {
      std::cout << "Error: NOX::Epetra::LinearSystemAztecOO::applyJacobianInverse() - "
	   << "Could not cast the Jacobian operator to an Epetra_RowMatrix!"
	   << "Please set the \"Write Linear System\" parameter to false."
	   << std::endl;
//.........这里部分代码省略.........

void
LOCA::Epetra::Group::printSolution(const NOX::Epetra::Vector& x_,
				      const double conParam) const
{
  userInterface->printSolution(x_.getEpetraVector(), conParam);
}

double NOX::Epetra::Vector::innerProduct(const NOX::Epetra::Vector& y) const
{
  return vectorSpace->innerProduct(*epetraVec, y.getEpetraVector());
}