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

//---------------------------------------------------------------------------//
TEUCHOS_UNIT_TEST( EpetraPointJacobiPreconditioner, tridiag_matrix )
{
    typedef Epetra_RowMatrix MatrixType;
    typedef Epetra_Vector VectorType;
    typedef MCLS::VectorTraits<VectorType> VT;
    typedef MCLS::MatrixTraits<VectorType,MatrixType> MT;

    Teuchos::RCP<const Teuchos::Comm<int> > comm = 
	Teuchos::DefaultComm<int>::getComm();
    Teuchos::RCP<Epetra_Comm> epetra_comm = getEpetraComm( comm );
    int comm_size = comm->getSize();

    int local_num_rows = 10;
    int global_num_rows = local_num_rows*comm_size;
    Teuchos::RCP<Epetra_Map> map = Teuchos::rcp(
	new Epetra_Map( global_num_rows, 0, *epetra_comm ) );

    Teuchos::RCP<Epetra_CrsMatrix> A = 
	Teuchos::rcp( new Epetra_CrsMatrix( Copy, *map, 0 ) );

    Teuchos::Array<int> global_columns( 3 );
    double diag_val = 2.0;
    Teuchos::Array<double> values( 3, diag_val );
    for ( int i = 1; i < global_num_rows-1; ++i )
    {
	global_columns[0] = i-1;
	global_columns[1] = i;
	global_columns[2] = i+1;
	A->InsertGlobalValues( i, global_columns().size(), 
			       &values[0], &global_columns[0] );
    }
    Teuchos::Array<int> single_col(1,0);
    Teuchos::Array<double> diag_elem(1,diag_val);
    A->InsertGlobalValues( 0, 1, diag_elem.getRawPtr(), single_col.getRawPtr() );
    single_col[0] = global_num_rows-1;
    A->InsertGlobalValues( global_num_rows-1, 1, diag_elem.getRawPtr(), 
			   single_col.getRawPtr() );
    A->FillComplete();

    // Build the preconditioner.
    Teuchos::RCP<MCLS::Preconditioner<MatrixType> > preconditioner = 
	Teuchos::rcp( new MCLS::EpetraPointJacobiPreconditioner() );
    preconditioner->setOperator( A );
    preconditioner->buildPreconditioner();
    Teuchos::RCP<const MatrixType> M = preconditioner->getLeftPreconditioner();

    // Check the preconditioner.
    Teuchos::RCP<VectorType> X = MT::cloneVectorFromMatrixRows(*A);
    MT::getLocalDiagCopy( *M, *X );
    Teuchos::ArrayRCP<const double> X_view = VT::view( *X );
    Teuchos::ArrayRCP<const double>::const_iterator view_iterator;
    for ( view_iterator = X_view.begin();
	  view_iterator != X_view.end();
	  ++view_iterator )
    {
	TEST_EQUALITY( *view_iterator, 1.0/diag_val );
    }
}

Teuchos::RCP<Epetra_CrsMatrix> buildMatrix(int nx, Epetra_Comm & comm)
{
   Epetra_Map map(nx*comm.NumProc(),0,comm);
   Teuchos::RCP<Epetra_CrsMatrix> mat = Teuchos::rcp(new Epetra_CrsMatrix(Copy,map,3));

   int offsets[3] = {-1, 0, 1 };
   double values[3] = { -1, 2, -1};
   int maxGid = map.MaxAllGID();
   for(int lid=0;lid<nx;lid++) {
      int gid = mat->GRID(lid);
      int numEntries = 3, offset = 0;
      int indices[3] = { gid+offsets[0],
                         gid+offsets[1],
                         gid+offsets[2] };
      if(gid==0) { // left end point
         numEntries = 2;
         offset = 1;
      }            // right end point
      else if(gid==maxGid)
         numEntries = 2;

      // insert rows
      mat->InsertGlobalValues(gid,numEntries,values+offset,indices+offset);
   }

   mat->FillComplete();
   return mat;
}

/// building the stencil
void PoissonSolver::StencilGeometry(Teuchos::RCP<Epetra_Vector> & RHS,
                                    Teuchos::RCP<Epetra_CrsMatrix> & A)
{
    const Epetra_BlockMap & MyMap = RHS->Map();
    int NumMyElements = MyMap.NumMyElements();
    int* MyGlobalElements = MyMap.MyGlobalElements();
    double * rhsvalues = RHS->Values();

    std::vector<double> Values(5);
    std::vector<int> Indices(5);

    const auto & inside = _bend.getInsideMask();

    for (int lid = 0; lid < NumMyElements; ++ lid) {
        size_t NumEntries = 0;

        const size_t & gid = MyGlobalElements[lid];

        cutoffStencil(Indices,
                      Values,
                      rhsvalues[lid],
                      NumEntries,
                      inside,
                      gid);
        A->InsertGlobalValues(gid, NumEntries, &Values[0], &Indices[0]);
    }

    A->FillComplete();
    A->OptimizeStorage();
}

