C++ MultiVector::GetVectorSpace方法代码示例

/*----------------------------------------------------------------------*
 |  apply multigrid linear preconditioner (private)          m.gee 03/06|
 *----------------------------------------------------------------------*/
int MOERTEL::Mortar_ML_Preconditioner::MultiLevelSA(
    const MultiVector& b1_f,
    const MultiVector& b2_f,
    MultiVector& x1_f,
    MultiVector& x2_f,
    int level) const
{
    MultiVector r1_f(b1_f.GetVectorSpace(),1,false);
    MultiVector z1_f(b1_f.GetVectorSpace(),1,false);

    // presmoothing
    x1_f = 0;
    G_.Apply(b1_f,x1_f);
    x2_f = mlapiMT_ * x1_f;
    x2_f.Scale(-1.0);

    // compute residual
    r1_f = b1_f - mlapiAhat11_ * x1_f;

    // postsmoothing
    z1_f = 0;
    G_.Apply(r1_f,z1_f);
    x1_f = x1_f + z1_f;
    x2_f = mlapiMT_ * x1_f;
    x2_f.Scale(-1.0);



    return 0;
}

// ====================================================================== 
Operator GetPtent1D(const MultiVector& D, const int offset = 0)
{
  if (D.GetNumVectors() != 1)
    ML_THROW("D.GetNumVectors() != 1", -1);

  int size = D.GetMyLength();
  if (size == 0)
    ML_THROW("empty diagonal vector in input", -1);

  double* diag = new double[size];
  for (int i = 0 ; i < size ; ++i)
    diag[i] = D(i);

  // creates the ML operator and store the diag pointer,
  // as well as the function pointers
  ML_Operator* MLDiag = ML_Operator_Create(GetML_Comm());

  int invec_leng = size / 3 + size % 3;
  int outvec_leng = size;

  MLDiag->invec_leng = invec_leng;
  MLDiag->outvec_leng = outvec_leng;
  MLDiag->data = (void*)diag;
  MLDiag->data_destroy = Ptent1D_destroy;
  MLDiag->matvec->func_ptr = Ptent1D_matvec;

  MLDiag->matvec->ML_id = ML_NONEMPTY;
  MLDiag->matvec->Nrows = outvec_leng;
  MLDiag->from_an_ml_operator = 0;

  MLDiag->getrow->func_ptr = Ptent1D_getrows;

  MLDiag->getrow->ML_id = ML_NONEMPTY;
  MLDiag->getrow->Nrows = outvec_leng;

  // creates the domain space
  vector<int> MyGlobalElements(invec_leng);
  for (int i = 0 ; i < invec_leng ; ++i) 
    MyGlobalElements[i] = D.GetVectorSpace()(i * 3) / 3;
  Space DomainSpace(invec_leng, -1, &MyGlobalElements[0]);
  Space RangeSpace = D.GetVectorSpace();

  // creates the MLAPI wrapper
  Operator Diag(DomainSpace,RangeSpace,MLDiag,true);
  return(Diag);
}

int 
InverseOperator::Apply(const MultiVector& x, MultiVector& y) const
{
  ResetTimer();
  StackPush();

  if (GetDomainSpace() != x.GetVectorSpace())
    ML_THROW("DomainSpace and x.GetVectorSpace() differ", -1);

  if (GetRangeSpace() != y.GetVectorSpace())
    ML_THROW("RangeSpace and y.GetVectorSpace() differ", -1);

  int x_nv = x.GetNumVectors();
  int y_nv = y.GetNumVectors();
  double FL = 0.0;
  if (RCPData_ != Teuchos::null)
    FL = RCPData_->ComputeFlops();

  if (x_nv != y_nv)
    ML_THROW("Number of vectors of x and y differ (" +
             GetString(x_nv) + " vs. " + GetString(x_nv), -1);

  for (int v = 0 ; v < x_nv ; ++v) {

    Epetra_Vector x_Epetra(View,RowMatrix()->OperatorDomainMap(),
                           (double*)&(x(0,v)));
    Epetra_Vector y_Epetra(View,RowMatrix()->OperatorRangeMap(),
                           (double*)&(y(0,v)));

    if (RCPData_ != Teuchos::null)
      RCPData_->ApplyInverse(x_Epetra,y_Epetra);
    else if (RCPMLPrec_ != Teuchos::null)
      RCPMLPrec_->ApplyInverse(x_Epetra,y_Epetra);
    else
      ML_THROW("Neither Ifpack nor ML smoother is properly set up", -1);
  }

  StackPop();
  if (RCPData_ != Teuchos::null)
    UpdateFlops(RCPData_->ComputeFlops() - FL);
  UpdateTime();

  return(0);
}

