C++ Array::clear方法代码示例

void LibmeshAdjacencies::getLibmeshAdjacencies<libMesh::Elem,libMesh::Elem>(
    const Teuchos::Ptr<libMesh::Elem>& /*entity*/,
    Teuchos::Array<Teuchos::Ptr<libMesh::Elem> >& adjacent_entities
    ) const
{
    adjacent_entities.clear();
}

//---------------------------------------------------------------------------//
// Given an entity, get the entities of the given type that are adjacent to
// it. 
void LibmeshEntitySet::getAdjacentEntities(
    const DataTransferKit::Entity& entity,
    const int adjacent_dimension,
    Teuchos::Array<DataTransferKit::Entity>& adjacent_entities ) const
{
    int entity_topo_dim = entity.topologicalDimension();
    
    if ( 0 == entity_topo_dim )
    {
	if ( 0 == adjacent_dimension )
	{
	    adjacent_entities.clear();
	}
	else
	{
	    getAdjacentEntitiesImpl<libMesh::Node,libMesh::Elem>(
		entity, adjacent_entities );
	}
    }
    else
    {
	if ( 0 == adjacent_dimension )
	{
	    getAdjacentEntitiesImpl<libMesh::Elem,libMesh::Node>(
		entity, adjacent_entities );
	}
	else
	{
	    getAdjacentEntitiesImpl<libMesh::Elem,libMesh::Elem>(
		entity, adjacent_entities );
	}
    }
}

void epetraFromThyra(
  const Teuchos::RCP<const Epetra_Comm> &comm,
  const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > &thyraResponses,
  const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > &thyraSensitivities,
  Teuchos::Array<Teuchos::RCP<const Epetra_Vector> > &responses,
  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > > &sensitivities)
{
  responses.clear();
  responses.reserve(thyraResponses.size());
  typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<double> > > ThyraResponseArray;
  for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
      it_end = thyraResponses.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    responses.push_back(epetraVectorFromThyra(comm, *it));
  }

  sensitivities.clear();
  sensitivities.reserve(thyraSensitivities.size());
  typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<double> > > > ThyraSensitivityArray;
  for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
      it_end = thyraSensitivities.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    ThyraSensitivityArray::const_reference sens_thyra = *it;
    Teuchos::Array<Teuchos::RCP<const Epetra_MultiVector> > sens;
    sens.reserve(sens_thyra.size());
    for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
        jt_end = sens_thyra.end();
        jt != jt_end;
        ++jt) {
        sens.push_back(epetraMultiVectorFromThyra(comm, *jt));
    }
    sensitivities.push_back(sens);
  }
}

void tpetraFromThyra(
  const Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<ST> > > &thyraResponses,
  const Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<ST> > > > &thyraSensitivities,
  Teuchos::Array<Teuchos::RCP<const Tpetra_Vector> > &responses,
  Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > > &sensitivities)
{
  responses.clear();
  responses.reserve(thyraResponses.size());
  typedef Teuchos::Array<Teuchos::RCP<const Thyra::VectorBase<ST> > > ThyraResponseArray;
  for (ThyraResponseArray::const_iterator it_begin = thyraResponses.begin(),
      it_end = thyraResponses.end(),
      it = it_begin;
      it != it_end;
      ++it) {
    responses.push_back(Teuchos::nonnull(*it) ? ConverterT::getConstTpetraVector(*it) : Teuchos::null);
  }

  sensitivities.clear();
  sensitivities.reserve(thyraSensitivities.size());
  typedef Teuchos::Array<Teuchos::Array<Teuchos::RCP<const Thyra::MultiVectorBase<ST> > > > ThyraSensitivityArray;
  for (ThyraSensitivityArray::const_iterator it_begin = thyraSensitivities.begin(),
      it_end = thyraSensitivities.end(),
      it = it_begin;
      it != it_end;
      ++it) {
  ThyraSensitivityArray::const_reference sens_thyra = *it;
    Teuchos::Array<Teuchos::RCP<const Tpetra_MultiVector> > sens;
    sens.reserve(sens_thyra.size());
    for (ThyraSensitivityArray::value_type::const_iterator jt = sens_thyra.begin(),
        jt_end = sens_thyra.end();
        jt != jt_end;
        ++jt) {
        sens.push_back(Teuchos::nonnull(*jt) ? ConverterT::getConstTpetraMultiVector(*jt) : Teuchos::null);
    }
    sensitivities.push_back(sens);
  }
}