Teuchos::RCP<Epetra_CrsMatrix>
buildH( const Teuchos::RCP<Epetra_CrsMatrix> &A)
{
    Teuchos::RCP<Epetra_CrsMatrix> H = Teuchos::rcp( 
	new Epetra_CrsMatrix(Copy, A->RowMap(), A->GlobalMaxNumEntries() ) );
    int N = A->NumGlobalRows();
    std::vector<double> A_values( N );
    std::vector<int> A_indices( N );
    int A_size = 0;
    double local_H;
    bool found_diag = false;
    for ( int i = 0; i < N; ++i )
    {
	A->ExtractGlobalRowCopy( i,
				 N, 
				 A_size, 
				 &A_values[0], 
				 &A_indices[0] );

	for ( int j = 0; j < A_size; ++j )
	{
	    if ( i == A_indices[j] )
	    {
		local_H = 1.0 - A_values[j];
		H->InsertGlobalValues( i, 1, &local_H, &A_indices[j] );
		found_diag = true;
	    }
	    else
	    {
		local_H = -A_values[j];
		H->InsertGlobalValues( i, 1, &local_H, &A_indices[j] );
	    }
	}
	if ( !found_diag )
	{
	    local_H = 1.0;
	    H->InsertGlobalValues( i, 1, &local_H, &i );
	}
    }
    H->FillComplete();
    return H;
}

TEUCHOS_UNIT_TEST( OperatorTools, spectral_radius_test)
{
    int problem_size = 100;

    Epetra_SerialComm comm;
    Epetra_Map map( problem_size, 0, comm );

    // Build A.
    Teuchos::RCP<Epetra_CrsMatrix> A = 
	Teuchos::rcp( new Epetra_CrsMatrix( Copy, map, problem_size ) );
    double lower_diag = -1.0;
    double diag = 2.0;
    double upper_diag = -1.0;
    int global_row = 0;
    int lower_row = 0;
    int upper_row = 0;
    for ( int i = 0; i < problem_size; ++i )
    {
	global_row = A->GRID(i);
	lower_row = i-1;
	upper_row = i+1;
	if ( lower_row > -1 )
	{
	    A->InsertGlobalValues( global_row, 1, &lower_diag, &lower_row );
	}
	A->InsertGlobalValues( global_row, 1, &diag, &global_row );
	if ( upper_row < problem_size )
	{
	    A->InsertGlobalValues( global_row, 1, &upper_diag, &upper_row );
	}
    }
    A->FillComplete();

    double spec_rad_A = HMCSA::OperatorTools::spectralRadius( A );
    std::cout << spec_rad_A << std::endl;
}

Teuchos::RCP<Epetra_CrsMatrix> Epetra_Operator_to_Epetra_Matrix::constructInverseMatrix(const Epetra_Operator &op, const Epetra_Map &map)
{
  int numEntriesPerRow = 0;
  Teuchos::RCP<Epetra_FECrsMatrix> matrix = Teuchos::rcp(new Epetra_FECrsMatrix(::Copy, map, numEntriesPerRow));

  int numRows = map.NumGlobalElements();

  Epetra_Vector X(map);
  Epetra_Vector Y(map);

  double tol = 1e-15; // values below this will be considered 0

  for (int rowIndex=0; rowIndex<numRows; rowIndex++)
  {
    int lid = map.LID(rowIndex);
    if (lid != -1)
    {
      X[lid] = 1.0;
    }
    op.ApplyInverse(X, Y);
    if (lid != -1)
    {
      X[lid] = 0.0;
    }

    std::vector<double> values;
    std::vector<int> indices;
    for (int i=0; i<map.NumMyElements(); i++)
    {
      if (abs(Y[i]) > tol)
      {
        values.push_back(Y[i]);
        indices.push_back(map.GID(i));
      }
    }

    matrix->InsertGlobalValues(rowIndex, values.size(), &values[0], &indices[0]);
  }

  matrix->GlobalAssemble();
  return matrix;
}

// ****************************************************************************
void DG_Prob::DG_EI_Epetra_Outflow(const EDGE border,
																	 Teuchos::RCP<Epetra_FECrsMatrix> A,
																	 Teuchos::RCP<Epetra_FEVector> RHS)
{
	MyElem e0;
	e0=el[border.elemento[0]];
	const int ntot =  e0.show_ptr_stdel(sat)->nn_val() +
	e0.show_ptr_stdel(pres)->nn_val();
	
  double  mx  [ntot*ntot];
  double  B   [ntot];
  int     indx[ntot];

	const int ntot1 = DG_EI_Outflow(border,mx,B,indx);
	
	if(ntot1 > ntot){
		cout<<"DG_EI_Epetra_Outflow: ntot1 > ntot\n";
		exit(1);
	}
  
	A->InsertGlobalValues(ntot1,&indx[0],&mx[0],Epetra_FECrsMatrix::ROW_MAJOR);
  RHS->SumIntoGlobalValues(ntot1,&indx[0],&B[0]);

};

