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

int Ifpack_IHSS::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{  
  if(!IsComputed_) return -1;
  Time_.ResetStartTime();
  bool initial_guess_is_zero=false;  

  // AztecOO gives X and Y pointing to the same memory location,
  // need to create an auxiliary vector, Xcopy
  Teuchos::RefCountPtr<const Epetra_MultiVector> Xcopy;
  Epetra_MultiVector Temp(X);
  if (X.Pointers()[0] == Y.Pointers()[0]){
    Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
    // Since the user didn't give us anything better, our initial guess is zero.
    Y.Scale(0.0);
    initial_guess_is_zero=true;
  }
  else
    Xcopy = Teuchos::rcp( &X, false );

  Epetra_MultiVector T1(Y),T2(*Xcopy);

  // Note: Since Aherm and Askew are actually (aI+H) and (aI+S) accordingly (to save memory), 
  // the application thereof needs to be a little different.
  // The published algorithm is:
  // temp = (aI+H)^{-1} [ (aI-S) y + x ]
  // y = (aI+S)^{-1} [ (aI-H) temp + x ]
  //
  // But we're doing:
  // temp = (aI+H)^{-1} [ 2a y - Shat y + x ]
  // y = (aI+S)^{-1} [ 2 a temp - Hhat temp + x ]

  for(int i=0;i<NumSweeps_;i++){
    // temp = (aI+H)^{-1} [ 2a y - Shat y + x ]
    if(!initial_guess_is_zero || i >0 ){
      Askew_->Apply(Y,T1);
      T2.Update(2*Alpha_,Y,-1,T1,1);
    }
    Pherm_->ApplyInverse(T2,Y);
    
    // y = (aI+S)^{-1} [ 2 a temp - Hhat temp + x ]
    Aherm_->Apply(Y,T1);
    T2.Scale(1.0,*Xcopy);
    T2.Update(2*Alpha_,Y,-1,T1,1.0);
    Pskew_->ApplyInverse(T2,Y);  
  }

  // Counter update
  NumApplyInverse_++;
  ApplyInverseTime_ += Time_.ElapsedTime();
  return 0;
}

void
Stokhos::EpetraMultiVectorOrthogPoly::
computeMean(Epetra_MultiVector& v) const
{
  v.Scale(1.0, *(coeff_[0]));
}

int Ifpack_SORa::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
  if(!IsComputed_) return -1;
  Time_.ResetStartTime();
  bool initial_guess_is_zero=false;
  const int lclNumRows = W_->NumMyRows();
  const int NumVectors = X.NumVectors();
  Epetra_MultiVector Temp(A_->RowMatrixRowMap(),NumVectors);

  double omega=GetOmega();

  // need to create an auxiliary vector, Xcopy
  Teuchos::RCP<const Epetra_MultiVector> Xcopy;
  if (X.Pointers()[0] == Y.Pointers()[0]){
    Xcopy = Teuchos::rcp( new Epetra_MultiVector(X) );
    // Since the user didn't give us anything better, our initial guess is zero.
    Y.Scale(0.0);
    initial_guess_is_zero=true;
  }
  else
    Xcopy = Teuchos::rcp( &X, false );

  Teuchos::RCP< Epetra_MultiVector > T2;
  // Note: Assuming that the matrix has an importer.  Ifpack_PointRelaxation doesn't do this, but given that
  // I have a CrsMatrix, I'm probably OK.
  // Note: This is the lazy man's version sacrificing a few extra flops for avoiding if statements to determine
  // if things are on or off processor.
  // Note: T2 must be zero'd out
  if (IsParallel_ && W_->Importer())  T2 = Teuchos::rcp( new Epetra_MultiVector(W_->Importer()->TargetMap(),NumVectors,true));
  else T2 = Teuchos::rcp( new Epetra_MultiVector(A_->RowMatrixRowMap(),NumVectors,true));

  // Pointer grabs
  int* rowptr,*colind;
  double *values;
  double **t_ptr,** y_ptr, ** t2_ptr, **x_ptr,*d_ptr;
  T2->ExtractView(&t2_ptr);
  Y.ExtractView(&y_ptr);
  Temp.ExtractView(&t_ptr);
  Xcopy->ExtractView(&x_ptr);
  Wdiag_->ExtractView(&d_ptr);
  IFPACK_CHK_ERR(W_->ExtractCrsDataPointers(rowptr,colind,values));


  for(int i=0; i<NumSweeps_; i++){
    // Calculate b-Ax
    if(!initial_guess_is_zero  || i > 0) {
      A_->Apply(Y,Temp);
      Temp.Update(1.0,*Xcopy,-1.0);
    }
    else
      Temp.Update(1.0,*Xcopy,0.0);

    // Note: The off-processor entries of T2 never get touched (they're always zero) and the other entries are updated
    // in this sweep before they are used, so we don't need to reset T2 to zero here.

    // Do backsolve & update
    // x = x  + W^{-1} (b - A x)
    for(int j=0; j<lclNumRows; j++){
      double diag=d_ptr[j];
      for (int m=0 ; m<NumVectors; m++) {
        double dtmp=0.0;
        // Note: Since the diagonal is in the matrix, we need to zero that entry of T2 here to make sure it doesn't contribute.
        t2_ptr[m][j]=0.0;
        for(int k=rowptr[j];k<rowptr[j+1];k++){
          dtmp+= values[k]*t2_ptr[m][colind[k]];
        }
        // Yes, we need to update both of these.
        t2_ptr[m][j] = (t_ptr[m][j]- dtmp)/diag;
        y_ptr[m][j] += omega*t2_ptr[m][j];
      }
    }
  }

  // Counter update
  NumApplyInverse_++;
  ApplyInverseTime_ += Time_.ElapsedTime();
  return 0;
}

