C++ teuchos::ArrayView类代码示例

  Teuchos::RCP<const Tpetra::Map<LO,GO,Node> >
  createMeshMap (const LO& blockSize, const Tpetra::Map<LO,GO,Node>& pointMap)
  {
    typedef Teuchos::OrdinalTraits<Tpetra::global_size_t> TOT;
    typedef Tpetra::Map<LO,GO,Node> map_type;

    //calculate mesh GIDs
    Teuchos::ArrayView<const GO> pointGids = pointMap.getNodeElementList();
    Teuchos::Array<GO> meshGids;
    GO indexBase = pointMap.getIndexBase();

    // Use hash table to track whether we've encountered this GID previously.  This will happen
    // when striding through the point DOFs in a block.  It should not happen otherwise.
    // I don't use sort/make unique because I don't want to change the ordering.
    meshGids.reserve(pointGids.size());
    Tpetra::Details::HashTable<GO,int> hashTable(pointGids.size());
    for (int i=0; i<pointGids.size(); ++i) {
      GO meshGid = (pointGids[i]-indexBase) / blockSize + indexBase;
      if (hashTable.get(meshGid) == -1) {
        hashTable.add(meshGid,1);   //(key,value)
        meshGids.push_back(meshGid);
      }
    }

    Teuchos::RCP<const map_type> meshMap = Teuchos::rcp( new map_type(TOT::invalid(), meshGids(), 0, pointMap.getComm()) );
    return meshMap;

  }

    // special constructor for generating a given subblock of a strided map
    static RCP<StridedMap> Build(const RCP<const StridedMap>& map, LocalOrdinal stridedBlockId) {
      TEUCHOS_TEST_FOR_EXCEPTION(stridedBlockId < 0, Exceptions::RuntimeError,
                                 "Xpetra::StridedMapFactory::Build: constructor expects stridedBlockId > -1.");
      TEUCHOS_TEST_FOR_EXCEPTION(map->getStridedBlockId() != -1, Exceptions::RuntimeError,
                                 "Xpetra::StridedMapFactory::Build: constructor expects a full map (stridedBlockId == -1).");

      std::vector<size_t> stridingInfo = map->getStridingData();

      Teuchos::ArrayView<const GlobalOrdinal> dofGids = map->getNodeElementList();
      // std::sort(dofGids.begin(),dofGids.end()); // TODO: do we need this?

      // determine nStridedOffset
      size_t nStridedOffset = 0;
      for (int j = 0; j < map->getStridedBlockId(); j++)
        nStridedOffset += stridingInfo[j];

      size_t numMyBlockDofs = (stridingInfo[stridedBlockId] * map->getNodeNumElements()) / map->getFixedBlockSize();
      std::vector<GlobalOrdinal> subBlockDofGids(numMyBlockDofs);

      // TODO fill vector with dofs
      LocalOrdinal ind = 0;
      for (typename Teuchos::ArrayView< const GlobalOrdinal >::iterator it = dofGids.begin(); it!=dofGids.end(); ++it)
        if (map->GID2StridingBlockId(*it) == Teuchos::as<size_t>(stridedBlockId))
          subBlockDofGids[ind++] = *it;

      const Teuchos::ArrayView<const LocalOrdinal> subBlockDofGids_view(&subBlockDofGids[0],subBlockDofGids.size());

      return rcp(new StridedMap(map->lib(), Teuchos::OrdinalTraits<global_size_t>::invalid(), subBlockDofGids_view, map->getIndexBase(), stridingInfo, map->getComm(), stridedBlockId, map->getNode()));
    }