void
Tpetra::Utils::generateMatrix(const Teuchos::RCP<Teuchos::ParameterList> &plist,
                              const Teuchos::RCP<const Teuchos::Comm<int> > &comm,
                              const Teuchos::RCP<Node> &node,
                              Teuchos::RCP< Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node,LocalMatOps> > &A)
{
    typedef Teuchos::ScalarTraits<Scalar> ST;
    using Teuchos::as;
    TEUCHOS_TEST_FOR_EXCEPTION( plist == Teuchos::null, std::runtime_error,
                                "Tpetra::Utils::generateMatrix(): ParameterList is null.");
    TEUCHOS_TEST_FOR_EXCEPTION( Teuchos::isParameterType<std::string>(*plist,"mat_type") == false, std::runtime_error,
                                "Tpetra::Utils::generateMatrix(): ParameterList did not contain string parameter ""mat_type"".");
    std::string mat_type = plist->get<std::string>("mat_type");
    if (mat_type == "Lap3D") {
        // 3D Laplacian, grid is a cube with dimension gridSize x gridSize x gridSize
        const GlobalOrdinal gridSize = as<GlobalOrdinal>(plist->get<int>("gridSize",100));
        const GlobalOrdinal gS2 = gridSize*gridSize;
        const GlobalOrdinal numRows = gS2*gridSize;
        Teuchos::RCP<Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap;
        rowMap = Teuchos::rcp(new Tpetra::Map<LocalOrdinal,GlobalOrdinal,Node>(as<global_size_t>(numRows),as<GlobalOrdinal>(0),comm,GloballyDistributed,node));
        A = rcp(new Tpetra::CrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(rowMap,7,Tpetra::StaticProfile));
        // fill matrix, one row at a time
        Teuchos::Array<GlobalOrdinal> neighbors;
        Teuchos::Array<Scalar> values(7, -ST::one());
        values[0] = (Scalar)6;
        for (GlobalOrdinal r = rowMap->getMinGlobalIndex(); r <= rowMap->getMaxGlobalIndex(); ++r) {
            neighbors.clear();
            neighbors.push_back(r); // add diagonal
            GlobalOrdinal ixy, iz, ix, iy;  // (x,y,z) coords and index in xy plane
            ixy = r%gS2;
            iz = (r - ixy)/gS2;
            ix = ixy%gridSize;
            iy = (ixy - ix)/gridSize;
            //
            if ( ix != 0 )          neighbors.push_back( r-1 );
            if ( ix != gridSize-1 ) neighbors.push_back( r+1 );
            if ( iy != 0 )          neighbors.push_back( r-gridSize );
            if ( iy != gridSize-1 ) neighbors.push_back( r+gridSize );
            if ( iz != 0 )          neighbors.push_back( r-gS2 );
            if ( iz != gridSize-1 ) neighbors.push_back( r+gS2 );
            A->insertGlobalValues( r, neighbors(), values(0,neighbors.size()) );
        }
        A->fillComplete();
    }
    else {
        TEUCHOS_TEST_FOR_EXCEPTION( true, std::runtime_error,
                                    "Tpetra::Utils::generateMatrix(): ParameterList specified unsupported ""mat_type"".");
    }
}

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

	// Partition in x direction.
	if ( partitioning_type == 0 )
	{
	    for ( int i = 0; i < comm_size; ++i )
	    {
		if ( elem_coords[0] >= min + width*(comm_rank)/comm_size - 1e-6 &&
		     elem_coords[0] <= min + width*(comm_rank+1)/comm_size + 1e-6 )
		{
		    elements.push_back( *global_elem_iterator );
		}
	    }
	}

	// Partition in y direction.
	else if ( partitioning_type == 1 )
	{
	    for ( int i = 0; i < comm_size; ++i )
	    {
		if ( elem_coords[1] >= min + width*(comm_rank)/comm_size - 1e-6 &&
		     elem_coords[1] <= min + width*(comm_rank+1)/comm_size + 1e-6 )
		{
		    elements.push_back( *global_elem_iterator );
		}

	    }
	}
	else
	{
	    throw std::logic_error( "Partitioning type not supported." );
	}
    }
    Teuchos::ArrayRCP<moab::EntityHandle> elements_arcp( elements.size() );
    std::copy( elements.begin(), elements.end(), elements_arcp.begin() );
    elements.clear();
    d_comm->barrier();

    // Get the nodes.
    std::vector<moab::EntityHandle> vertices;
    error = d_moab->get_connectivity( &elements_arcp[0],
				      elements_arcp.size(),
				      vertices );
    assert( error == moab::MB_SUCCESS );
    d_vertices = Teuchos::ArrayRCP<moab::EntityHandle>( vertices.size() );
    std::copy( vertices.begin(), vertices.end(), d_vertices.begin() );
    vertices.clear();

    // Get the node coordinates.
    Teuchos::ArrayRCP<double> coords( node_dim * d_vertices.size() );
    std::vector<double> interleaved_coords( 3*d_vertices.size() );
    error = d_moab->get_coords( &d_vertices[0], d_vertices.size(),
				&interleaved_coords[0] );
    assert( error == moab::MB_SUCCESS );

    for ( int n = 0; n < (int) d_vertices.size(); ++n )
    {
	for ( int d = 0; d < (int) node_dim; ++d )
	{
	    coords[ d*d_vertices.size() + n ] =
		interleaved_coords[ n*3 + d ];
	}
    }
    interleaved_coords.clear();

    // Get the connectivity.
    int connectivity_size = elements_arcp.size() * nodes_per_element;
    Teuchos::ArrayRCP<moab::EntityHandle> connectivity( connectivity_size );
    std::vector<moab::EntityHandle> elem_conn;
    for ( int i = 0; i < (int) elements_arcp.size(); ++i )
    {
	error = d_moab->get_connectivity( &elements_arcp[i], 1, elem_conn );

	assert( error == moab::MB_SUCCESS );
	assert( elem_conn.size() ==
		Teuchos::as<std::vector<moab::EntityHandle>::size_type>(nodes_per_element) );

	for ( int n = 0; n < (int) elem_conn.size(); ++n )
	{
	    connectivity[ n*elements_arcp.size() + i ] = elem_conn[n];
	}
    }

    // Get the permutation vector.
    Teuchos::ArrayRCP<int> permutation_list( nodes_per_element );
    for ( int i = 0; i < (int) nodes_per_element; ++i )
    {
	permutation_list[i] = i;
    }

    // Create the mesh container.
    d_mesh_container = Teuchos::rcp(
	new Container( node_dim,
		       d_vertices,
		       coords,
		       getTopology(block_topology),
		       nodes_per_element,
		       elements_arcp,
		       connectivity,
		       permutation_list ) );
}

