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

Epetra_CrsMatrix *
NOX::Epetra::DebugTools::compute_matrix_using_operator( const Epetra_Operator * op)
{

  const Epetra_Map & rowMap  =     op->OperatorDomainMap();
  Epetra_CrsMatrix * p_mat   = new Epetra_CrsMatrix(Copy, rowMap, 0);
  Epetra_Vector *    tempVec = new Epetra_Vector(rowMap);
  Epetra_Vector *    tempRes = new Epetra_Vector(rowMap);

  // Show progress on what could be a long operation
  if( 0 == op->Comm().MyPID() )
  {
    std::cout << "****************  CREATING MATRIX FROM OPERATOR ************************ "
          << std::endl;
    std::cout << NOX::Utils::fill(72) << std::endl;
  }
  int totalPerturbations = tempVec->GlobalLength();
  int outFreq = totalPerturbations / 71;
  if( 1 > outFreq )
    outFreq = 1;

  for( int col = 0; col < tempVec->GlobalLength(); ++col )
  {
    tempVec->PutScalar(0.0);
    if( rowMap.MyGID(col) )
      (*tempVec)[col] = 1.0;

    op->Apply(*tempVec, *tempRes);

    // Fill in only the nonzero entries from the apply
    for( int row = 0; row < p_mat->NumMyRows(); ++row)
    {
      if( fabs( (*tempRes)[row] ) > 1.e-12 )
      {
    int ierr = p_mat->InsertGlobalValues( rowMap.GID(row), 1, &(*tempRes)[row], &col );
    if( ierr < 0 )
    {
          std::string msg = //"ERROR (" + ierr + ") : "
           "NOX::Epetra::DebugTools::compute_matrix_using_operator crsMatrix.ExtractGlobalRowView(...) failed for row : ";// + row;
          throw msg;
    }
      }
    }
    if( (0 == op->Comm().MyPID()) && (0 == (col % outFreq)) )
      std::cout << "-" << std::flush;
  }
  if( 0 == op->Comm().MyPID() )
    std::cout << "*" << std::endl;

  p_mat->FillComplete();

  delete tempRes;
  delete tempVec;

  return p_mat;
}

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

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

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

//=============================================================================
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 */