C++ Epetra_Vector类代码示例

 void observeSolution(
   const Epetra_Vector& solution)
   {
     double norm; solution.Norm2(&norm);
     if (solution.Comm().MyPID()==0)
       std::cout << "ObserveSolution: Norm = " << norm << std::endl;
   }

int power_method(bool TransA, Epetra_RowMatrix& A, Epetra_Vector& q, Epetra_Vector& z0,
		 Epetra_Vector& resid, double* lambda, int niters, double tolerance, bool verbose)
{
	
  // Fill z with random Numbers
  Epetra_Vector z(z0);
	
  // variable needed for iteration
  double normz, residual;

  int ierr = 1;
	
  for(int iter = 0; iter < niters; iter++) {
		z.Norm2(&normz); // Compute 2-norm of z
		q.Scale(1.0/normz, z);
		A.Multiply(TransA, q, z); // Compute z = A*q // SEGFAULT HAPPENS HERE
		q.Dot(z, lambda); // Approximate maximum eigenvaluE
		if(iter%100==0 || iter+1==niters) {
			resid.Update(1.0, z, -(*lambda), q, 0.0); // Compute A*q - lambda*q
			resid.Norm2(&residual);
			if(verbose) cout << "Iter = " << iter << "  Lambda = " << *lambda
											 << "  Residual of A*q - lambda*q = " << residual << endl;
		}
		if(residual < tolerance) {
			ierr = 0;
			break;
		}
	}
  return(ierr);
}

//=============================================================================
double Epetra_MsrMatrix::NormOne() const {

  if (NormOne_>-1.0) return(NormOne_);

  if (!Filled()) EPETRA_CHK_ERR(-1); // Matrix must be filled.

  Epetra_Vector * x = new Epetra_Vector(RowMatrixRowMap()); // Need temp vector for column sums
  
  Epetra_Vector * xp = 0;
  Epetra_Vector * x_tmp = 0;
  

  // If we have a non-trivial importer, we must export elements that are permuted or belong to other processors
  if (RowMatrixImporter()!=0) {
    x_tmp = new Epetra_Vector(RowMatrixColMap()); // Create temporary import vector if needed
    xp = x_tmp;
  }
  int i, j;

  for (i=0; i < NumMyCols_; i++) (*xp)[i] = 0.0;

  for (i=0; i < NumMyRows_; i++) {
    int NumEntries = GetRow(i);
    for (j=0; j < NumEntries; j++) (*xp)[Indices_[j]] += fabs(Values_[j]);
  }
  if (RowMatrixImporter()!=0) x->Export(*x_tmp, *RowMatrixImporter(), Add); // Fill x with Values from temp vector
  x->MaxValue(&NormOne_); // Find max
  if (x_tmp!=0) delete x_tmp;
  delete x;
  UpdateFlops(NumGlobalNonzeros());
  return(NormOne_);
}

void
Albany::SolutionMaxValueResponseFunction::
computeMaxValue(const Epetra_Vector& x, double& global_max, int& global_index)
{
  double my_max = -Epetra_MaxDouble;
  int my_index = -1, index;
  
  // Loop over nodes to find max value for equation eq
  int num_my_nodes = x.MyLength() / neq;
  for (int node=0; node<num_my_nodes; node++) {
    if (interleavedOrdering)  index = node*neq+eq;
    else                      index = node + eq*num_my_nodes;
    if (x[index] > my_max) {
      my_max = x[index];
      my_index = index;
    }
  }

  // Get max value across all proc's
  x.Comm().MaxAll(&my_max, &global_max, 1);

  // Compute min of all global indices equal to max value
  if (my_max == global_max)
    my_index = x.Map().GID(my_index);
  else
    my_index = x.GlobalLength();
  x.Comm().MinAll(&my_index, &global_index, 1);
}

/*----------------------------------------------------------------------*
 |  restrict from this to next coarser level (public)         m.gee 3/06|
 *----------------------------------------------------------------------*/
