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

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

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 observeSolution(
   const Epetra_Vector& solution)
   {
     double norm; solution.Norm2(&norm);
     if (solution.Comm().MyPID()==0)
       std::cout << "ObserveSolution: Norm = " << norm << std::endl;
   }

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

bool FiniteDifferenceColoring::computeJacobian(const Epetra_Vector& x, Epetra_Operator& Jac)
{
  // First check to make sure Jac is a
  // NOX::Epetra::FiniteDifferenceColoring object
  FiniteDifferenceColoring* testMatrix =
         dynamic_cast<FiniteDifferenceColoring*>(&Jac);
  if (testMatrix == 0) {
    cout << "ERROR: NOX::Epetra::FiniteDifferenceColoring::computeJacobian() - "
	 << "Jacobian to evaluate is not a FiniteDifferenceColoring object!"
         << endl;
    throw "NOX Error";
  }

  const Epetra_BlockMap& map = fo.Map();

  // Create a timer for performance
  Epetra_Time fillTimer(x.Comm());

  // We need the Epetra_CrsMatrix inside the FiniteDifferenceColoring object
  // for the correct insertion commands.
  Epetra_CrsMatrix& jac = *testMatrix->jacobian;

  // Zero out Jacobian
  jac.PutScalar(0.0);

  // Create an extra perturbed residual vector pointer if needed
  if ( diffType == Centered )
    if ( Teuchos::is_null(fmPtr) )
      fmPtr = Teuchos::rcp(new Epetra_Vector(x));

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

  int myMin = map.MinMyGID(); // Minimum locally owned GID 
  int myMax = map.MaxMyGID(); // Maximum locally owned GID

  // We need to loop over the largest number of colors on a processor

  // Use the overlap (column-space) version of the solution
  xCol_perturb->Import(x, *rowColImporter, Insert);

  // Compute the RHS at the initial solution
  if( coloringType == NOX_SERIAL )
    computeF(*xCol_perturb, fo, NOX::Epetra::Interface::Required::FD_Res);
  else {
//    x_perturb.Export(*xCol_perturb, *rowColImporter, Insert);
    for (int i=0; i<x_perturb.MyLength(); i++)
      x_perturb[i] = (*xCol_perturb)[columnMap->LID(map.GID(i))];
    computeF(x_perturb, fo, NOX::Epetra::Interface::Required::FD_Res);
  }

  // loop over each color in the colorGraph
  list<int>::iterator allBegin = listOfAllColors.begin(),
                      allEnd   = listOfAllColors.end(),
                      allIter;

  std::map<int, int>::iterator mapEnd = colorToNumMap.end(),
                               myMapIter;
  for ( allIter = allBegin; allIter != allEnd; ++allIter ) {
    bool skipIt = true;       // Used to screen colors this proc does not have
    int k = -1;               // index in colorList of color
    myMapIter = colorToNumMap.find( *allIter );
    if( myMapIter != mapEnd ) {
      skipIt = false;
      k = (*myMapIter).second;
    }

    // Perturb the solution vector using coloring

    // ----- First create color map and assoc perturbation vectors
    int color = -1;
    if( !skipIt ) color = colorList[k];
    cMap = colorMap->GenerateMap(color);
    colorVect = new Epetra_Vector(*cMap);
    betaColorVect = new Epetra_Vector(*cMap);
    betaColorVect->PutScalar(beta);

    // ----- Fill colorVect with computed perturbation values
    // NOTE that we do the mapping ourselves here instead of using an
    // Epetra_Import object.  This is to ensure local mapping only since
    // Import/Export operations involve off-processor transfers, which we
    // wish to avoid.
    for (int i=0; i<colorVect->MyLength(); i++)
      (*colorVect)[i] = (*xCol_perturb)[columnMap->LID(cMap->GID(i))];
    colorVect->Abs(*colorVect);
    colorVect->Update(1.0, *betaColorVect, alpha);

    // ----- Map perturbation vector to original index space
    // ----- and use it
    mappedColorVect->PutScalar(0.0);
    // Here again we do the mapping ourselves to avoid off-processor data
    // transfers that would accompany use of Epetra_Import/Export objects.
    for (int i=0; i<colorVect->MyLength(); i++)
      (*mappedColorVect)[columnMap->LID(cMap->GID(i))] = (*colorVect)[i];
    xCol_perturb->Update(scaleFactor, *mappedColorVect, 1.0);

    // Compute the perturbed RHS
    if( coloringType == NOX_SERIAL )
      computeF(*xCol_perturb, fp, NOX::Epetra::Interface::Required::FD_Res);
//.........这里部分代码省略.........

bool FiniteDifferenceColoringWithUpdate::differenceProbe(const Epetra_Vector& x, Epetra_CrsMatrix& jac,const Epetra_MapColoring& colors){

  // Allocate space for perturbation, get column version of x for scaling
  Epetra_Vector xp(x);
  Epetra_Vector *xcol;
  int N=jac.NumMyRows();

  if(jac.ColMap().SameAs(x.Map()))
     xcol=const_cast<Epetra_Vector*>(&x);
  else{
    xcol=new Epetra_Vector(jac.ColMap(),true);//zeros out by default
    xcol->Import(x,*jac.Importer(),InsertAdd);
  }

  // Counters for probing diagnostics
  double tmp,probing_error_lower_bound=0.0,jc_norm=0.0;

  // Grab coloring info (being very careful to ignore color 0)
  int Ncolors=colors.MaxNumColors()+1;
  int num_c0_global,num_c0_local=colors.NumElementsWithColor(0);
  colors.Comm().MaxAll(&num_c0_local,&num_c0_global,1);
  if(num_c0_global>0) Ncolors--;

  if(Ncolors==0) return false;

  // Pointers for Matrix Info
  int entries, *indices;
  double *values;

  // NTS: Fix me
  if ( diffType == Centered ) exit(1);

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

  // Compute RHS at initial solution
  computeF(x,fo,NOX::Epetra::Interface::Required::FD_Res);

  /* Probing, vector by vector since computeF does not have a MultiVector interface */
  // Assume that anything with Color 0 gets ignored.
  for(int j=1;j<Ncolors;j++){
    xp=x;
    for(int i=0;i<N;i++){
      if(colors[i]==j)
    xp[i] += scaleFactor*(alpha*abs(x[i])+beta);
    }

    computeF(xp, fp, NOX::Epetra::Interface::Required::FD_Res);

    // Do the subtraction to estimate the Jacobian (w/o including step length)
    Jc.Update(1.0, fp, -1.0, fo, 0.0);

    // Relative error in probing
     if(use_probing_diags){
       Jc.Norm2(&tmp);
       jc_norm+=tmp*tmp;
     }

    for(int i=0;i<N;i++){
      // Skip for uncolored row/columns, else update entries
      if(colors[i]==0) continue;

      jac.ExtractMyRowView(i,entries,values,indices);
      for(int k=0;k<jac.NumMyEntries(i);k++){
    if(colors[indices[k]]==j){
      values[k]=Jc[i] / (scaleFactor*(alpha*abs((*xcol)[indices[k]])+beta));
      // If probing diagnostics are on, zero out the entries as they are used
      if(use_probing_diags) Jc[i]=0.0;
      break;// Only one value per row...
    }
      }
    }
    if(use_probing_diags){
      Jc.Norm2(&tmp);
      probing_error_lower_bound+=tmp*tmp;
    }
  }

  // If diagnostics are requested, output Frobenius norm lower bound
  if(use_probing_diags && !x.Comm().MyPID()) printf("Probing Error Lower Bound (Frobenius) abs = %6.4e rel = %6.4e\n",sqrt(probing_error_lower_bound),sqrt(probing_error_lower_bound)/sqrt(jc_norm));

  // Cleanup
  if(!jac.ColMap().SameAs(x.Map()))
    delete xcol;

  return true;
}