Teuchos::Array<bool> createResponseTable(
    int count,
    const std::string selectionType,
    int index,
    const Teuchos::ArrayView<const int> &list)
{
  Teuchos::Array<bool> result;

  if (count > 0) {
    if (selectionType == "All") {
      result.resize(count, true);
    } else if (selectionType == "Last") {
      result = createResponseTableFromIndex(count - 1, count);
    } else if (selectionType == "AllButLast") {
      result.reserve(count);
      result.resize(count - 1, true);
      result.push_back(false);
    } else if (selectionType == "Index") {
      result = createResponseTableFromIndex(index, count);
    } else if (selectionType == "List") {
      result.resize(count, false);
      for (Teuchos::ArrayView<const int>::const_iterator it = list.begin(), it_end = list.end(); it != it_end; ++it) {
        result.at(*it) = true;
      }
    } else {
      TEUCHOS_TEST_FOR_EXCEPT(false);
    }
  }

  return result;
}

  void
  MatrixAdapter<Matrix>::do_getCcs(const Teuchos::ArrayView<scalar_t> nzval,
				   const Teuchos::ArrayView<global_ordinal_t> rowind,
				   const Teuchos::ArrayView<typename MatrixAdapter<Matrix>::global_size_t> colptr,
				   typename MatrixAdapter<Matrix>::global_size_t& nnz,
				   const Teuchos::Ptr<const Tpetra::Map<local_ordinal_t,global_ordinal_t,node_t> > colmap,
		       EDistribution distribution,
				   EStorage_Ordering ordering,
				   col_access ca) const
  {
    using Teuchos::RCP;
    using Teuchos::ArrayView;
    using Teuchos::OrdinalTraits;
    
    RCP<const type> get_mat;
    if( *colmap == *this->col_map_ ){
      // No need to redistribute
      get_mat = rcp(this,false); // non-owning
    } else {
      get_mat = get(colmap, distribution);
    }

    // If all is well and good, then colmap == cmap
    RCP<const Tpetra::Map<scalar_t,local_ordinal_t,global_ordinal_t> > cmap = get_mat->getColMap();
    TEUCHOS_ASSERT( *colmap == *cmap );

    ArrayView<global_ordinal_t> node_elements = cmap->getNodeElementList();
    if( node_elements.size() == 0 ) return; // no more contribution
    
    typename ArrayView<global_ordinal_t>::iterator col_it, col_end;
    col_end = node_elements.end();

    size_t colptr_ind = OrdinalTraits<size_t>::zero();
    global_ordinal_t colInd = OrdinalTraits<global_ordinal_t>::zero();
    for( col_it = node_elements.begin(); col_it != col_end; ++col_it ){
      colptr[colptr_ind++] = colInd;
      size_t colNNZ = getGlobalColNNZ(*col_it);
      size_t nnzRet = 0;
      ArrayView<global_ordinal_t> rowind_view = rowind.view(colInd,colNNZ);
      ArrayView<scalar_t> nzval_view = nzval.view(colInd,colNNZ);
      getGlobalColCopy(*col_it, rowind_view, nzval_view, nnzRet);
      
      // It was suggested that instead of sorting each row's indices
      // individually, that we instead do a double-transpose at the
      // end, which would also lead to the indices being sorted.
      if( ordering == SORTED_INDICES ){
        Tpetra::sort2(rowind_view.begin(), rowind_view.end(), nzval_view.begin());
      }
      
      TEUCHOS_TEST_FOR_EXCEPTION( colNNZ != nnzRet,
			  std::runtime_error,
			  "Number of values returned different from "
                          "number of values reported");
      colInd += colNNZ;
    }
    colptr[colptr_ind] = nnz = colInd;
  }

 /// \brief Constructor.
 ///
 /// \brief matrix [in] The original input matrix.  This Container
 ///   will construct a local diagonal block from the rows given by
 ///   <tt>localRows</tt>.
 ///
 /// \param localRows [in] The set of (local) rows assigned to this
 ///   container.  <tt>localRows[i] == j</tt>, where i (from 0 to
 ///   <tt>getNumRows() - 1</tt>) indicates the Container's row, and
 ///   j indicates the local row in the calling process.  Subclasses
 ///   must always pass along these indices to the base class.
 Container (const Teuchos::RCP<const row_matrix_type>& matrix,
            const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
   inputMatrix_ (matrix),
   localRows_ (localRows.begin (), localRows.end ())
 {
   TEUCHOS_TEST_FOR_EXCEPTION(
     matrix.is_null (), std::invalid_argument, "Ifpack2::Container: "
     "The constructor's input matrix must be non-null.");
 }

  void
  RowGraph<LocalOrdinal,GlobalOrdinal,Node>::
  pack (const Teuchos::ArrayView<const LocalOrdinal>& exportLIDs,
	Teuchos::Array<GlobalOrdinal>& exports,
	const Teuchos::ArrayView<size_t>& numPacketsPerLID,
	size_t& constantNumPackets,
	Distributor& distor) const
  {
    using Teuchos::Array;
    typedef LocalOrdinal LO;
    typedef GlobalOrdinal GO;
    typedef Map<LO, GO, Node> map_type;
    const char tfecfFuncName[] = "packAndPrepare";
    (void) distor; // forestall "unused argument" compiler warning

    TEUCHOS_TEST_FOR_EXCEPTION_CLASS_FUNC(
      exportLIDs.size() != numPacketsPerLID.size(), std::runtime_error,
      ": exportLIDs and numPacketsPerLID must have the same size.");

    const map_type& srcMap = * (this->getRowMap ());
    constantNumPackets = 0;

    // Set numPacketsPerLID[i] to the number of entries owned by the
    // calling process in (local) row exportLIDs[i] of the graph, that
    // the caller wants us to send out.  Compute the total number of
    // packets (that is, entries) owned by this process in all the
    // rows that the caller wants us to send out.
    size_t totalNumPackets = 0;
    Array<GO> row;
    for (LO i = 0; i < exportLIDs.size (); ++i) {
      const GO GID = srcMap.getGlobalElement (exportLIDs[i]);
      size_t row_length = this->getNumEntriesInGlobalRow (GID);
      numPacketsPerLID[i] = row_length;
      totalNumPackets += row_length;
    }

    exports.resize (totalNumPackets);

    // Loop again over the rows to export, and pack rows of indices
    // into the output buffer.
    size_t exportsOffset = 0;
    for (LO i = 0; i < exportLIDs.size (); ++i) {
      const GO GID = srcMap.getGlobalElement (exportLIDs[i]);
      size_t row_length = this->getNumEntriesInGlobalRow (GID);
      row.resize (row_length);
      size_t check_row_length = 0;
      this->getGlobalRowCopy (GID, row (), check_row_length);
      typename Array<GO>::const_iterator row_iter = row.begin();
      typename Array<GO>::const_iterator row_end = row.end();
      size_t j = 0;
      for (; row_iter != row_end; ++row_iter, ++j) {
        exports[exportsOffset+j] = *row_iter;
      }
      exportsOffset += row.size ();
    }
  }

void
DTKInterpolationAdapter::update_variable_values(std::string var_name, Teuchos::ArrayView<GlobalOrdinal> missed_points)
{
  MPI_Comm old_comm = Moose::swapLibMeshComm(*comm->getRawMpiComm());

  System * sys = find_sys(var_name);
  unsigned int var_num = sys->variable_number(var_name);

  bool is_nodal = sys->variable_type(var_num).family == LAGRANGE;

  Teuchos::RCP<FieldContainerType> values = values_to_fill[var_name]->field();

  // Create a vector containing true or false for each point saying whether it was missed or not
  // We're only going to update values for points that were not missed
  std::vector<bool> missed(values->size(), false);

  for (Teuchos::ArrayView<const GlobalOrdinal>::const_iterator i=missed_points.begin();
      i != missed_points.end();
      ++i)
    missed[*i] = true;

  unsigned int i=0;
  // Loop over the values (one for each node) and assign the value of this variable at each node
  for (FieldContainerType::iterator it=values->begin(); it != values->end(); ++it)
  {
    // If this point "missed" then skip it
    if (missed[i])
    {
      i++;
      continue;
    }

    const DofObject * dof_object = NULL;

    if (is_nodal)
      dof_object = mesh.node_ptr(vertices[i]);
    else
      dof_object = mesh.elem(elements[i]);

    if (dof_object->processor_id() == mesh.processor_id())
    {
      // The 0 is for the component... this only works for LAGRANGE!
      dof_id_type dof = dof_object->dof_number(sys->number(), var_num, 0);
      sys->solution->set(dof, *it);
    }

    i++;
  }

  sys->solution->close();

  // Swap back
  Moose::swapLibMeshComm(old_comm);
}

//---------------------------------------------------------------------------//
// Map a reference point to the physical space of an entity.
void MoabEntityLocalMap::mapToPhysicalFrame( 
    const Entity& entity,
    const Teuchos::ArrayView<const double>& reference_point,
    const Teuchos::ArrayView<double>& physical_point ) const
{
    cacheEntity( entity );

    DTK_CHECK_ERROR_CODE(
	d_moab_evaluator->eval( reference_point.getRawPtr(),
				physical_point.getRawPtr() )
	);
}

CrsWrapper_GraphBuilder<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
CrsWrapper_GraphBuilder (const Teuchos::RCP<const Map<LocalOrdinal, GlobalOrdinal, Node> >& map)
 : graph_(),
   rowmap_(map),
   max_row_length_(0)
{
  Teuchos::ArrayView<const GlobalOrdinal> rows = map->getNodeElementList ();
  const LocalOrdinal numRows = static_cast<LocalOrdinal> (rows.size ());
  for (LocalOrdinal i = 0; i < numRows; ++i) {
    graph_[rows[i]] = new std::set<GlobalOrdinal>;
  }
}

Teuchos::Array< int >
computeCommIndexes(int rank,
                   const Teuchos::ArrayView< int > & commStrides)
{
  Teuchos::Array< int > result(commStrides.size());
  for (int axis = 0; axis < commStrides.size(); ++axis)
  {
    result[axis] = rank / commStrides[axis];
    rank         = rank % commStrides[axis];
  }
  return result;
}

void CUDANode::readyBuffers(Teuchos::ArrayView<Teuchos::ArrayRCP<const char> > buffers, Teuchos::ArrayView<Teuchos::ArrayRCP<char> > ncBuffers)
{
#ifdef HAVE_KOKKOS_DEBUG
  for (size_t i=0; i < buffers.size(); ++i) {
    CHECK_COMPUTE_BUFFER(buffers[i]);
  }
  for (size_t i=0; i < ncBuffers.size(); ++i) {
    CHECK_COMPUTE_BUFFER(ncBuffers[i]);
  }
#endif
  TEST_FOR_EXCEPT(true);
}

DenseContainer<MatrixType, LocalScalarType>::
DenseContainer (const Teuchos::RCP<const row_matrix_type>& matrix,
                const Teuchos::ArrayView<const local_ordinal_type>& localRows) :
    Container<MatrixType> (matrix, localRows),
    numRows_ (localRows.size ()),
    diagBlock_ (numRows_, numRows_),
    ipiv_ (numRows_, 0)
{
    using Teuchos::Array;
    using Teuchos::ArrayView;
    using Teuchos::RCP;
    using Teuchos::rcp;
    using Teuchos::toString;
    typedef Tpetra::Map<local_ordinal_type, global_ordinal_type, node_type> map_type;
    typedef typename ArrayView<const local_ordinal_type>::size_type size_type;
    TEUCHOS_TEST_FOR_EXCEPTION(
        ! matrix->hasColMap (), std::invalid_argument, "Ifpack2::DenseContainer: "
        "The constructor's input matrix must have a column Map.");

    // Check whether the input set of local row indices is correct.
    const map_type& rowMap = * (matrix->getRowMap ());
    const size_type numRows = localRows.size ();
    bool rowIndicesValid = true;
    Array<local_ordinal_type> invalidLocalRowIndices;
    for (size_type i = 0; i < numRows; ++i) {
        if (! rowMap.isNodeLocalElement (localRows[i])) {
            rowIndicesValid = false;
            invalidLocalRowIndices.push_back (localRows[i]);
            break;
        }
    }
    TEUCHOS_TEST_FOR_EXCEPTION(
        ! rowIndicesValid, std::invalid_argument, "Ifpack2::DenseContainer: "
        "On process " << rowMap.getComm ()->getRank () << " of "
        << rowMap.getComm ()->getSize () << ", in the given set of local row "
        "indices localRows = " << toString (localRows) << ", the following "
        "entries are not valid local row indices on the calling process: "
        << toString (invalidLocalRowIndices) << ".");

#ifdef HAVE_MPI
    RCP<const Teuchos::Comm<int> > localComm =
        rcp (new Teuchos::MpiComm<int> (MPI_COMM_SELF));
#else
    RCP<const Teuchos::Comm<int> > localComm =
        rcp (new Teuchos::SerialComm<int> ());
#endif // HAVE_MPI

    // FIXME (mfh 25 Aug 2013) What if the matrix's row Map has a
    // different index base than zero?
    const global_ordinal_type indexBase = 0;
    localMap_ = rcp (new map_type (numRows_, indexBase, localComm));
}

void CrsMatrixWrapper<ST>::ypAx(const Teuchos::ArrayView<ST>& y,
                                const Teuchos::ArrayView<const ST>& x) const
{
    RCP<VectorType<ST> > X = rcp(new VectorType<ST>(mat.getRowMap(), x, x.size(), 1));
    RCP<VectorType<ST> > Y = rcp(new VectorType<ST>(mat.getRowMap(), y, y.size(), 1));

    const ST alpha = Teuchos::ScalarTraits<ST>::one();
    const ST beta = Teuchos::ScalarTraits<ST>::one();

    // Y = beta*Y + alpha*A*X
    mat.apply(*X, *Y, Teuchos::NO_TRANS, alpha, beta);
    Y->get1dCopy(y, y.size());
}

void DropNegativeEntriesFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& currentLevel) const {
    FactoryMonitor m(*this, "Matrix filtering (springs)", currentLevel);

    RCP<Matrix> Ain = Get< RCP<Matrix> >(currentLevel, "A");

    LocalOrdinal nDofsPerNode = Ain->GetFixedBlockSize();

    // create new empty Operator
    Teuchos::RCP<Matrix> Aout = MatrixFactory::Build(Ain->getRowMap(), Ain->getGlobalMaxNumRowEntries(), Xpetra::StaticProfile);

    size_t numLocalRows = Ain->getNodeNumRows();
    for(size_t row=0; row<numLocalRows; row++) {
        GlobalOrdinal grid = Ain->getRowMap()->getGlobalElement(row);

        int rDofID = Teuchos::as<int>(grid % nDofsPerNode);

        // extract row information from input matrix
        Teuchos::ArrayView<const LocalOrdinal> indices;
        Teuchos::ArrayView<const Scalar> vals;
        Ain->getLocalRowView(row, indices, vals);

        // just copy all values in output
        Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
        Teuchos::ArrayRCP<Scalar>        valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());

        size_t nNonzeros = 0;
        for(size_t i=0; i<(size_t)indices.size(); i++) {
            GlobalOrdinal gcid = Ain->getColMap()->getGlobalElement(indices[i]); // global column id

            int cDofID = Teuchos::as<int>(gcid % nDofsPerNode);
            if(rDofID == cDofID && Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) >= Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero())) {
                indout [nNonzeros] = gcid;
                valout [nNonzeros] = vals[i];
                nNonzeros++;
            }
        }
        indout.resize(nNonzeros);
        valout.resize(nNonzeros);

        Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
    }

    Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());

    // copy block size information
    Aout->SetFixedBlockSize(nDofsPerNode);

    GetOStream(Statistics0, 0) << "Nonzeros in A (input): " << Ain->getGlobalNumEntries() << ", Nonzeros after filtering A: " << Aout->getGlobalNumEntries() << std::endl;

    Set(currentLevel, "A", Aout);
}