Epetra_Vector* NLNML::NLNML_CoarseLevelNoxInterface::restrict_to_next_coarser_level(
                                                  Epetra_Vector* thisvec,
                                                  int current, int next)
{
  Epetra_Vector* xfineP = 0;
  if (current==0)
  {
    xfineP = new Epetra_Vector((*P_)[1]->RowMap(),false);
    if (thisvec->MyLength() != xfineP->MyLength() || thisvec->GlobalLength() != xfineP->GlobalLength())
    {
      cout << "**ERR**: NLNML::NLNML_CoarseLevelNoxInterface::restrict_to_next_coarser_level:\n"
           << "**ERR**: mismatch in dimension of thisvec and xfineP\n"
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
    }
    const int mylength = thisvec->MyLength();
    for (int i=0; i<mylength; i++)
      (*xfineP)[i] = (*thisvec)[i];
  }
  else
  {
    xfineP = thisvec;
  }
  Epetra_Vector* cvec = new Epetra_Vector((*P_)[next]->OperatorDomainMap(),false);
  (*P_)[next]->Multiply(true,*xfineP,*cvec);
  if (current==0) delete xfineP;
  return cvec;
}

/*!
 * \brief Solve.
 */
void JacobiSolver::iterate( const int max_iters, const double tolerance )
{
    // Extract the linear problem.
    Epetra_CrsMatrix *A = 
	dynamic_cast<Epetra_CrsMatrix*>( d_linear_problem->GetMatrix() );
    Epetra_Vector *x = 
	dynamic_cast<Epetra_Vector*>( d_linear_problem->GetLHS() );
    const Epetra_Vector *b = 
	dynamic_cast<Epetra_Vector*>( d_linear_problem->GetRHS() );

    // Setup the residual.
    Epetra_Map row_map = A->RowMap();
    Epetra_Vector residual( row_map );

    // Iterate.
    Epetra_CrsMatrix H = buildH( A );
    Epetra_Vector temp_vec( row_map );
    d_num_iters = 0;
    double residual_norm = 1.0;
    double b_norm;
    b->NormInf( &b_norm );
    double conv_crit = b_norm*tolerance;
    while ( residual_norm > conv_crit && d_num_iters < max_iters )
    {
	H.Apply( *x, temp_vec );
	x->Update( 1.0, temp_vec, 1.0, *b, 0.0 );

	A->Apply( *x, temp_vec );
	residual.Update( 1.0, *b, -1.0, temp_vec, 0.0 );

	residual.NormInf( &residual_norm );
	++d_num_iters;
    }
}

/*----------------------------------------------------------------------*
 |  evaluate nonlinear function (public, derived)             m.gee 3/06|
 *----------------------------------------------------------------------*/
bool NLNML::NLNML_CoarseLevelNoxInterface::computeF(
                                 const Epetra_Vector& x, Epetra_Vector& F,
               const FillType fillFlag)
{
  bool err;
  if (!Level())
  {
    err = fineinterface_->computeF(x,F,fillFlag);
    if (err==false)
    {
      cout << "**ERR**: NLNML::NLNML_CoarseLevelNoxInterface::computeF:\n"
           << "**ERR**: call to fine-userinterface returned false on level " << level_ << "\n"
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
    }
  }
  else
  {
    RefCountPtr<Epetra_Vector> Ffine = rcp(new Epetra_Vector(fineinterface_->getGraph()->RowMap(),false));
    RefCountPtr<Epetra_Vector> xfine = rcp(prolong_this_to_fine(x));
    err = fineinterface_->computeF(*xfine,*Ffine,fillFlag);
    if (err==false)
    {
      cout << "**ERR**: NLNML::NLNML_CoarseLevelNoxInterface::computeF:\n"
           << "**ERR**: call to fine-userinterface returned false on level " << level_ << "\n"
           << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
    }
    RefCountPtr<Epetra_Vector> Fcoarse = rcp(restrict_fine_to_this(*Ffine));
    F.Scale(1.0,*Fcoarse);
  }
  if (isFAS())
    F.Update(-1.0,*fxbar_,1.0,*fbar_,1.0);
  return err;
}

