C++ RCP::ExtractMyRowView方法代码示例

void KfieldBC<PHAL::AlbanyTraits::Jacobian, Traits>::
evaluateFields(typename Traits::EvalData dirichletWorkset)
{

  Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
  Teuchos::RCP<Epetra_CrsMatrix> jac = dirichletWorkset.Jac;
  Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;

  RealType time = dirichletWorkset.current_time;

  const RealType j_coeff = dirichletWorkset.j_coeff;
  const std::vector<std::vector<int> >& nsNodes = 
    dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
  const std::vector<double*>& nsNodeCoords = 
    dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;

  RealType* matrixEntries;
  int*    matrixIndices;
  int     numEntries;
  RealType diag=j_coeff;
  bool fillResid = (f != Teuchos::null);

  int xlunk, ylunk; // local indicies into unknown vector
  double* coord;
  ScalarT Xval, Yval; 
  for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
  {
    xlunk = nsNodes[inode][0];
    ylunk = nsNodes[inode][1];
    coord = nsNodeCoords[inode];
    
    this->computeBCs(coord, Xval, Yval, time);
    
    // replace jac values for the X dof 
    jac->ExtractMyRowView(xlunk, numEntries, matrixEntries, matrixIndices);
    for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
    jac->ReplaceMyValues(xlunk, 1, &diag, &xlunk);

    // replace jac values for the y dof
    jac->ExtractMyRowView(ylunk, numEntries, matrixEntries, matrixIndices);
    for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
    jac->ReplaceMyValues(ylunk, 1, &diag, &ylunk);
    
    if (fillResid)
    {
      (*f)[xlunk] = ((*x)[xlunk] - Xval.val());
      (*f)[ylunk] = ((*x)[ylunk] - Yval.val());
    } 
  }
}

//.........这里部分代码省略.........
  myGlobalIndices[1] = globalIndices[2];
  broydenWorkGraph.InsertGlobalIndices( globalIndices[2], 2, &myGlobalIndices[0] );

  broydenWorkGraph.FillComplete();

  Teuchos::RCP<Epetra_CrsMatrix> broydenWorkMatrix =
    Teuchos::rcp( new Epetra_CrsMatrix( Copy, broydenWorkGraph ) );

  // Create an identity matrix
  broydenWorkVec.PutScalar(1.0);
  broydenWorkMatrix->ReplaceDiagonalValues(broydenWorkVec);

  NOX::Epetra::BroydenOperator broydenOp( noxParams, printing, broydenWorkVec, broydenWorkMatrix, true );

  broydenWorkVec[0] =  1.0;
  broydenWorkVec[1] = -1.0;
  broydenWorkVec[2] =  2.0;
  broydenOp.setStepVector( broydenWorkVec );

  broydenWorkVec[0] =  2.0;
  broydenWorkVec[1] =  1.0;
  broydenWorkVec[2] =  3.0;
  broydenOp.setYieldVector( broydenWorkVec );

  broydenOp.computeSparseBroydenUpdate();

  // Create the gold matrix for comparison
  Teuchos::RCP<Epetra_CrsMatrix> goldMatrix = Teuchos::rcp( new Epetra_CrsMatrix( Copy, broydenWorkGraph ) );

  int      numCols ;
  double * values  ;

  // Row 1 answers
  goldMatrix->ExtractMyRowView( 0, numCols, values );
  values[0] =  6.0 ;
  values[1] =  2.0 ;
  // Row 2 answers
  goldMatrix->ExtractMyRowView( 1, numCols, values );
  values[0] =  5.0 ;
  values[1] =  0.0 ;
  // Row 3 structure
  goldMatrix->ExtractMyRowView( 2, numCols, values );
  values[0] = -1.0 ;
  values[1] =  7.0 ;

  goldMatrix->Scale(0.2);

  status += tester.testCrsMatrices( broydenOp.getBroydenMatrix(), *goldMatrix, reltol, abstol,
                              "Broyden Sparse Operator Update Test" );


  // Now try a dense Broyden Update
  Epetra_CrsGraph broydenWorkGraph2( Copy, broydenRowMap, 0 );

  myGlobalIndices.resize(3);

  // All Rowsstructure
  myGlobalIndices[0] = globalIndices[0];
  myGlobalIndices[1] = globalIndices[1];
  myGlobalIndices[2] = globalIndices[2];
  broydenWorkGraph2.InsertGlobalIndices( globalIndices[0], 3, &myGlobalIndices[0] );
  broydenWorkGraph2.InsertGlobalIndices( globalIndices[1], 3, &myGlobalIndices[0] );
  broydenWorkGraph2.InsertGlobalIndices( globalIndices[2], 3, &myGlobalIndices[0] );

  broydenWorkGraph2.FillComplete();