void
createEpetraProblem (const Teuchos::RCP<const Epetra_Comm>& epetraComm,
		     const std::string& filename,
		     Teuchos::RCP<Epetra_Map>& rowMap,
		     Teuchos::RCP<Epetra_CrsMatrix>& A,
		     Teuchos::RCP<Epetra_MultiVector>& B,
		     Teuchos::RCP<Epetra_MultiVector>& X,
		     int& numRHS)
{
  using Teuchos::inOutArg;
  using Teuchos::ptr;
  using Teuchos::RCP;
  using Teuchos::rcp;
  using Teuchos::set_extra_data;

  const int MyPID = epetraComm->MyPID();

  int n_nonzeros, N_update;
  int *bindx = NULL, *update = NULL, *col_inds = NULL;
  double *val = NULL, *row_vals = NULL;
  double *xguess = NULL, *b = NULL, *xexact = NULL;

  //
  // Set up the problem to be solved
  //
  int NumGlobalElements;  // total # of rows in matrix

  try {
    // Read in matrix from HB file
    Trilinos_Util_read_hb (const_cast<char *> (filename.c_str()), MyPID, 
			   &NumGlobalElements, &n_nonzeros, &val, &bindx, 
			   &xguess, &b, &xexact);

    // Distribute data among processors
    Trilinos_Util_distrib_msr_matrix (*epetraComm, &NumGlobalElements, 
				      &n_nonzeros, &N_update, &update, &val, 
				      &bindx, &xguess, &b, &xexact);
    //
    // Construct the matrix
    //
    int NumMyElements = N_update; // # local rows of matrix on processor
    //
    // Create an int array NumNz that is used to build the Petra Matrix.
    // NumNz[i] is the number of OFF-DIAGONAL terms for the i-th global
    // equation on this processor.
    //
    std::vector<int> NumNz (NumMyElements);
    for (int i = 0; i < NumMyElements; ++i) {
      NumNz[i] = bindx[i+1] - bindx[i] + 1;
    }
    rowMap = rcp (new Epetra_Map (NumGlobalElements, NumMyElements, 
				  update, 0, *epetraComm));
    //set_extra_data (epetraComm, "Map::Comm", inOutArg (rowMap));

    // Create an Epetra sparse matrix.
    if (NumMyElements == 0) {
      A = rcp (new Epetra_CrsMatrix (Copy, *rowMap, static_cast<int*>(NULL)));
    } else {
      A = rcp (new Epetra_CrsMatrix (Copy, *rowMap, &NumNz[0]));
    }
    //set_extra_data (rowMap, "Operator::Map", ptr (A));

    // Add rows to the sparse matrix one at a time.
    int NumEntries;
    for (int i = 0; i < NumMyElements; ++i) {
      row_vals = val + bindx[i];
      col_inds = bindx + bindx[i];
      NumEntries = bindx[i+1] - bindx[i];
      int info = A->InsertGlobalValues (update[i], NumEntries, row_vals, col_inds);
      TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				  "Failed to insert global value into A." );
      info = A->InsertGlobalValues (update[i], 1, val+i, update+i);
      TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				  "Failed to insert global value into A." );
    }
    // Finish initializing the sparse matrix.
    int info = A->FillComplete();
    TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				"FillComplete() failed on the sparse matrix A.");
    info = A->OptimizeStorage();
    TEUCHOS_TEST_FOR_EXCEPTION( info != 0, std::logic_error,
				"OptimizeStorage() failed on the sparse matrix A.");
    A->SetTracebackMode (1); // Shut down Epetra warning tracebacks
    //
    // Construct the right-hand side and solution multivectors.
    //
    if (false && b != NULL) {
      B = rcp (new Epetra_MultiVector (::Copy, *rowMap, b, NumMyElements, 1));
      numRHS = 1;
    } else {
      B = rcp (new Epetra_MultiVector (*rowMap, numRHS));
      B->Random ();
    }
    X = rcp (new Epetra_MultiVector (*rowMap, numRHS));
    X->PutScalar (0.0);

    //set_extra_data (rowMap, "X::Map", Teuchos::ptr (X));
    //set_extra_data (rowMap, "B::Map", Teuchos::ptr (B));
  } 
  catch (std::exception& e) {
//.........这里部分代码省略.........

//---------------------------------------------------------------------------//
TEUCHOS_UNIT_TEST( UniformForwardSource, nh_set )
{
    typedef Epetra_Vector VectorType;
    typedef MCLS::VectorTraits<VectorType> VT;
    typedef Epetra_RowMatrix MatrixType;
    typedef MCLS::MatrixTraits<VectorType,MatrixType> MT;
    typedef MCLS::ForwardHistory<int> HistoryType;
    typedef std::mt19937 rng_type;
    typedef MCLS::ForwardDomain<VectorType,MatrixType,rng_type> DomainType;

    Teuchos::RCP<const Teuchos::Comm<int> > comm = 
	Teuchos::DefaultComm<int>::getComm();
    Teuchos::RCP<Epetra_Comm> epetra_comm = getEpetraComm( comm );
    int comm_size = comm->getSize();

    int local_num_rows = 10;
    int global_num_rows = local_num_rows*comm_size;
    Teuchos::RCP<Epetra_Map> map = Teuchos::rcp(
	new Epetra_Map( global_num_rows, 0, *epetra_comm ) );

    // Build the linear system.
    Teuchos::RCP<Epetra_CrsMatrix> A = 	
	Teuchos::rcp( new Epetra_CrsMatrix( Copy, *map, 0 ) );
    Teuchos::Array<int> global_columns( 1 );
    Teuchos::Array<double> values( 1 );
    for ( int i = 1; i < global_num_rows; ++i )
    {
	global_columns[0] = i-1;
	values[0] = -0.5;
	A->InsertGlobalValues( i, global_columns().size(), 
			       &values[0], &global_columns[0] );
    }
    global_columns[0] = global_num_rows-1;
    values[0] = -0.5;
    A->InsertGlobalValues( global_num_rows-1, global_columns().size(),
			   &values[0], &global_columns[0] );
    A->FillComplete();

    Teuchos::RCP<MatrixType> A_T = MT::copyTranspose(*A);
    Teuchos::RCP<VectorType> x = MT::cloneVectorFromMatrixRows( *A );
    Teuchos::RCP<VectorType> b = MT::cloneVectorFromMatrixRows( *A );
    VT::putScalar( *b, -1.0 );

    // Build the forward domain.
    Teuchos::ParameterList plist;
    plist.set<int>( "Overlap Size", 0 );
    Teuchos::RCP<DomainType> domain = Teuchos::rcp( new DomainType( A_T, x, plist ) );

    // History setup.
    HistoryType::setByteSize();

    // Create the forward source with a set number of histories.
    int mult = 10;
    double cutoff = 1.0e-8;
    plist.set<double>("Sample Ratio",mult);
    plist.set<double>("Weight Cutoff", cutoff);
    MCLS::UniformForwardSource<DomainType> 
	source( b, domain, comm, comm->getSize(), comm->getRank(), plist );
    TEST_ASSERT( source.empty() );
    TEST_EQUALITY( source.numToTransport(), 0 );
    TEST_EQUALITY( source.numToTransportInSet(), mult*global_num_rows );
    TEST_EQUALITY( source.numRequested(), mult*global_num_rows );
    TEST_EQUALITY( source.numLeft(), 0 );
    TEST_EQUALITY( source.numEmitted(), 0 );

    // Build the source.
    source.buildSource();
    TEST_ASSERT( !source.empty() );
    TEST_EQUALITY( source.numToTransport(), mult*local_num_rows );
    TEST_EQUALITY( source.numToTransportInSet(), mult*global_num_rows );
    TEST_EQUALITY( source.numRequested(), mult*global_num_rows );
    TEST_EQUALITY( source.numLeft(), mult*local_num_rows );
    TEST_EQUALITY( source.numEmitted(), 0 );

    // Sample the source.
    Teuchos::RCP<MCLS::PRNG<rng_type> > rng = Teuchos::rcp(
	new MCLS::PRNG<rng_type>(comm->getRank()) );
    source.setRNG( rng );
    for ( int i = 0; i < mult*local_num_rows; ++i )
    {
	TEST_ASSERT( !source.empty() );
	TEST_EQUALITY( source.numLeft(), mult*local_num_rows-i );
	TEST_EQUALITY( source.numEmitted(), i );

	Teuchos::RCP<HistoryType> history = source.getHistory();

	TEST_EQUALITY( history->weight(), 1.0 );
	TEST_ASSERT( domain->isGlobalState( history->globalState() ) );
	TEST_ASSERT( history->alive() );
	TEST_ASSERT( VT::isGlobalRow( *x, history->globalState() ) );
    }
    TEST_ASSERT( source.empty() );
    TEST_EQUALITY( source.numLeft(), 0 );
    TEST_EQUALITY( source.numEmitted(), mult*local_num_rows );
}

int main(int argc, char** argv) {
#if defined(HAVE_MPI) && defined(HAVE_EPETRA)

  int numProcs = 1;
  int localProc = 0;

  //first, set up our MPI environment...
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &localProc);
  MPI_Comm_size(MPI_COMM_WORLD, &numProcs);

// This program can only run on 3 processors, to make things
// work out easy.
//
  if (numProcs != 3) {
    std::cout << "num-procs="<<numProcs<<". This program can only "
      << "run on 3 procs. Exiting."<<std::endl;
    MPI_Finalize();
    return(0);
  }

//Consider the following mesh of 4 2-D quad elements:
//
//  *-------*-------*
// 8|      7|      6|
//  |  E2   |  E3   |
//  *-------*-------*
// 3|      2|      5|
//  |  E0   |  E1   |
//  *-------*-------*
// 0       1       4
//
// Node-ids are to the lower-left of each node (*).
//
// Mimicing a finite-element application, we will say that
// each node has 1 scalar degree-of-freedom, and assemble
// a matrix which would have 9 global rows and columns.
//
// Each processor will have 3 rows. We'll set up a strange
// initial map, where nodes are distributed as follows:
//
// proc 0: nodes 0,3,8,
// proc 1: nodes 1,2,7
// proc 2: nodes 4,5,6.
//
// After we assemble our matrix, we'll create another matrix
// and populate it with graph edge weights such that the
// partitioner repartitions the problem so that nodes are
// laid out as follows:
//
// proc 0: nodes 0, 1, 4
// proc 1: nodes 3, 2, 5
// proc 2: nodes 8, 7, 6
//

  int nodesPerElem = 4;
  int global_n = 9;

  //First, set up the initial map:

  std::vector<int> mynodes(3);
  if (localProc == 0) {
    mynodes[0] = 0; mynodes[1] = 3; mynodes[2] = 8;
  }
  if (localProc == 1) {
    mynodes[0] = 1; mynodes[1] = 2; mynodes[2] = 7;
  }
  if (localProc == 2) {
    mynodes[0] = 4; mynodes[1] = 5; mynodes[2] = 6;
  }

  Epetra_MpiComm comm(MPI_COMM_WORLD);
  Epetra_Map origmap(global_n, 3, &mynodes[0], 0, comm);

  Teuchos::RCP<Epetra_FECrsMatrix> matrix =
    Teuchos::rcp(new Epetra_FECrsMatrix(Copy, origmap, 0));

  //We'll assemble elements E0 and E1 on proc 0,
  //               element E2 or proc 1,
  //               element E3 on proc 2.

  std::vector<int> indices(nodesPerElem);
  std::vector<double> coefs(nodesPerElem*nodesPerElem,2.0);

  if (localProc == 0) {
    //element E0:
    indices[0] = 0; indices[1] = 1; indices[2] = 2; indices[3] = 3;
    matrix->InsertGlobalValues(nodesPerElem, &indices[0], &coefs[0]);

    //element E1:
    indices[0] = 1; indices[1] = 4; indices[2] = 5; indices[3] = 2;
    matrix->InsertGlobalValues(nodesPerElem, &indices[0], &coefs[0]);
  }
  else if (localProc == 1) {
    //element E2:
    indices[0] = 3; indices[1] = 2; indices[2] = 7; indices[3] = 8;
    matrix->InsertGlobalValues(nodesPerElem, &indices[0], &coefs[0]);
  }
  else { //localProc==2
    //element E3:
//.........这里部分代码省略.........

int main(int argc, char *argv[])
{
  
  using Teuchos::rcp_implicit_cast;

  int i, ierr, gerr;
  gerr = 0;

#ifdef HAVE_MPI
  // Initialize MPI and setup an Epetra communicator
  MPI_Init(&argc,&argv);
  Teuchos::RCP<Epetra_MpiComm> Comm = Teuchos::rcp( new Epetra_MpiComm(MPI_COMM_WORLD) );
#else
  // If we aren't using MPI, then setup a serial communicator.
  Teuchos::RCP<Epetra_SerialComm> Comm = Teuchos::rcp( new Epetra_SerialComm() );
#endif


   // number of global elements
  int dim = 100;
  int blockSize = 3;

  // PID info
  int MyPID = Comm->MyPID();
  bool verbose = 0;

  if (argc>1) {
    if (argv[1][0]=='-' && argv[1][1]=='v') {
      verbose = true;
    }
  }

  // Construct a Map that puts approximately the same number of 
  // equations on each processor.
  Teuchos::RCP<Epetra_Map> Map = Teuchos::rcp( new Epetra_Map(dim, 0, *Comm) );
  
  // Get update list and number of local equations from newly created Map.
  int NumMyElements = Map->NumMyElements();
  std::vector<int> MyGlobalElements(NumMyElements);
  Map->MyGlobalElements(&MyGlobalElements[0]);

  // Create an integer std::vector NumNz that is used to build the Petra Matrix.
  // NumNz[i] is the Number of OFF-DIAGONAL term for the ith global equation 
  // on this processor
  std::vector<int> NumNz(NumMyElements);

  // We are building a tridiagonal matrix where each row has (-1 2 -1)
  // So we need 2 off-diagonal terms (except for the first and last equation)
  for (i=0; i<NumMyElements; i++) {
    if (MyGlobalElements[i]==0 || MyGlobalElements[i] == dim-1) {
      NumNz[i] = 2;
    }
    else {
      NumNz[i] = 3;
    }
  }

  // Create an Epetra_Matrix
  Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, *Map, &NumNz[0]) );
   
  // Add  rows one-at-a-time
  // Need some vectors to help
  // Off diagonal Values will always be -1
  std::vector<double> Values(2);
  Values[0] = -1.0; Values[1] = -1.0;
  std::vector<int> Indices(2);
  double two = 2.0;
  int NumEntries;
  for (i=0; i<NumMyElements; i++) {
    if (MyGlobalElements[i]==0) {
      Indices[0] = 1;
      NumEntries = 1;
    }
    else if (MyGlobalElements[i] == dim-1) {
      Indices[0] = dim-2;
      NumEntries = 1;
    }
    else {
      Indices[0] = MyGlobalElements[i]-1;
      Indices[1] = MyGlobalElements[i]+1;
      NumEntries = 2;
    }
    ierr = A->InsertGlobalValues(MyGlobalElements[i],NumEntries,&Values[0],&Indices[0]);
    assert(ierr==0);
    // Put in the diagonal entry
    ierr = A->InsertGlobalValues(MyGlobalElements[i],1,&two,&MyGlobalElements[i]);
    assert(ierr==0);
  }
   
  // Finish building the epetra matrix A
  ierr = A->FillComplete();
  assert(ierr==0);

  // Create an Belos::EpetraOp from this Epetra_CrsMatrix
  Teuchos::RCP<Belos::EpetraOp> op = Teuchos::rcp(new Belos::EpetraOp(A));

  // Issue several useful typedefs;
  typedef Belos::MultiVec<double> EMV;
  typedef Belos::Operator<double> EOP;

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

//.........这里部分代码省略.........
    }
    else {
      NumNz[i] = 5;
    }
  }

  // Create an Epetra_Matrix

  Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, Map, &NumNz[0]) );

  // Diffusion coefficient, can be set by user.
  // When rho*h/2 <= 1, the discrete convection-diffusion operator has real eigenvalues.
  // When rho*h/2 > 1, the operator has complex eigenvalues.
  double rho = 2*(nx+1);

  // Compute coefficients for discrete convection-diffution operator
  const double one = 1.0;
  std::vector<double> Values(4);
  std::vector<int> Indices(4);
  double h = one /(nx+1);
  double h2 = h*h;
  double c = 5.0e-01*rho/ h;
  Values[0] = -one/h2 - c; Values[1] = -one/h2 + c; Values[2] = -one/h2; Values[3]= -one/h2;
  double diag = 4.0 / h2;
  int NumEntries, info;

  for (int i=0; i<NumMyElements; i++)
  {
    if (MyGlobalElements[i]==0)
    {
      Indices[0] = 1;
      Indices[1] = nx;
      NumEntries = 2;
      info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
      assert( info==0 );
    }
    else if (MyGlobalElements[i] == nx*(nx-1))
    {
      Indices[0] = nx*(nx-1)+1;
      Indices[1] = nx*(nx-2);
      NumEntries = 2;
      info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[1], &Indices[0]);
      assert( info==0 );
    }
    else if (MyGlobalElements[i] == nx-1)
    {
      Indices[0] = nx-2;
      NumEntries = 1;
      info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
      assert( info==0 );
      Indices[0] = 2*nx-1;
      info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
      assert( info==0 );
    }
    else if (MyGlobalElements[i] == NumGlobalElements-1)
    {
      Indices[0] = NumGlobalElements-2;
      NumEntries = 1;
      info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
      assert( info==0 );
      Indices[0] = nx*(nx-1)-1;
      info = A->InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[2], &Indices[0]);
      assert( info==0 );
    }
    else if (MyGlobalElements[i] < nx)
    {

int extract_matrices
(
    Epetra_CrsMatrix *A,    // i/p: A matrix
    shylu_symbolic *ssym,   // symbolic structure
    shylu_data *data,       // numeric structure, TODO: Required ?
    shylu_config *config,   // i/p: library configuration
    bool insertValues       // true implies values will be inserted and fill
                            // complete will be called. false implies values
                            // will be replaced.
)
{
    Teuchos::RCP<Epetra_CrsMatrix> D = ssym->D;
    Teuchos::RCP<Epetra_CrsMatrix> C = ssym->C;
    Teuchos::RCP<Epetra_CrsMatrix> R = ssym->R;
    Teuchos::RCP<Epetra_CrsMatrix> G = ssym->G;
    Teuchos::RCP<Epetra_CrsGraph> Sg = ssym->Sg;
    int *DColElems = data->DColElems;
    int *gvals = data->gvals;
    double Sdiagfactor = config->Sdiagfactor;

    int *LeftIndex = new int[data->lmax];
    double *LeftValues = new double[data->lmax];
    int *RightIndex = new int[data->rmax];
    double *RightValues = new double[data->rmax];
    int err;
    int lcnt, rcnt ;
    int gcid;
    int gid;
    int *Ai;
    double *Ax;

    int nrows = A->RowMap().NumMyElements();
    int *rows = A->RowMap().MyGlobalElements();

    for (int i = 0; i < nrows ; i++)
    {
        int NumEntries;
        err = A->ExtractMyRowView(i, NumEntries, Ax, Ai);

        lcnt = 0; rcnt = 0;
        // Place the entry in the correct sub matrix, Works only for sym
        gid = rows[i];
        int lcid;
        for (int j = 0 ; j < NumEntries ; j++)
        { // O(nnz) ! Careful what you do inside
            // Row permutation does not matter here 
            gcid = A->GCID(Ai[j]);
            assert(gcid != -1);
            //Either in D or R 
            if ((gvals[gid] != 1 && gvals[gcid] == 1)
               || (gvals[gid] == 1 && A->LRID(gcid) != -1 && gvals[gcid] == 1))
            {
                assert(lcnt < data->lmax);
                if (insertValues)
                    LeftIndex[lcnt] = gcid;
                else
                {
                    //local column id
                    lcid = (gvals[gid] == 1 ? D->LCID(gcid) : R->LCID(gcid));
                    assert(lcid != -1);
                    LeftIndex[lcnt] = lcid;
                }
                LeftValues[lcnt++] = Ax[j];
            }
            else
            {
                assert(rcnt < data->rmax);
                if (insertValues)
                    RightIndex[rcnt] = gcid;
                else
                {
                    //local column id
                    lcid = (gvals[gid] == 1 ? C->LCID(gcid) : G->LCID(gcid));
                    assert(lcid != -1);
                    RightIndex[rcnt] = lcid;
                }
                RightValues[rcnt++] = Ax[j];
            }
        }

        if (gvals[gid] == 1)
        { // D or C row
            if (insertValues)
            {
                err = D->InsertGlobalValues(gid, lcnt, LeftValues, LeftIndex);
                assert(err == 0);
                err = C->InsertGlobalValues(gid, rcnt, RightValues, RightIndex);
                assert(err == 0);
            }
            else
            {
                err = D->ReplaceMyValues(D->LRID(gid), lcnt, LeftValues,
                                    LeftIndex);
                assert(err == 0);
                err = C->ReplaceMyValues(C->LRID(gid), rcnt, RightValues,
                                    RightIndex);
                assert(err == 0);
            }
        }
        else
//.........这里部分代码省略.........

//---------------------------------------------------------------------------//
TEUCHOS_UNIT_TEST( DomainTransporter, Boundary )
{
    typedef Epetra_Vector VectorType;
    typedef Epetra_RowMatrix MatrixType;
    typedef MCLS::MatrixTraits<VectorType,MatrixType> MT;
    typedef MCLS::AdjointHistory<int> HistoryType;
    typedef std::mt19937 rng_type;
    typedef MCLS::AdjointDomain<VectorType,MatrixType,rng_type> DomainType;

    Teuchos::RCP<const Teuchos::Comm<int> > comm = 
	Teuchos::DefaultComm<int>::getComm();
    Teuchos::RCP<Epetra_Comm> epetra_comm = getEpetraComm( comm );
    int comm_size = comm->getSize();
    int comm_rank = comm->getRank();

    // This test really needs a decomposed domain such that we can check
    // hitting the local domain boundary.
    if ( comm_size > 1 )
    {
	int local_num_rows = 10;
	int global_num_rows = local_num_rows*comm_size;
	Teuchos::RCP<Epetra_Map> map = Teuchos::rcp(
	    new Epetra_Map( global_num_rows, 0, *epetra_comm ) );

	// Build the linear operator and solution vector. This operator will
	// be assymetric so we quickly move the histories out of the domain
	// before they hit the low weight cutoff.
	Teuchos::RCP<Epetra_CrsMatrix> A = 	
	    Teuchos::rcp( new Epetra_CrsMatrix( Copy, *map, 0 ) );
	Teuchos::Array<int> global_columns( 3 );
	Teuchos::Array<double> values( 3 );

	global_columns[0] = 0;
	global_columns[1] = 1;
	global_columns[2] = 2;
	values[0] = 1.0;
	values[1] = 0.49;
	values[2] = 0.49;
	A->InsertGlobalValues( 0, global_columns.size(),
			       &values[0], &global_columns[0] );
	for ( int i = 1; i < global_num_rows-1; ++i )
	{
	    global_columns[0] = i-1;
	    global_columns[1] = i;
	    global_columns[2] = i+1;
	    values[0] = 0.49;
	    values[1] = 1.0;
	    values[2] = 0.49;
	    A->InsertGlobalValues( i, global_columns.size(),
				   &values[0], &global_columns[0] );
	}
	global_columns[0] = global_num_rows-3;
	global_columns[1] = global_num_rows-2;
	global_columns[2] = global_num_rows-1;
	values[0] = 0.49;
	values[1] = 0.49;
	values[2] = 1.0;
	A->InsertGlobalValues( global_num_rows-1, global_columns().size(), 
			       &values[0], &global_columns[0] );
	A->FillComplete();

	Teuchos::RCP<MatrixType> B = MT::copyTranspose(*A);
	Teuchos::RCP<VectorType> x = MT::cloneVectorFromMatrixRows( *A );

	// Build the adjoint domain.
	Teuchos::ParameterList plist;
	plist.set<int>( "Overlap Size", 2 );
	Teuchos::RCP<DomainType> domain = Teuchos::rcp( new DomainType( B, x, plist ) );
	Teuchos::RCP<MCLS::PRNG<rng_type> > rng = Teuchos::rcp(
	    new MCLS::PRNG<rng_type>( comm->getRank() ) );
	domain->setRNG( rng );

	// Build the domain transporter.
	MCLS::DomainTransporter<DomainType> transporter( domain, plist );
	domain->setCutoff( 1.0e-12 );

	// Transport histories through the domain until they hit a boundary.
	double weight = 3.0; 
	for ( int i = 0; i < global_num_rows-1; ++i )
	{
	    if ( comm_rank == comm_size - 1 )
	    {
		if ( i >= local_num_rows*comm_rank && i < local_num_rows*(comm_rank+1) )
		{
		    HistoryType history( i, i, weight );
		    history.live();
		    transporter.transport( history );

		    TEST_ASSERT( history.event() == 
				 MCLS::Event::BOUNDARY || MCLS::Event::CUTOFF );
		    TEST_ASSERT( !history.alive() );
		}
	    }
	    else
	    {
		if ( i >= local_num_rows*comm_rank && i < 2+local_num_rows*(comm_rank+1) )
		{
		    HistoryType history( i, i, weight );
		    history.live();
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
  int i, j, info;
  const double one = 1.0;
  const double zero = 0.0;
  Teuchos::LAPACK<int,double> lapack;

#ifdef EPETRA_MPI
  // Initialize MPI
  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  int MyPID = Comm.MyPID();

  //  Dimension of the matrix
  int m = 500;
  int n = 100;

  // Construct a Map that puts approximately the same number of
  // equations on each processor.

  Epetra_Map RowMap(m, 0, Comm);
  Epetra_Map ColMap(n, 0, Comm);

  // Get update list and number of local equations from newly created Map.

  int NumMyRowElements = RowMap.NumMyElements();
  
  std::vector<int> MyGlobalRowElements(NumMyRowElements);
  RowMap.MyGlobalElements(&MyGlobalRowElements[0]);

  /* We are building an m by n matrix with entries
    
              A(i,j) = k*(si)*(tj - 1) if i <= j
               = k*(tj)*(si - 1) if i  > j
  
     where si = i/(m+1) and tj = j/(n+1) and k = 1/(n+1).
  */

  // Create an Epetra_Matrix
  Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Copy, RowMap, n) );

  // Compute coefficients for discrete integral operator
  std::vector<double> Values(n);
  std::vector<int> Indices(n);
  double inv_mp1 = one/(m+1);
  double inv_np1 = one/(n+1);
  for (i=0; i<n; i++) { Indices[i] = i; }
  
  for (i=0; i<NumMyRowElements; i++) {
    //
    for (j=0; j<n; j++) {
      //
      if ( MyGlobalRowElements[i] <= j ) {
        Values[j] = inv_np1 * ( (MyGlobalRowElements[i]+one)*inv_mp1 ) * ( (j+one)*inv_np1 - one );  // k*(si)*(tj-1)
      }
      else {
        Values[j] = inv_np1 * ( (j+one)*inv_np1 ) * ( (MyGlobalRowElements[i]+one)*inv_mp1 - one );  // k*(tj)*(si-1)
      }
    }
    info = A->InsertGlobalValues(MyGlobalRowElements[i], n, &Values[0], &Indices[0]);
    assert( info==0 );
  }

  // Finish up
  info = A->FillComplete(ColMap, RowMap);
  assert( info==0 );
  info = A->OptimizeStorage();
  assert( info==0 );
  A->SetTracebackMode(1); // Shutdown Epetra Warning tracebacks

  //************************************
  // Start the block Arnoldi iteration
  //***********************************
  //
  //  Variables used for the Block Arnoldi Method
  //
  int nev = 4;
  int blockSize = 1;
  int numBlocks = 10;
  int maxRestarts = 20;
  int verbosity = Anasazi::Errors + Anasazi::Warnings + Anasazi::FinalSummary;
  double tol = lapack.LAMCH('E');
  std::string which = "LM";
  //
  // Create parameter list to pass into solver
  //
  Teuchos::ParameterList MyPL;
  MyPL.set( "Verbosity", verbosity );
  MyPL.set( "Which", which );
  MyPL.set( "Block Size", blockSize );
  MyPL.set( "Num Blocks", numBlocks );
  MyPL.set( "Maximum Restarts", maxRestarts );
  MyPL.set( "Convergence Tolerance", tol );

  typedef Anasazi::MultiVec<double> MV;
  typedef Anasazi::Operator<double> OP;

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