C++ Epetra_MultiVector::NumVectors方法代码示例

//=============================================================================
int Epetra_MsrMatrix::Multiply(bool TransA,
                               const Epetra_MultiVector& X,
                               Epetra_MultiVector& Y) const
{
  (void)TransA;
  int NumVectors = X.NumVectors();
  if (NumVectors!=Y.NumVectors()) EPETRA_CHK_ERR(-1);  // X and Y must have same number of vectors

  double ** xptrs;
  double ** yptrs;
  X.ExtractView(&xptrs);
  Y.ExtractView(&yptrs);
  if (RowMatrixImporter()!=0) {
    if (ImportVector_!=0) {
      if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
    }
    if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(RowMatrixColMap(),NumVectors); // Create import vector if needed
    ImportVector_->Import(X, *RowMatrixImporter(), Insert);
    ImportVector_->ExtractView(&xptrs);
  }
  for (int i=0; i<NumVectors; i++)
    Amat_->matvec(xptrs[i], yptrs[i], Amat_, proc_config_);
  
  double flops = NumGlobalNonzeros();
  flops *= 2.0;
  flops *= (double) NumVectors;
  UpdateFlops(flops);
  return(0);
}

//=============================================================================
// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS
int Ifpack_SPARSKIT::ApplyInverse(const Epetra_MultiVector& X, 
                                  Epetra_MultiVector& Y) const
{
  if (!IsComputed())
    IFPACK_CHK_ERR(-3); // compute preconditioner first

  if (X.NumVectors() != Y.NumVectors()) 
    IFPACK_CHK_ERR(-2); // Return error: X and Y not the same size

  int n = Matrix().NumMyRows();

  for (int i = 0 ; i < X.NumVectors() ; ++i)
    F77_LUSOL(&n, (double*)X(i)->Values(), Y(i)->Values(), (double*)&alu_[0], 
                (int*)&jlu_[0], (int*)&ju_[0]);

  // still need to fix support for permutation
  if (Type_ == "ILUTP" || Type_ == "ILUDP")
  {
    vector<double> tmp(n);
    for (int j = 0 ; j < n ; ++j)
      tmp[iperm_[j]] = Y[0][j];
    for (int j = 0 ; j < n ; ++j)
      Y[0][j] = tmp[j];
  }

  ++NumApplyInverse_;
  return(0);

}

// ============================================================================ 
int ML_Epetra::MatrixFreePreconditioner::
ApplyInvBlockDiag(const double alpha, Epetra_MultiVector& X,
                  const double beta, const Epetra_MultiVector& B) const
{
  assert (X.NumVectors() == 1);
  int NumPDEEqns2 = NumPDEEqns_ * NumPDEEqns_;

  char trans = 'N';
  int NumVectorsX = X.NumVectors();
  std::vector<double> tmp(NumPDEEqns_);

  size_t len = sizeof(double) * NumPDEEqns_;
  for (int i = 0; i < NumMyBlockRows_; ++i)
  {
    memcpy(&tmp[0], &(B[0][i * NumPDEEqns_]), len);

    int offset = i * NumPDEEqns2;

    DGEMM_F77(&trans, &trans, (int*)&NumPDEEqns_, &NumVectorsX, (int*)&NumPDEEqns_,
              (double*)&alpha, (double*)&InvBlockDiag_[offset], (int*)&NumPDEEqns_, 
              &tmp[0], (int*)&NumPDEEqns_, (double*)&beta, 
              (double*)&X[0][i * NumPDEEqns_], (int*)&NumPDEEqns_);
  }

  return(0);
}

Epetra_OskiMultiVector::Epetra_OskiMultiVector(const Epetra_MultiVector& Source) 
  : Epetra_MultiVector(Source),
  Epetra_View_(&Source), 
  Copy_Created_(false) {
    double* A;
    double** Aptr;
    int LDA;
    int* LDAptr;
    LDAptr = new int[1];
    Aptr = new double*[1];
    if(Source.ConstantStride() || (Source.NumVectors() == 1)) {
      if(Source.ExtractView(Aptr, LDAptr))
        std::cerr << "Extract view failed\n";
      else
        Oski_View_ = oski_CreateMultiVecView(*Aptr, Source.MyLength(), Source.NumVectors(), LAYOUT_COLMAJ, *LDAptr);
    }
    else {
      Copy_Created_ = true;
      LDA = Source.MyLength();
      A = new double[LDA*Source.NumVectors()];
      if(Source.ExtractCopy(A, LDA))
        std::cerr << "Extract copy failed\n";
      else
        Oski_View_ = oski_CreateMultiVecView(A, Source.MyLength(), Source.NumVectors(), LAYOUT_COLMAJ, LDA);
    }
    delete [] LDAptr;
    delete [] Aptr;
}