/*----------------------------------------------------------------------*
 |  prolongate from this level to fine (public)               m.gee 3/06|
 *----------------------------------------------------------------------*/
Epetra_Vector* NLNML::NLNML_CoarseLevelNoxInterface::prolong_this_to_fine(
                                                const Epetra_Vector& xcoarse)
{
  if (!Level())
    return new Epetra_Vector(xcoarse);
  else
  {
    Epetra_Vector* cvec = const_cast<Epetra_Vector*>(&xcoarse);
    for (int i=Level(); i>0; i--)
    {
       // allocate a vector matching level i-1
       Epetra_Vector* fvec = wvec_[i-1].get();
       // multiply
       (*P_)[i]->Multiply(false,*cvec,*fvec);
       cvec = fvec;
    }
    // Note that the GIDs in cvec do NOT match those of the fineinterface as
    // they match the P_[1]->RangeMap.
    // The LIDs match, so we have to copy cvec to xfine_fineinterface
    // using LIDs
    Epetra_Vector* xfine_fineinterface = new Epetra_Vector(fineinterface_->getGraph()->RowMap(),false);
    if (cvec->MyLength() != xfine_fineinterface->MyLength() ||
        cvec->GlobalLength() != xfine_fineinterface->GlobalLength())
    {
        cout << "**ERR**: NLNML::NLNML_CoarseLevelNoxInterface::prolong_this_to_fine:\n"
             << "**ERR**: mismatch in dimension of cvec and xfine_fineinterface\n"
             << "**ERR**: file/line: " << __FILE__ << "/" << __LINE__ << "\n"; throw -1;
    }
    const int mylength = cvec->MyLength();
    for (int i=0; i<mylength; i++)
       (*xfine_fineinterface)[i] = (*cvec)[i];
    return xfine_fineinterface;
  }
}

void
Albany::SolutionTwoNormResponseFunction::
evaluateGradient(const double current_time,
		 const Epetra_Vector* xdot,
		 const Epetra_Vector& x,
		 const Teuchos::Array<ParamVec>& p,
		 ParamVec* deriv_p,
		 Epetra_Vector* g,
		 Epetra_MultiVector* dg_dx,
		 Epetra_MultiVector* dg_dxdot,
		 Epetra_MultiVector* dg_dp)
{

  // Evaluate response g
  if (g != NULL)
    x.Norm2(&(*g)[0]);

  // Evaluate dg/dx
  if (dg_dx != NULL) {
    double nrm;
    if (g != NULL)
      nrm = (*g)[0];
    else
      x.Norm2(&nrm);
    dg_dx->Scale(1.0/nrm,x);
  }

  // Evaluate dg/dxdot
  if (dg_dxdot != NULL)
    dg_dxdot->PutScalar(0.0);

  // Evaluate dg/dp
  if (dg_dp != NULL)
    dg_dp->PutScalar(0.0);
}

// B here is the "reduced" matrix.  Square matrices w/ Row=Domain=Range only.
double test_with_matvec_reduced_maps(const Epetra_CrsMatrix &A, const Epetra_CrsMatrix &B, const Epetra_Map & Bfullmap){
  const Epetra_Map & Amap  = A.DomainMap();
  Epetra_Vector Xa(Amap), Ya(Amap), Diff(Amap);
  const Epetra_Map *Bmap  = Bfullmap.NumMyElements() > 0 ? &B.DomainMap() : 0;
  Epetra_Vector *Xb = Bmap ? new Epetra_Vector(*Bmap) : 0;
  Epetra_Vector *Yb = Bmap ? new Epetra_Vector(*Bmap) : 0;

  Epetra_Vector Xb_alias(View,Bfullmap, Bmap ? Xb->Values(): 0);
  Epetra_Vector Yb_alias(View,Bfullmap, Bmap ? Yb->Values(): 0);

  Epetra_Import Ximport(Bfullmap,Amap);

  // Set the input vector
  Xa.SetSeed(24601);
  Xa.Random();
  Xb_alias.Import(Xa,Ximport,Insert);

  // Do the multiplies
  A.Apply(Xa,Ya);
  if(Bmap) B.Apply(*Xb,*Yb);

  // Check solution
  Epetra_Import Yimport(Amap,Bfullmap);
  Diff.Import(Yb_alias,Yimport,Insert);


  Diff.Update(-1.0,Ya,1.0);
  double norm;
  Diff.Norm2(&norm);

  delete Xb; delete Yb;
  return norm;
}