// ======================================================================
MultiVector Extract(const MultiVector& y, const int v)
{
  if ((v < 0) || v >= y.GetNumVectors())
    ML_THROW("Wrong input parameter v (" +
             GetString(v) + ")", -1);

  MultiVector x(y.GetVectorSpace(), y.GetRCPValues(v));

  return(x);
}

MultiVector LoadBalanceInverseOperator::operator()(const MultiVector& LHS)
{
  StackPush();

  MultiVector RHS(LHS.GetVectorSpace());
  RHS = 0.0;
  Apply(LHS,RHS);

  StackPop();

  return(RHS);
}

// ======================================================================
MultiVector Duplicate(const MultiVector& y, const int v)
{
  if ((v < 0) || v >= y.GetNumVectors())
    ML_THROW("Wrong input parameter v (" +
             GetString(v) + ")", -1);

  // FIXME: use Extract
  MultiVector x(y.GetVectorSpace(), 1);
  for (int i = 0 ; i < x.GetMyLength() ; ++i)
    x(i) = y(i,v);

  return(x);
}

// ====================================================================== 
Operator GetDiagonal(const MultiVector& D)
{
  if (D.GetNumVectors() != 1)
    ML_THROW("D.GetNumVectors() != 1", -1);

  int size = D.GetMyLength();
  if (size == 0)
    ML_THROW("empty diagonal vector in input", -1);

  double* diag = new double[size];
  for (int i = 0 ; i < size ; ++i) 
    diag[i] = D(i);

  // creates the ML operator and store the diag pointer,
  // as well as the function pointers
  ML_Operator* MLDiag = ML_Operator_Create(GetML_Comm());

  MLDiag->invec_leng = size;
  MLDiag->outvec_leng = size;
  MLDiag->data = (void*)diag;
  MLDiag->matvec->func_ptr = diag_matvec;

  MLDiag->matvec->ML_id = ML_NONEMPTY;
  MLDiag->matvec->Nrows = size;
  MLDiag->from_an_ml_operator = 0;
  MLDiag->data_destroy = diag_destroy;

  MLDiag->getrow->func_ptr = diag_getrows;

  MLDiag->getrow->ML_id = ML_NONEMPTY;
  MLDiag->getrow->Nrows = size;

  // creates the MLAPI wrapper
  Operator Diag(D.GetVectorSpace(),D.GetVectorSpace(),MLDiag,true);
  return(Diag);
}

/*----------------------------------------------------------------------*
 |  apply multigrid linear preconditioner (private)          m.gee 03/06|
 *----------------------------------------------------------------------*/
int MOERTEL::Mortar_ML_Preconditioner::MultiLevelSA(const MultiVector& b_f,
        MultiVector& x_f,
        int level) const
{
    if (level == maxlevels_-1)
    {
        x_f = 0;
        S(level).Apply(b_f,x_f);
        x_f = ImWBT(level) * x_f;
        return 0;
    }

    MultiVector r_f(b_f.GetVectorSpace(),1,false);
    MultiVector z_f(b_f.GetVectorSpace(),1,false);
    MultiVector r_c(P(level).GetDomainSpace(),1,false);
    MultiVector z_c(P(level).GetDomainSpace(),1,false);

    // presmoothing
    x_f = 0;
    S(level).Apply(b_f,x_f);
    x_f = ImWBT(level) * x_f;


    // compute residual (different operator)
    r_f = b_f - Ahat(level) * x_f;

    // restrict
    r_c = R(level) * r_f;

    // solve coarser problem
    MultiLevelSA(r_c,z_c,level+1);

    // prolongate
    x_f = x_f + P(level) * z_c;
    x_f = ImWBT(level) * x_f;

    // recompute residual using a different operator
    r_f = b_f - Ahat(level) * x_f;

    // postsmoothing
    z_f = 0;
    S(level).Apply(r_f,z_f);
    z_f = ImWBT(level) * z_f;
    x_f = x_f + z_f;

    return 0;
}