void
LOCA::Epetra::CompactWYOp::applyCompactWY(const Epetra_MultiVector& x, 
					  Epetra_MultiVector& result_x, 
					  Epetra_MultiVector& result_p) const
{
  // Compute Y_x^T*x
  result_p.Multiply('T', 'N', 1.0, *Y_x, x, 0.0);

  // Compute T*(Y_x^T*x)
  dblas.TRMM(Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS, 
	     Teuchos::NON_UNIT_DIAG, result_p.MyLength(), 
	     result_p.NumVectors(), 1.0, T.Values(), T.MyLength(), 
	     result_p.Values(), result_p.MyLength());

  // Compute x = x + Y_x*T*(Y_x^T*x)
  result_x = x;
  result_x.Multiply('N', 'N', 1.0, *Y_x, result_p, 1.0);

  // Compute result_p = Y_p*T*(Y_x^T*x)
  dblas.TRMM(Teuchos::LEFT_SIDE, Teuchos::LOWER_TRI, Teuchos::NO_TRANS, 
	     Teuchos::UNIT_DIAG, result_p.MyLength(), 
	     result_p.NumVectors(), 1.0, Y_p.Values(), Y_p.MyLength(), 
	     result_p.Values(), result_p.MyLength());

}

//==============================================================================
int Ifpack_ML::
ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{

  if (IsComputed() == false)
    IFPACK_CHK_ERR(-1);

  int numVectors = X.NumVectors();
  if (numVectors != Y.NumVectors())
    IFPACK_CHK_ERR(-1); // wrong input

  
  Time_->ResetStartTime();

  // AztecOO gives X and Y pointing to the same memory location,
  // need to create an auxiliary vector, Xcopy
  Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
  if (X.Pointers()[0] == Y.Pointers()[0])
    Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
  else
    Xcopy = Teuchos::rcp( &X, false );

  for (int i=0; i<numVectors; i++) {
    
  IFPACK_CHK_ERR(MLPrec_->ApplyInverse(*Xcopy,Y));
  }

  ++NumApplyInverse_;
  ApplyInverseTime_ += Time_->ElapsedTime();

  return(0);
}

void PeridigmNS::State::copyLocallyOwnedMultiVectorData(Epetra_MultiVector& source, Epetra_MultiVector& target)
{
  TEUCHOS_TEST_FOR_EXCEPTION(source.NumVectors() != target.NumVectors(), std::runtime_error,
                     "PeridigmNS::State::copyLocallyOwnedMultiVectorData() called with incompatible MultiVectors.\n");
  int numVectors = target.NumVectors();
  const Epetra_BlockMap& sourceMap = source.Map();
  const Epetra_BlockMap& targetMap = target.Map();
  for(int iVec=0 ; iVec<numVectors ; ++iVec){
    Epetra_Vector& sourceVector = *source(iVec);
    Epetra_Vector& targetVector = *target(iVec);
    for(int targetLID=0 ; targetLID<targetMap.NumMyElements() ; ++targetLID){
      int GID = targetMap.GID(targetLID);
      int sourceLID = sourceMap.LID(GID);
      TEUCHOS_TEST_FOR_EXCEPTION(sourceLID == -1, std::range_error,
                         "PeridigmNS::State::copyLocallyOwnedMultiVectorData() called with incompatible MultiVectors.\n");
      TEUCHOS_TEST_FOR_EXCEPTION(sourceMap.ElementSize(sourceLID) != targetMap.ElementSize(targetLID), std::range_error,
                         "PeridigmNS::State::copyLocallyOwnedMultiVectorData() called with incompatible MultiVectors.\n");
      int elementSize = targetMap.ElementSize(targetLID);
      int sourceFirstPointInElement = sourceMap.FirstPointInElement(sourceLID);
      int targetFirstPointInElement = targetMap.FirstPointInElement(targetLID);
      for(int i=0 ; i<elementSize ; ++i){
        targetVector[targetFirstPointInElement+i] = sourceVector[sourceFirstPointInElement+i];
      }
    }
  }
}

