C++ Epetra_Vector::MyLength方法代码示例

/*----------------------------------------------------------------------*
 |  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;
  }
}

//----------------------------------------------------------------------------
//
// 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;
}

int Ifpack_AnalyzeVectorElements(const Epetra_Vector& Diagonal,
                                 const bool abs, const int steps)
{

  bool verbose = (Diagonal.Comm().MyPID() == 0);
  double min_val =  DBL_MAX;
  double max_val = -DBL_MAX;

  for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
    double v = Diagonal[i];
    if (abs)
      if (v < 0) v = -v;
    if (v > max_val)
      max_val = v;
    if (v < min_val)
      min_val = v;
  }

  if (verbose) {
    cout << endl;
    Ifpack_PrintLine();
    cout << "Vector label = " << Diagonal.Label() << endl;
    cout << endl;
  }

  double delta = (max_val - min_val) / steps;
  for (int k = 0 ; k < steps ; ++k) {

    double below = delta * k + min_val;
    double above = below + delta;
    int MyBelow = 0, GlobalBelow;

    for (int i = 0 ; i < Diagonal.MyLength() ; ++i) {
      double v = Diagonal[i];
      if (v < 0) v = -v;
      if (v >= below && v < above) MyBelow++;
    }

    Diagonal.Comm().SumAll(&MyBelow, &GlobalBelow, 1);

    if (verbose) {
      printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
             below, above, GlobalBelow,
                         100.0 * GlobalBelow / Diagonal.GlobalLength64());
    }
  }

  if (verbose) {
    Ifpack_PrintLine();
    cout << endl;
  }

  return(0);
}

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;
}

 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;
 }

//=============================================================================
int Epetra_PETScAIJMatrix::RightScale(const Epetra_Vector& X) {
//
// This function scales the jth row of A by x[j].
//
  double *xptr;
  X.ExtractView(&xptr);
  Vec petscX;
# ifdef HAVE_MPI
  int ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
# else //FIXME  untested
  int ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptr,&petscX); CHKERRQ(ierr);
# endif

  MatDiagonalScale(Amat_, PETSC_NULL, petscX);

  ierr=VecDestroy(petscX); CHKERRQ(ierr);
  return(0);
} //RightScale()

bool Vector2MATLAB( const Epetra_Vector & v,
		    const Epetra_Map & Map)
{
  
  int MyPID = Map.Comm().MyPID(); 
  int NumProc = Map.Comm().NumProc();
  int MyLength = v.MyLength();
  int GlobalLength = v.GlobalLength();
  
  // print out on cout if no filename is provided

  // write on file the dimension of the matrix

  if( MyPID == 0 ) cout << "v = zeros(" << GlobalLength << ")\n";

  // get update list
  int * MyGlobalElements = Map.MyGlobalElements( );
  
  int Row;

  for( int Proc=0 ; Proc<NumProc ; ++Proc ) {

    if( MyPID == Proc ) {

      cout << "% On proc " << Proc << ": ";
      cout << MyLength << " rows of ";
      cout << GlobalLength << " elements\n";

      for( Row=0 ; Row<MyLength ; ++Row ) {
	cout << "b(" << MyGlobalElements[Row]
	     << ") = " << v[Row] << ";\n";
      }
      
      if( MyPID == NumProc-1  ) {
	cout << "% End of vector\n";
      }
      
    }
      
    Map.Comm().Barrier();
  }

  return true;

} /* Vector2MATLAB */

/******************************************************************
  Compute weight balance
******************************************************************/
int compute_balance(const Epetra_Vector &wgts, double myGoalWeight,
                      double &min, double &max, double &avg)
{
  if ((myGoalWeight < 0) || (myGoalWeight > 1.0)){
    std::cerr << "compute_balance: Goal weight should be in the range [0, 1]" << std::endl;
    return -1;
  }

  double weightTotal;
  wgts.Norm1(&weightTotal);

  double weightLocal = 0.0;
  for (int i=0; i < wgts.MyLength(); i++){
    weightLocal += wgts[i];
  }

  /* My degree of imbalance. 
   * If myGoalWeight is zero, I'm in perfect balance since I got what I wanted.
   */
  double goalWeight = myGoalWeight * weightTotal;
  double imbalance = 1.0;

  if (myGoalWeight > 0.0){
    if (weightLocal >= goalWeight)
      imbalance += (weightLocal - goalWeight) / goalWeight;
    else
      imbalance += (goalWeight - weightLocal) / goalWeight;
  }

  const Epetra_Comm &comm = wgts.Comm();

  comm.MaxAll(&imbalance, &max, 1);
  comm.MinAll(&imbalance, &min, 1);
  comm.SumAll(&imbalance, &avg, 1);

  avg /= comm.NumProc();

  return 0;
}