int Amesos_Scalapack::Solve() { 
  
  if( debug_ == 1 ) std::cout << "Entering `Solve()'" << std::endl;
  
  NumSolve_++;
  
  Epetra_MultiVector   *vecX = Problem_->GetLHS() ; 
  Epetra_MultiVector   *vecB = Problem_->GetRHS() ; 
  
  //
  //  Compute the number of right hands sides 
  //  (and check that X and B have the same shape) 
  //
  int nrhs; 
  if ( vecX == 0 ) { 
    nrhs = 0 ;
    EPETRA_CHK_ERR( vecB != 0 ) ; 
  } else { 
    nrhs = vecX->NumVectors() ; 
    EPETRA_CHK_ERR( vecB->NumVectors() != nrhs ) ; 
  }
  
  Epetra_MultiVector *ScalapackB =0;
  Epetra_MultiVector *ScalapackX =0;
  //
  //  Extract Scalapack versions of X and B 
  //
  double *ScalapackXvalues ;
  
  Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
  Time_->ResetStartTime(); // track time to broadcast vectors
  //
  //  Copy B to the scalapack version of B
  //
  const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap();
  Epetra_MultiVector *ScalapackXextract = new Epetra_MultiVector( *VectorMap_, nrhs ) ; 
  Epetra_MultiVector *ScalapackBextract = new Epetra_MultiVector( *VectorMap_, nrhs ) ; 
  
  Epetra_Import ImportToScalapack( *VectorMap_, OriginalMap );
  ScalapackBextract->Import( *vecB, ImportToScalapack, Insert ) ;
  ScalapackB = ScalapackBextract ; 
  ScalapackX = ScalapackXextract ; 
  
  VecTime_ += Time_->ElapsedTime();
  
  //
  //  Call SCALAPACKs PDGETRS to perform the solve
  //
  
  int DescX[10];  
  
  ScalapackX->Scale(1.0, *ScalapackB) ;  
  
  int ScalapackXlda ; 
  
  Time_->ResetStartTime(); // tract time to solve
  
  //
  //  Setup DescX 
  //
  
  if( nrhs > nb_ ) {
    EPETRA_CHK_ERR( -2 );  
  }
  
  int Ierr[1] ; 
  Ierr[0] = 0 ; 
  const int zero = 0 ; 
  const int one = 1 ; 
  if ( iam_ < nprow_ * npcol_ ) {
    assert( ScalapackX->ExtractView( &ScalapackXvalues, &ScalapackXlda ) == 0 ) ; 
    
    if ( false ) std::cout << "Amesos_Scalapack.cpp: " << __LINE__ << " ScalapackXlda = "  <<  ScalapackXlda 
		      << " lda_ = "  << lda_ 
		      << " nprow_ = "  << nprow_ 
		      << " npcol_ = "  << npcol_ 
		      << " myprow_ = "  << myprow_ 
		      << " mypcol_ = "  << mypcol_ 
		      << " iam_ = "  << iam_ << std::endl ;
    if (  TwoD_distribution_ )    assert( mypcol_ >0 || EPETRA_MAX(ScalapackXlda,1) == lda_ ) ; 
    
    DESCINIT_F77(DescX, 
		 &NumGlobalElements_, 
		 &nrhs, 
		 &nb_,
		 &nb_,
		 &zero,
		 &zero,
		 &ictxt_,
		 &lda_,
		 Ierr ) ;
    assert( Ierr[0] == 0 ) ; 
    
    //
    //  For the 1D data distribution, we factor the transposed 
    //  matrix, hence we must invert the sense of the transposition
    //
    char trans = 'N';
    if ( TwoD_distribution_ ) {
      if ( UseTranspose() ) trans = 'T' ;
//.........这里部分代码省略.........

 void operator () (const Epetra_MultiVector &x, Epetra_MultiVector &y)
 {
   y.Scale( v, x );
 };

int MatrixFree::Apply(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{

  // Use a directional derivative to compute y = Jx
  /*
   * eta = scalar perturbation
   * u = solution vector used to evaluate f
   * f = function evaluation (RHS)
   * x = vector that J is applied to
   *
   *        f(u+eta*x) - f(u)
   * Jx =   -----------------
   *               eta
   */

  // Convert X and Y from an Epetra_MultiVector to a Epetra_Vectors
  // and NOX::epetra::Vectors.  This is done so we use a consistent
  // vector space for norms and inner products.
  Teuchos::RCP<Epetra_Vector> wrappedX =
    Teuchos::rcp(new Epetra_Vector(View, X, 0));
  Teuchos::RCP<Epetra_Vector> wrappedY =
    Teuchos::rcp(new Epetra_Vector(View, Y, 0));
  NOX::Epetra::Vector nevX(wrappedX, NOX::Epetra::Vector::CreateView);
  NOX::Epetra::Vector nevY(wrappedY, NOX::Epetra::Vector::CreateView);

  // Compute perturbation constant, eta
  // Taken from LOCA v1.0 manual SAND2002-0396 p. 28 eqn. 2.43
  // eta = lambda*(lambda + 2norm(u)/2norm(x))
  double solutionNorm = 1.0;
  double vectorNorm = 1.0;

  solutionNorm = currentX.norm();
  vectorNorm = currentX.getVectorSpace()->norm(*wrappedX);

  // Make sure the norm is not zero, otherwise we can get an inf perturbation
  if (vectorNorm == 0.0) {
    //utils.out(Utils::Warning) << "Warning: NOX::Epetra::MatrixFree::Apply() - vectorNorm is zero" << std::endl;
    vectorNorm = 1.0;
    wrappedY->PutScalar(0.0);
    return 0;
  }

  // Create an extra perturbed residual vector pointer if needed
  if ( diffType == Centered )
    if ( Teuchos::is_null(fmPtr) )
      fmPtr = Teuchos::rcp(new NOX::Epetra::Vector(fo));

  double scaleFactor = 1.0;
  if ( diffType == Backward )
  scaleFactor = -1.0;

  if (computeEta) {
    if (useNewPerturbation) {
      double dotprod = currentX.getVectorSpace()->
    innerProduct(currentX.getEpetraVector(), *wrappedX);
      if (dotprod==0.0)
    dotprod = 1.0e-12;
      eta = lambda*(1.0e-12/lambda + fabs(dotprod)/(vectorNorm * vectorNorm))
    * dotprod/fabs(dotprod);
    }
    else
      eta = lambda*(lambda + solutionNorm/vectorNorm);
  }
  else
    eta = userEta;

  // Compute the perturbed RHS
  perturbX = currentX;
  Y = X;
  Y.Scale(eta);
  perturbX.update(1.0,nevY,1.0);

  if (!useGroupForComputeF)
      interface->computeF(perturbX.getEpetraVector(), fp.getEpetraVector(),
              NOX::Epetra::Interface::Required::MF_Res);
  else{
    groupPtr->setX(perturbX);
    groupPtr->computeF();
    fp = dynamic_cast<const NOX::Epetra::Vector&>
      (groupPtr->getF());
  }

  if ( diffType == Centered ) {
    Y.Scale(-2.0);
    perturbX.update(scaleFactor,nevY,1.0);
    if (!useGroupForComputeF)
      interface->computeF(perturbX.getEpetraVector(), fmPtr->getEpetraVector(),
              NOX::Epetra::Interface::Required::MF_Res);
    else{
      groupPtr->setX(perturbX);
      groupPtr->computeF();
      *fmPtr = dynamic_cast<const NOX::Epetra::Vector&>
        (groupPtr->getF());
    }
  }

  // Compute the directional derivative
  if ( diffType != Centered ) {
    nevY.update(1.0, fp, -1.0, fo, 0.0);
    nevY.scale( 1.0/(scaleFactor * eta) );
//.........这里部分代码省略.........

int FSIExactJacobian::Epetra_ExactJacobian::Apply (const Epetra_MultiVector& X, Epetra_MultiVector& Y) const
{

    LifeChrono chronoFluid, chronoSolid, chronoInterface;

    M_comm->Barrier();

    double xnorm = 0.;
    X.NormInf (&xnorm);

    Epetra_FEVector  dz (Y.Map() );

    if (M_ej->isSolid() )
    {
        std::cout << "\n ***** norm (z)= " << xnorm << std::endl << std::endl;
    }


    if ( xnorm == 0.0 )
    {
        Y.Scale (0.);
        dz.Scale (0.);
    }
    else
    {
        vector_Type const z (X,  M_ej->solidInterfaceMap(), Unique);

        M_ej->displayer().leaderPrint ( "NormInf res   " , z.normInf(), "\n" );

        //M_ej->solid().residual() *= 0.;
        //M_ej->transferInterfaceOnSolid(z, M_ej->solid().residual());

        //std::cout << "NormInf res_d " << M_ej->solid().residual().NormInf() << std::endl;

        //if (M_ej->isSolid())
        //    M_ej->solid().postProcess();

        M_ej->setLambdaFluid (z);

        //M_ej->transferInterfaceOnSolid(z, M_ej->solid().disp());

        chronoInterface.start();
        vector_Type sigmaFluidUnique (M_ej->sigmaFluid(), Unique);
        chronoInterface.stop();

        M_comm->Barrier();
        chronoFluid.start();

        if (M_ej->isFluid() )
        {

            //to be used when we correct the other lines
            if (true || ( !this->M_ej->dataFluid()->isSemiImplicit() /*|| this->M_ej->dataFluid().semiImplicit()==-1*/) )
            {
                M_ej->meshMotion().iterate (*M_ej->BCh_harmonicExtension() );
                //std::cout<<" mesh motion iterated!!!"<<std::endl;
            }

            M_ej->displayer().leaderPrint ( " norm inf dx = " , M_ej->meshMotion().disp().normInf(), "\n" );

            M_ej->solveLinearFluid();

            M_ej->transferFluidOnInterface (M_ej->fluid().residual(), sigmaFluidUnique);

            //M_ej->fluidPostProcess();
        }

        M_comm->Barrier();
        chronoFluid.stop();
        M_ej->displayer().leaderPrintMax ( "Fluid linear solution: total time : ", chronoFluid.diff() );


        chronoInterface.start();
        // M_ej->setSigmaFluid(sigmaFluidUnique);
        M_ej->setSigmaSolid (sigmaFluidUnique);


        vector_Type lambdaSolidUnique (M_ej->lambdaSolid(), Unique);
        chronoInterface.stop();

        M_comm->Barrier();
        chronoFluid.start();

        if (M_ej->isSolid() )
        {
            M_ej->solveLinearSolid();
            M_ej->transferSolidOnInterface (M_ej->solid().displacement(), lambdaSolidUnique);
        }

        M_comm->Barrier();
        chronoSolid.stop();
        M_ej->displayer().leaderPrintMax ( "Solid linear solution: total time : " , chronoSolid.diff() );

        chronoInterface.start();
        M_ej->setLambdaSolid (lambdaSolidUnique);

        chronoInterface.stop();
        M_ej->displayer().leaderPrintMax ( "Interface linear transfer: total time : " , chronoInterface.diffCumul() );

        dz = lambdaSolidUnique.epetraVector();
//.........这里部分代码省略.........