//=========================================================================
int Ifpack_CrsRiluk::GenerateXY(bool Trans, 
				const Epetra_MultiVector& Xin, const Epetra_MultiVector& Yin,
				Teuchos::RefCountPtr<Epetra_MultiVector>* Xout, 
				Teuchos::RefCountPtr<Epetra_MultiVector>* Yout) const {

  // Generate an X and Y suitable for performing Solve() and Multiply() methods

  if (Xin.NumVectors()!=Yin.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size

  //cout << "Xin = " << Xin << endl;
  (*Xout) = Teuchos::rcp( (Epetra_MultiVector *) &Xin, false );
  (*Yout) = Teuchos::rcp( (Epetra_MultiVector *) &Yin, false );
  if (!IsOverlapped_ && UserMatrixIsCrs_) return(0); // Nothing more to do

  if (UserMatrixIsVbr_) {
    if (VbrX_!=Teuchos::null) {
      if (VbrX_->NumVectors()!=Xin.NumVectors()) {
	VbrX_ = Teuchos::null;
	VbrY_ = Teuchos::null;
      }
    }
    if (VbrX_==Teuchos::null) { // Need to allocate space for overlap X and Y
      VbrX_ = Teuchos::rcp( new Epetra_MultiVector(View, *U_DomainMap_, (*Xout)->Pointers(), (*Xout)->NumVectors()) );
      VbrY_ = Teuchos::rcp( new Epetra_MultiVector(View, *L_RangeMap_, (*Yout)->Pointers(), (*Yout)->NumVectors()) );
    }
    else {
      EPETRA_CHK_ERR(VbrX_->ResetView((*Xout)->Pointers()));
      EPETRA_CHK_ERR(VbrY_->ResetView((*Yout)->Pointers()));
    }
    (*Xout) = VbrX_;
    (*Yout) = VbrY_;
  }
    
  if (IsOverlapped_) {
    // Make sure the number of vectors in the multivector is the same as before.
    if (OverlapX_!=Teuchos::null) {
      if (OverlapX_->NumVectors()!=Xin.NumVectors()) {
	OverlapX_ = Teuchos::null;
	OverlapY_ = Teuchos::null;
      }
    }
    if (OverlapX_==Teuchos::null) { // Need to allocate space for overlap X and Y
      OverlapX_ = Teuchos::rcp( new Epetra_MultiVector(U_->RowMatrixColMap(), (*Xout)->NumVectors()) );
      OverlapY_ = Teuchos::rcp( new Epetra_MultiVector(L_->RowMatrixRowMap(), (*Yout)->NumVectors()) );
    }
    if (!Trans) {
      EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*U_->Importer(), Insert)); // Import X values for solve
    }
    else {
      EPETRA_CHK_ERR(OverlapX_->Import(*(*Xout),*L_->Exporter(), Insert)); // Import X values for solve
    }
    (*Xout) = OverlapX_;
    (*Yout) = OverlapY_; // Set pointers for Xout and Yout to point to overlap space
    //cout << "OverlapX_ = " << *OverlapX_ << endl;
  }
  
  return(0);
}