void AZOO_iterate(double * xsolve, double * b, 
		  int * options, double * params, 
		  double * status, int *proc_config,
		  AZ_MATRIX * Amat,
		  AZ_PRECOND *precond, struct AZ_SCALING *scaling)
{
  (void)precond;
  (void)scaling;
  bool verbose = (options[AZ_output]!=AZ_none); // Print info unless all output is turned off

  Epetra_Comm * comm;
  Epetra_BlockMap * map;
  Epetra_RowMatrix * A;
  Epetra_Vector * px;
  Epetra_Vector * pb;
  int * global_indices;

  int ierr = Aztec2Petra(proc_config, Amat, xsolve, b, comm, map, A, px, pb, &global_indices);
  if (ierr!=0) {
    cerr << "Error detected in Aztec2Petra. Value = " << ierr << endl;
    exit(1);
  }


  Epetra_LinearProblem problem(A, px, pb);

  Epetra_Vector * leftScaleVec = 0;
  Epetra_Vector * rightScaleVec = 0;
  bool doRowScaling = false;
  bool doColScaling = false;
  
  if ((options[AZ_scaling]==AZ_Jacobi) || options[AZ_scaling]==AZ_BJacobi) {
    doRowScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->ExtractDiagonalCopy(*leftScaleVec); // Extract diagonal of matrix
    leftScaleVec->Reciprocal(*leftScaleVec); // invert it
  }

  else if (options[AZ_scaling]==AZ_row_sum) {
    doRowScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->InvRowSums(*leftScaleVec);
  }
  else if (options[AZ_scaling]==AZ_sym_diag) {
    doRowScaling = true;
    doColScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->ExtractDiagonalCopy(*leftScaleVec); // Extract diagonal of matrix

    int length = leftScaleVec->MyLength();
    for (int i=0; i<length; i++) (*leftScaleVec)[i] = sqrt(fabs((*leftScaleVec)[i])); // Take its sqrt

    rightScaleVec = leftScaleVec; // symmetric, so left and right the same
    leftScaleVec->Reciprocal(*leftScaleVec); // invert it
  }
  else if (options[AZ_scaling]==AZ_sym_row_sum) {
    doRowScaling = true;
    doColScaling = true;
    leftScaleVec = new Epetra_Vector(*map);
    A->InvRowSums(*leftScaleVec);
    int length = leftScaleVec->MyLength();
    for (int i=0; i<length; i++) (*leftScaleVec)[i] = sqrt(fabs((*leftScaleVec)[i])); // Take its sqrt

    rightScaleVec = leftScaleVec; // symmetric, so left and right the same
  }
  if ((doRowScaling || doColScaling) && verbose) {
    double norminf = A->NormInf();
    double normone = A->NormOne();
    if (comm->MyPID()==0) 
      cout << "\n Inf-norm of A before scaling = " << norminf 
	   << "\n One-norm of A before scaling = " << normone<< endl << endl;
  }
  if (doRowScaling) problem.LeftScale(*leftScaleVec);
  if (doColScaling) problem.RightScale(*rightScaleVec);

  if ((doRowScaling || doColScaling) && verbose) {
    double norminf = A->NormInf();
    double normone = A->NormOne();
    if (comm->MyPID()==0) 
      cout << "\n Inf-norm of A after  scaling = " << norminf  
	   << "\n One-norm of A after  scaling = " << normone << endl << endl;
  }



  AztecOO solver(problem);

  solver.SetAllAztecParams(params); // set all AztecOO params with user-provided params
  solver.SetAllAztecOptions(options); // set all AztecOO options with user-provided options

  solver.CheckInput();
  solver.SetAztecOption(AZ_scaling, AZ_none); // Always must have scaling off
  solver.Iterate(options[AZ_max_iter], params[AZ_tol]);
  solver.GetAllAztecStatus(status);
  
  if (doColScaling) {
    rightScaleVec->Reciprocal(*rightScaleVec);
    problem.RightScale(*rightScaleVec);
  }
  if (doRowScaling) {
//.........这里部分代码省略.........

void 
HMX_PDE::computeSourceTerm(map<string, Epetra_Vector*> fields, 
                                Epetra_Vector& result)
{
  // All dependent variables should be in place before now

  // We assume all other registered dependent problems are species which
  // affect the reaction rate of this specie

  Epetra_Vector *TvecPtr = 0;
  map<string, Epetra_Vector*>::iterator iter = fields.find(tempFieldName);
  if( iter == fields.end() ) 
  {
    std::cout << "ERROR: Cannot find Temperature field \"" << tempFieldName 
         << "\" for use in computeSourceTerm for problem \""
         << getName() << std::endl;
    throw "HMX_PDE ERROR";
  }
  else
    TvecPtr = (*iter).second;

  Epetra_Vector &T = *TvecPtr;

  // If this problem is the temperature equation, don't compute a source
  // term.  This would be where a volumetric heating term would go.
  if( getName() == tempFieldName ) 
  {
    result.PutScalar(0.0);
    return;
  }
  else 
  {
    double rateK;
    
    map<string, double>::iterator requiredFieldIter;
    map<string, double>::iterator requiredFieldEnd = SrcTermExponent.end();

    for( int i = 0; i<result.MyLength(); i++ ) {

      rateK = pow(T[i],StericCoef) * PreExp_A * 
              exp( -ActEnergy / (Const_R * T[i]) );

      result[i] = 1.0;  // Start point for product

      // Loop over all required fields and contribute to product
      for( requiredFieldIter = SrcTermExponent.begin();
           requiredFieldIter != requiredFieldEnd; requiredFieldIter++) 
      {
        iter = fields.find( (*requiredFieldIter).first );
        if( iter == fields.end() ) 
        {
          std::cout << "ERROR: Cannot find required field \"" 
               << (*requiredFieldIter).first 
               << "\" for use in computeSourceTerm for problem \""
               << getName() << std::endl;
          throw "HMX_PDE ERROR";
        }
	Epetra_Vector &reqFieldVec = *((*iter).second);

        result[i] *= pow( reqFieldVec[i], (*requiredFieldIter).second );
      }
      result[i] *= rateK;
    }
  }
}

 virtual int ExtractDiagonalCopy(Epetra_Vector & Diagonal) const
 {
   for (int i = 0 ; i < Diagonal.MyLength() ; ++i)
     Diagonal[i] = 7.0;
   return(0);
 }