void
Albany::SolutionAverageResponseFunction::
evaluateGradient(const double current_time,
		 const Epetra_Vector* xdot,
		 const Epetra_Vector& x,
		 const Teuchos::Array<ParamVec>& p,
		 ParamVec* deriv_p,
		 Epetra_Vector* g,
		 Epetra_MultiVector* dg_dx,
		 Epetra_MultiVector* dg_dxdot,
		 Epetra_MultiVector* dg_dp)
{

  // Evaluate response g
  if (g != NULL)
    x.MeanValue(&(*g)[0]);

  // Evaluate dg/dx
  if (dg_dx != NULL)
    dg_dx->PutScalar(1.0 / x.GlobalLength());

  // Evaluate dg/dxdot
  if (dg_dxdot != NULL)
    dg_dxdot->PutScalar(0.0);

  // Evaluate dg/dp
  if (dg_dp != NULL)
    dg_dp->PutScalar(0.0);
}

void
Albany::SolutionAverageResponseFunction::
evaluateTangent(const double alpha, 
		const double beta,
		const double current_time,
		bool sum_derivs,
		const Epetra_Vector* xdot,
		const Epetra_Vector& x,
		const Teuchos::Array<ParamVec>& p,
		ParamVec* deriv_p,
		const Epetra_MultiVector* Vxdot,
		const Epetra_MultiVector* Vx,
		const Epetra_MultiVector* Vp,
		Epetra_Vector* g,
		Epetra_MultiVector* gx,
		Epetra_MultiVector* gp)
{
  // Evaluate response g
  if (g != NULL)
    x.MeanValue(&(*g)[0]);

  // Evaluate tangent of g = dg/dx*Vx + dg/dxdot*Vxdot + dg/dp*Vp
  // If Vx == NULL, Vx is the identity
  if (gx != NULL) {
    if (Vx != NULL)
      for (int j=0; j<Vx->NumVectors(); j++)
	(*Vx)(j)->MeanValue(&(*gx)[j][0]);
    else
      gx->PutScalar(1.0/x.GlobalLength());
    gx->Scale(alpha);
  }
  if (gp != NULL)
    gp->PutScalar(0.0);
}

 double AztecOO_StatusTestResNorm::ComputeNorm(const Epetra_Vector & vec, NormType typeofnorm) {

   double result = 0.0;
   if (typeofnorm==TwoNorm) vec.Norm2(&result);
   else if (typeofnorm==OneNorm) vec.Norm1(&result);
   else vec.NormInf(&result);
   return(result);
 }

 EpetraVector<Scalar>::EpetraVector(const Epetra_Vector &v)
 {
   this->vec = (Epetra_Vector *)&v;
   this->vec_im = nullptr;
   this->std_map = (Epetra_BlockMap *)&v.Map();
   this->size = v.MyLength();
   this->owner = false;
 }

//----------------------------------------------------------------------------
//
// Transfer solution between meshes.
//
void
AAdapt::CopyRemesh::
solutionTransfer(const Epetra_Vector& oldSolution,
                 Epetra_Vector& newSolution) {

  TEUCHOS_TEST_FOR_EXCEPT(oldSolution.MyLength() != newSolution.MyLength());
  newSolution = oldSolution;
}