//=======================================================
int EpetraExt_HypreIJMatrix::Multiply(bool TransA,
                               const Epetra_MultiVector& X,
                               Epetra_MultiVector& Y) const
{
  
  //printf("Proc[%d], Row start: %d, Row End: %d\n", Comm().MyPID(), MyRowStart_, MyRowEnd_);
  bool SameVectors = false; 
  int NumVectors = X.NumVectors();
  if (NumVectors != Y.NumVectors()) return -1;  // X and Y must have same number of vectors
  if(X.Pointers() == Y.Pointers()){
    SameVectors = true;
  }
  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
    //Get values for current vector in multivector.
    double * x_values;
    double * y_values;
    EPETRA_CHK_ERR((*X(VecNum)).ExtractView(&x_values));
    double *x_temp = x_local->data; 
    double *y_temp = y_local->data;
    if(!SameVectors){
      EPETRA_CHK_ERR((*Y(VecNum)).ExtractView(&y_values));
    } else {
      y_values = new double[X.MyLength()];
    }
    y_local->data = y_values;
    EPETRA_CHK_ERR(HYPRE_ParVectorSetConstantValues(par_y,0.0));
    // Temporarily make a pointer to data in Hypre for end
    // Replace data in Hypre vectors with epetra values
    x_local->data = x_values;
    // Do actual computation.
    if(TransA) {
      // Use transpose of A in multiply
      EPETRA_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, par_x, 1.0, par_y));
    } else {
      EPETRA_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, par_x, 1.0, par_y));
    }
    if(SameVectors){
      int NumEntries = Y.MyLength();
      std::vector<double> new_values; new_values.resize(NumEntries);
      std::vector<int> new_indices; new_indices.resize(NumEntries);
      for(int i = 0; i < NumEntries; i++){
        new_values[i] = y_values[i];
        new_indices[i] = i;
      }
      EPETRA_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
      delete[] y_values;
    }
    x_local->data = x_temp;
    y_local->data = y_temp;
  }
  double flops = (double) NumVectors * (double) NumGlobalNonzeros();
  UpdateFlops(flops);
  return 0;
} //Multiply() 

// ============================================================================ 
 int Ifpack_DiagPreconditioner::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{
  if (X.NumVectors() != Y.NumVectors())
    IFPACK_CHK_ERR(-1); 

  for (int v = 0; v < X.NumVectors(); ++v)
    for (int i = 0; i < X.MyLength(); ++i)
      Y[v][i] = diag_[i] * X[v][i];
  ///Y.ReciprocalMultiply(1.0, diag_, X, 0.0);

  return(0);
}

int BlockPCGSolver::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const {

  int xcol = X.NumVectors();
  int info = 0;

  if (Y.NumVectors() < xcol)
    return -1;

  // Use block PCG for multiple right-hand sides
  info = (xcol == 1) ? Solve(X, Y) : Solve(X, Y, xcol);

  return info;
}

//==============================================================================
int Ifpack_Hypre::Multiply(bool TransA, const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
  if(IsInitialized() == false){
    IFPACK_CHK_ERR(-1);
  }
  bool SameVectors = false;
  int NumVectors = X.NumVectors();
  if (NumVectors != Y.NumVectors()) IFPACK_CHK_ERR(-1);  // X and Y must have same number of vectors
  if(X.Pointers() == Y.Pointers()){
    SameVectors = true;
  }
  for(int VecNum = 0; VecNum < NumVectors; VecNum++) {
    //Get values for current vector in multivector.
    double * XValues;
    double * YValues;
    IFPACK_CHK_ERR((*X(VecNum)).ExtractView(&XValues));
    double *XTemp = XLocal_->data;
    double *YTemp = YLocal_->data;
    if(!SameVectors){
      IFPACK_CHK_ERR((*Y(VecNum)).ExtractView(&YValues));
    } else {
      YValues = new double[X.MyLength()];
    }
    YLocal_->data = YValues;
    IFPACK_CHK_ERR(HYPRE_ParVectorSetConstantValues(ParY_,0.0));
    // Temporarily make a pointer to data in Hypre for end
    // Replace data in Hypre vectors with epetra values
    XLocal_->data = XValues;
    // Do actual computation.
    if(TransA) {
      // Use transpose of A in multiply
      IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvecT(1.0, ParMatrix_, ParX_, 1.0, ParY_));
    } else {
      IFPACK_CHK_ERR(HYPRE_ParCSRMatrixMatvec(1.0, ParMatrix_, ParX_, 1.0, ParY_));
    }
    if(SameVectors){
      int NumEntries = Y.MyLength();
      std::vector<double> new_values; new_values.resize(NumEntries);
      std::vector<int> new_indices; new_indices.resize(NumEntries);
      for(int i = 0; i < NumEntries; i++){
        new_values[i] = YValues[i];
        new_indices[i] = i;
      }
      IFPACK_CHK_ERR((*Y(VecNum)).ReplaceMyValues(NumEntries, &new_values[0], &new_indices[0]));
      delete[] YValues;
    }
    XLocal_->data = XTemp;
    YLocal_->data = YTemp;
  }
  return 0;
} //Multiply()