Epetra_CrsMatrix*
read_matrix_mm(const std::string& mm_file,
               const Epetra_Comm& comm)
{
  int my_proc = comm.MyPID();

  long long num_global_rows = 0;
  int nnz_per_row = 0;

  std::ifstream* infile = NULL;
  infile = new std::ifstream(mm_file.c_str());
  if (infile == NULL || !*infile) {
    throw std::runtime_error("Failed to open file "+mm_file);
  }

  std::ifstream& in = *infile;

  //first skip over the file header, which has
  //lines beginning with '%'.
  std::string line;
  do {
    getline(in, line);
  } while(line[0] == '%');

  //now get the matrix dimensions.

  int numrows, numcols, nnz;
  std::istringstream isstr(line);
  isstr >> numrows >> numcols >> nnz;

  //make sure we successfully read the three ints from that line.
  if (isstr.fail()) {
    throw std::runtime_error("Failed to parse matrix-market header.");
  }

  if (my_proc == 0) {
    num_global_rows = numrows;
    nnz_per_row = nnz/numrows;
  }

  comm.Broadcast(&num_global_rows, 1, 0);
  comm.Broadcast(&nnz_per_row, 1, 0);

  const int indexBase = 0;
  Epetra_Map rowmap(num_global_rows, indexBase, comm);

  Epetra_CrsMatrix* A = new Epetra_CrsMatrix(Copy, rowmap, nnz_per_row);

  Teuchos::Array<long long> col;
  Teuchos::Array<double> coef;

  int irow=0, icol=0;
  int g_row=-1, last_row=-1;
  double val=0;

  while(!in.eof()) {
    getline(in, line);
    std::istringstream isstr(line);
    isstr >> irow >> icol >> val;
  
    if (isstr.fail()) continue;
    if (!rowmap.MyGID(irow-1)) continue;

    g_row = irow-1;
    if (g_row != last_row) {
      if (col.size() > 0) {
        A->InsertGlobalValues(last_row, col.size(), &coef[0], &col[0] );
        col.clear();
        coef.clear();
      }
      last_row = g_row;
    }
    col.push_back(icol-1);
    coef.push_back(val);
  }

  if (col.size() > 0) {
    A->InsertGlobalValues(g_row, col.size(), &coef[0], &col[0]);
  }

  A->FillComplete();

  return A;
}