//.........这里部分代码省略.........

                int basisId = basisIdMap[basis];

                if (checkApplyBC_)
                    if (!applyBC_[fieldIndex](worksetCellIndex,basisId))
                        continue;

                // zero out matrix row
                for(int blockColIndex=0; blockColIndex<numFieldBlocks; blockColIndex++) {
                    int start = blockOffsets[blockColIndex];
                    int end = blockOffsets[blockColIndex+1];

                    if(end-start<=0)
                        continue;

                    // check hash table for jacobian sub block
                    std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
                    Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];

                    // if you didn't find one before, add it to the hash table
                    if(subJac==Teuchos::null) {
                        Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac->getNonconstBlock(blockIndex.first,blockIndex.second);

                        // block operator is null, don't do anything (it is excluded)
                        if(Teuchos::is_null(tOp))
                            continue;

                        Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
                        subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
                        jacEpetraBlocks[blockIndex] = subJac;
                    }

                    int numEntries = 0;
                    int * rowIndices = 0;
                    double * rowValues = 0;

                    subJac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);

                    for(int i=0; i<numEntries; i++)
                        rowValues[i] = 0.0;
                }

                const ScalarT scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);

                local_r[lid] = scatterField.val();
                local_dc[lid] = 1.0; // mark row as dirichlet

                // loop over the sensitivity indices: all DOFs on a cell
                std::vector<double> jacRow(scatterField.size(),0.0);

                for(int sensIndex=0; sensIndex<scatterField.size(); ++sensIndex)
                    jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
                TEUCHOS_ASSERT(jacRow.size()==GIDs.size());

                for(int blockColIndex=0; blockColIndex<numFieldBlocks; blockColIndex++) {
                    int start = blockOffsets[blockColIndex];
                    int end = blockOffsets[blockColIndex+1];

                    if(end-start<=0)
                        continue;

                    // check hash table for jacobian sub block
                    std::pair<int,int> blockIndex = std::make_pair(blockRowIndex,blockColIndex);
                    Teuchos::RCP<Epetra_CrsMatrix> subJac = jacEpetraBlocks[blockIndex];

                    // if you didn't find one before, add it to the hash table
                    if(subJac==Teuchos::null) {
                        Teuchos::RCP<Thyra::LinearOpBase<double> > tOp = Jac->getNonconstBlock(blockIndex.first,blockIndex.second);

                        // block operator is null, don't do anything (it is excluded)
                        if(Teuchos::is_null(tOp))
                            continue;

                        Teuchos::RCP<Epetra_Operator> eOp = Thyra::get_Epetra_Operator(*tOp);
                        subJac = rcp_dynamic_cast<Epetra_CrsMatrix>(eOp,true);
                        jacEpetraBlocks[blockIndex] = subJac;
                    }

                    // Sum Jacobian
                    int err = subJac->ReplaceMyValues(lid, end-start, &jacRow[start],&LIDs[start]);
                    if(err!=0) {
                        std::stringstream ss;
                        ss << "Failed inserting row: " << GIDs[offset].second << " (" << lid << "): ";
                        for(int i=start; i<end; i++)
                            ss << GIDs[i].second << " (" << LIDs[i] << ") ";
                        ss << std::endl;
                        ss << "Into block " << blockRowIndex << ", " << blockColIndex << std::endl;

                        ss << "scatter field = ";
                        scatterFields_[fieldIndex].print(ss);
                        ss << std::endl;

                        TEUCHOS_TEST_FOR_EXCEPTION(err!=0,std::runtime_error,ss.str());
                    }

                }
            }
        }
    }
}