int Epetra_PETScAIJMatrix::Multiply(bool TransA,
                               const Epetra_MultiVector& X,
                               Epetra_MultiVector& Y) const
{
  (void)TransA;
  int NumVectors = X.NumVectors();
  if (NumVectors!=Y.NumVectors()) EPETRA_CHK_ERR(-1);  // X and Y must have same number of vectors

  double ** xptrs;
  double ** yptrs;
  X.ExtractView(&xptrs);
  Y.ExtractView(&yptrs);
  if (RowMatrixImporter()!=0) {
    if (ImportVector_!=0) {
      if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
    }
    if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(RowMatrixColMap(),NumVectors);
    ImportVector_->Import(X, *RowMatrixImporter(), Insert);
    ImportVector_->ExtractView(&xptrs);
  }

  double *vals=0;
  int length;
  Vec petscX, petscY;
  int ierr;
  for (int i=0; i<NumVectors; i++) {
#   ifdef HAVE_MPI
    ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
    ierr=VecCreateMPIWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
#   else //FIXME  untested
    ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
    ierr=VecCreateSeqWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
#   endif

    ierr = MatMult(Amat_,petscX,petscY);CHKERRQ(ierr);

    ierr = VecGetArray(petscY,&vals);CHKERRQ(ierr);
    ierr = VecGetLocalSize(petscY,&length);CHKERRQ(ierr);
    for (int j=0; j<length; j++) yptrs[i][j] = vals[j];
    ierr = VecRestoreArray(petscY,&vals);CHKERRQ(ierr);
  }

  VecDestroy(petscX); VecDestroy(petscY);
  
  double flops = NumGlobalNonzeros();
  flops *= 2.0;
  flops *= (double) NumVectors;
  UpdateFlops(flops);
  return(0);
} //Multiply()

//=========================================================================  
// Returns the result of a Epetra_Operator inverse applied to an Epetra_MultiVector X in Y.
int EpetraExt_PointToBlockDiagPermute::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
  // Stuff borrowed from Epetra_CrsMatrix
  int NumVectors = X.NumVectors();
  if (NumVectors!=Y.NumVectors()) {
    EPETRA_CHK_ERR(-2); // Need same number of vectors in each MV
  }

  const Epetra_MultiVector *Xp=&X;
  Epetra_MultiVector *Yp=&Y;

  // Allocate temp workspace if X==Y and there are no imports or exports
  Epetra_MultiVector * Xcopy = 0;
  if (&X==&Y && Importer_==0 && Exporter_==0) {
    Xcopy = new Epetra_MultiVector(X);
    Xp=Xcopy;
  }
  
  UpdateImportVector(NumVectors); // Make sure Import and Export Vectors are compatible
  UpdateExportVector(NumVectors);

  // If we have a non-trivial importer, we must import elements that are permuted or are on other processors
  if (Importer_){
    EPETRA_CHK_ERR(ImportVector_->Import(X, *Importer_, Insert));
    Xp=ImportVector_;
  }
  
  // If we have a non-trivial exporter, we must export elements that are permuted or belong to other processors
  if (Exporter_) {
    Yp=ExportVector_;
  }
  
  // Do the matvec 
  BDMat_->ApplyInverse(*Xp,*Yp);

  // Export if needed
  if (Exporter_) {
    Y.PutScalar(0.0);  // Make sure target is zero
    Y.Export(*ExportVector_, *Exporter_, Add); // Fill Y with Values from export vector
  }
  
  // Cleanup
  if(Xcopy) {
    delete Xcopy;
    EPETRA_CHK_ERR(1); // Return positive code to alert the user about needing extra copy of X
    return 1;
  }

  return 0;
}

void
LOCA::Epetra::AugmentedOp::init(const Epetra_MultiVector& x)
{
  if (importedInput == NULL || importedInput->NumVectors() != x.NumVectors()) {
    if (importedInput != NULL) {
      delete importedInput;
      delete result_y;
      delete tmp;
    }
    importedInput = new Epetra_MultiVector(*extendedImportMapPtr, 
					   x.NumVectors());
    result_y = new Epetra_MultiVector(localMap, x.NumVectors());
    tmp = new Epetra_MultiVector(underlyingMap, x.NumVectors());
  }
}