void RTOpC::apply_op_impl(
    const Teuchos::ArrayView<const ConstSubVectorView<Scalar> > &sub_vecs,
    const Teuchos::ArrayView<const SubVectorView<Scalar> > &targ_sub_vecs,
    const Teuchos::Ptr<ReductTarget> &_reduct_obj
  ) const
{

  using Teuchos::Workspace;
  Teuchos::WorkspaceStore* wss =
    Teuchos::get_default_workspace_store().get();

  const int num_vecs = sub_vecs.size();
  const int num_targ_vecs = targ_sub_vecs.size();

  RTOp_ReductTarget reduct_obj = RTOp_REDUCT_OBJ_NULL;
  if(!is_null(_reduct_obj))
    reduct_obj = (*this)(*_reduct_obj);

  int k;
  Workspace<RTOp_SubVector>        c_sub_vecs(wss,num_vecs,false);
  for( k = 0; k < num_vecs; ++k ) {
    const SubVector& v = sub_vecs[k];
    RTOp_sub_vector(v.globalOffset(),v.subDim(),v.values(),v.stride(),&c_sub_vecs[k]);
  }
  Workspace<RTOp_MutableSubVector>  c_targ_sub_vecs(wss,num_targ_vecs,false);
  for( k = 0; k < num_targ_vecs; ++k ) {
    const MutableSubVector& v = targ_sub_vecs[k];
    RTOp_mutable_sub_vector(v.globalOffset(),v.subDim(),v.values(),v.stride(),&c_targ_sub_vecs[k]);
  }

  const int err = RTOp_apply_op(
    &op_
    ,num_vecs,       num_vecs       ? &c_sub_vecs[0]      : (RTOp_SubVector*)NULL
    ,num_targ_vecs,  num_targ_vecs  ? &c_targ_sub_vecs[0] : (RTOp_MutableSubVector*)NULL
    ,reduct_obj
    );
  TEUCHOS_TEST_FOR_EXCEPTION(
    err==RTOp_ERR_INVALID_NUM_VECS, InvalidNumVecs
    ,"RTOpC::apply_op(...): Error, "
    "RTOp_apply_op(...) returned RTOp_ERR_INVALID_NUM_VECS" );
  TEUCHOS_TEST_FOR_EXCEPTION(
    err==RTOp_ERR_INVALID_NUM_TARG_VECS, InvalidNumTargVecs
    ,"RTOpC::apply_op(...): Error, "
    "RTOp_apply_op(...) returned RTOp_ERR_INVALID_NUM_TARG_VECS" );
  TEUCHOS_TEST_FOR_EXCEPTION(
    err!=0, UnknownError
    ,"RTOpC::apply_op(...): Error, "
    "RTOp_apply_op(...) returned != 0 with unknown meaning" );

}