//.........这里部分代码省略.........

  Isorropia::Epetra::CostDescriber costs;

  Teuchos::RCP<Epetra_Vector> vptr;

  Teuchos::RCP<Epetra_CrsMatrix> eptr;

  Teuchos::RCP<Epetra_Vector> hyperEdgeWeights;

  if (edgeWeightType != NO_APPLICATION_SUPPLIED_WEIGHTS){

    if (partitioningType == GRAPH_PARTITIONING){

      // Create graph edge weights.

      eptr = Teuchos::rcp(new Epetra_CrsMatrix(*matrix));

      if (vertexWeightType == SUPPLY_EQUAL_WEIGHTS){
	eptr->PutScalar(1.0);   // set all nonzeros to 1.0
      }
      else{
	int maxRowSize = eptr->MaxNumEntries();
	double *newVal = NULL;
	if (maxRowSize > 0){
	  newVal = new double [maxRowSize];
	  for (int j=0; j<maxRowSize; j++){
	    newVal[j] = localProc + 1 + j;
	  }
	}
	int numEntries;
	int *idx;
	double *val;
	for (int i=0; i<nMyRows; i++){
	  rc = eptr->ExtractMyRowView(i, numEntries, val, idx);
	  for (int j=0; j<numEntries; j++){
	    val[j] = newVal[j];
	  }
	}
	if (newVal) delete [] newVal;
      }

      eptr->FillComplete(sourceDomainMap, sourceRangeMap);

      costs.setGraphEdgeWeights(eptr);
    }
    else{
      // Create hyperedge weights.  (Note that the list of hyperedges that a
      // process provides weights for has no relation to the columns
      // that it has non-zeroes for, or the rows that is has.  Hypergraphs
      // in general are not square.  Also more than one process can provide
      // a weight for the same edge.  Zoltan combines the weights according
      // to the value of the PHG_EDGE_WEIGHT_OPERATION parameter.  The default
      // for this parameter is to use the maximum edge weight provided by any
      // process for a given hyperedge.)

      Epetra_Map hyperEdgeMap(numCols, base, Comm);

      hyperEdgeWeights = Teuchos::rcp(new Epetra_Vector(hyperEdgeMap));

      int *edgeGIDs = NULL;
      double *weights = NULL;
      int numHEweights = hyperEdgeMap.NumMyElements();

      if (numHEweights){
	edgeGIDs = new int [numHEweights];
	weights = new double [numHEweights];

void panzer::ScatterDirichletResidual_Epetra<panzer::Traits::Jacobian, Traits,LO,GO>::
evaluateFields(typename Traits::EvalData workset)
{ 
   std::vector<GO> GIDs;
   std::vector<int> LIDs;
 
   // for convenience pull out some objects from workset
   std::string blockId = workset.block_id;
   const std::vector<std::size_t> & localCellIds = workset.cell_local_ids;

   Teuchos::RCP<const EpetraLinearObjContainer> epetraContainer = epetraContainer_;
   TEUCHOS_ASSERT(epetraContainer!=Teuchos::null);
   Teuchos::RCP<Epetra_Vector> r = epetraContainer->get_f(); 
   Teuchos::RCP<Epetra_CrsMatrix> Jac = epetraContainer->get_A();

   // NOTE: A reordering of these loops will likely improve performance
   //       The "getGIDFieldOffsets may be expensive.  However the
   //       "getElementGIDs" can be cheaper. However the lookup for LIDs
   //       may be more expensive!

   // scatter operation for each cell in workset
   for(std::size_t worksetCellIndex=0;worksetCellIndex<localCellIds.size();++worksetCellIndex) {
      std::size_t cellLocalId = localCellIds[worksetCellIndex];

      globalIndexer_->getElementGIDs(cellLocalId,GIDs); 

      if(r!=Teuchos::null) {
         // caculate the local IDs for this element
         LIDs.resize(GIDs.size());
         for(std::size_t i=0;i<GIDs.size();i++)
            LIDs[i] = r->Map().LID(GIDs[i]);
      }
      else {
         // caculate the local IDs for this element
         LIDs.resize(GIDs.size());
         for(std::size_t i=0;i<GIDs.size();i++)
            LIDs[i] = Jac->RowMap().LID(GIDs[i]);
      }

      // loop over each field to be scattered
      for(std::size_t fieldIndex = 0; fieldIndex < scatterFields_.size(); fieldIndex++) {
         int fieldNum = fieldIds_[fieldIndex];
   
         // this call "should" get the right ordering accordint to the Intrepid basis
         const std::pair<std::vector<int>,std::vector<int> > & indicePair 
               = globalIndexer_->getGIDFieldOffsets_closure(blockId,fieldNum, side_subcell_dim_, local_side_id_);
         const std::vector<int> & elmtOffset = indicePair.first;
         const std::vector<int> & basisIdMap = indicePair.second;
   
         // loop over basis functions
         for(std::size_t basis=0;basis<elmtOffset.size();basis++) {
            int offset = elmtOffset[basis];
            int lid = LIDs[offset];
            if(lid<0) // not on this processor
               continue;

            // zero out matrix row
            {
               int numEntries = 0;
               int * rowIndices = 0;
               double * rowValues = 0;

               Jac->ExtractMyRowView(lid,numEntries,rowValues,rowIndices);

               for(int i=0;i<numEntries;i++) {
                  if(preserveDiagonal_) {
                    if(lid!=rowIndices[i])
                      rowValues[i] = 0.0;
                  }
                  else
                    rowValues[i] = 0.0;
               }
            }
 
            int basisId = basisIdMap[basis];
            int gid = GIDs[offset];
            const ScalarT & scatterField = (scatterFields_[fieldIndex])(worksetCellIndex,basisId);
          
            if(r!=Teuchos::null) 
              (*r)[lid] = scatterField.val();
            if(dirichletCounter_!=Teuchos::null)
              (*dirichletCounter_)[lid] = 1.0; // mark row as dirichlet
    
            // loop over the sensitivity indices: all DOFs on a cell
            std::vector<double> jacRow(scatterField.size(),0.0);
    
            if(!preserveDiagonal_) {
              // this is the default case
              for(int sensIndex=0;sensIndex<scatterField.size();++sensIndex)
                 jacRow[sensIndex] = scatterField.fastAccessDx(sensIndex);
              TEUCHOS_ASSERT(jacRow.size()==GIDs.size());
    
              int err = Jac->ReplaceGlobalValues(gid, scatterField.size(), &jacRow[0],&GIDs[0]);
              TEUCHOS_ASSERT(err==0); 
            }
         }
      }
   }
}

// helper routines
bool SplitMatrix2x2(Teuchos::RCP<const Epetra_CrsMatrix> A,
    const Epetra_Map& A11rowmap,
    const Epetra_Map& A22rowmap,
    Teuchos::RCP<Epetra_CrsMatrix>& A11,
    Teuchos::RCP<Epetra_CrsMatrix>& A12,
    Teuchos::RCP<Epetra_CrsMatrix>& A21,
    Teuchos::RCP<Epetra_CrsMatrix>& A22)
{
  if (A==Teuchos::null)
  {
    std::cout << "ERROR: SplitMatrix2x2: A==null on entry" << std::endl;
    return false;
  }

  const Epetra_Comm& Comm   = A->Comm();
  const Epetra_Map&  A22map = A22rowmap;
  const Epetra_Map&  A11map = A11rowmap;

  //----------------------------- create a parallel redundant map of A22map
  std::map<int,int> a22gmap;
  {
    std::vector<int> a22global(A22map.NumGlobalElements());
    int count=0;
    for (int proc=0; proc<Comm.NumProc(); ++proc)
    {
      int length = 0;
      if (proc==Comm.MyPID())
      {
        for (int i=0; i<A22map.NumMyElements(); ++i)
        {
          a22global[count+length] = A22map.GID(i);
          ++length;
        }
      }
      Comm.Broadcast(&length,1,proc);
      Comm.Broadcast(&a22global[count],length,proc);
      count += length;
    }
    if (count != A22map.NumGlobalElements())
    {
      std::cout << "ERROR SplitMatrix2x2: mismatch in dimensions" << std::endl;
      return false;
    }

    // create the map
    for (int i=0; i<count; ++i)
      a22gmap[a22global[i]] = 1;
    a22global.clear();
  }

  //--------------------------------------------------- create matrix A22
  A22 = Teuchos::rcp(new Epetra_CrsMatrix(Copy,A22map,100));
  {
    std::vector<int>    a22gcindices(100);
    std::vector<double> a22values(100);
    for (int i=0; i<A->NumMyRows(); ++i)
    {
      const int grid = A->GRID(i);
      if (A22map.MyGID(grid)==false)
        continue;
      int     numentries;
      double* values;
      int*    cindices;
      int err = A->ExtractMyRowView(i,numentries,values,cindices);
      if (err)
      {
        std::cout << "ERROR: SplitMatrix2x2: A->ExtractMyRowView returned " << err << std::endl;
        return false;
      }
      if (numentries>(int)a22gcindices.size())
      {
        a22gcindices.resize(numentries);
        a22values.resize(numentries);
      }
      int count=0;
      for (int j=0; j<numentries; ++j)
      {
        const int gcid = A->ColMap().GID(cindices[j]);
        // see whether we have gcid in a22gmap
        std::map<int,int>::iterator curr = a22gmap.find(gcid);
        if (curr==a22gmap.end()) continue;
        //std::cout << gcid << " ";
        a22gcindices[count] = gcid;
        a22values[count]    = values[j];
        ++count;
      }
      //std::cout << std::endl; fflush(stdout);
      // add this filtered row to A22
      err = A22->InsertGlobalValues(grid,count,&a22values[0],&a22gcindices[0]);
      if (err<0)
      {
        std::cout << "ERROR: SplitMatrix2x2: A->InsertGlobalValues returned " << err << std::endl;
        return false;
      }

    } //for (int i=0; i<A->NumMyRows(); ++i)
    a22gcindices.clear();
    a22values.clear();
  }
//.........这里部分代码省略.........