void convert_snp_affymetrix_C(char **dirname_, char **filelist, unsigned *files_amount_, char **map_filename_, char **outfilename_, unsigned *skipaffym, char **alleleID_names, char *alleleID, unsigned *alleleID_amount)
{

char *outfilename = *outfilename_;
char *dirname = *dirname_;
char *map_filename = *map_filename_;


unsigned files_amount=*files_amount_;



std::map<std::string, char> coding;
for(unsigned i=0 ; i<*alleleID_amount ; i++)
	{
	coding[alleleID_names[i]] = alleleID[i];
	}


Rprintf("reading map...\n");
//std::cout<<"reading map...\n";
AffymetrixChipMap Map(map_filename, 2, 0, 2, 4, 5, 3, 9, 10, 6);
//std::cout<<"map is read...\n";
Rprintf("map is read...\n");

if(Map.get_exclude_amount() != 0) 
	{
	Rprintf("%i SNPs excluded from annotation because of absent enough information annotation file\n", Map.get_exclude_amount());			
	}



std::vector<ChipData *> ids_chip;
for(unsigned i=0 ; i<files_amount ; i++)
	{
	std::string file = (std::string(dirname) + "/" + std::string(filelist[i]));
	Rprintf("%i: opening file %s\n", i+1, file.c_str());
	ids_chip.push_back(new affymetrix_chip_data(file, 0, 1, *skipaffym));
	}





unsigned id_amount=ids_chip.size(); 


std::ofstream outfile(outfilename);
if(!outfile.is_open()){error("Can not open file \"\"\n",outfilename);}


Rprintf("Save to file %s\n", outfilename);


outfile << "#GenABEL raw data version 0.1\n";

//save IDs
Rprintf("saving Id names...\n");

std::string filelist_str;
for(unsigned id=0 ; id<files_amount ; id++)
	{
   filelist_str = filelist[id];
   for(unsigned i=0 ; i<filelist_str.length() ; i++)
      {
      if(filelist_str[i]==' ') filelist_str[i]='_';
      }

	outfile<<filelist_str<<" ";
	}



outfile<<"\n";

std::string snpname;


unsigned long snp_excludet_from_output_data=0;


//save snpnames
Rprintf("saving SNP names...\n");
unsigned snp_amount=ids_chip[0]->get_snp_amount();
for(unsigned snp=0 ; snp<snp_amount ; snp++)
	{
	snpname = ids_chip[0]->get_snp_name(snp);
	if(Map.is_snp_in_map(snpname)){outfile<<Map.recode_snp(snpname.c_str())<<" ";}
	else{snp_excludet_from_output_data++;}
	}
outfile<<"\n";

//save chromosome 
Rprintf("saving chromosome data...\n");
for(unsigned snp=0 ; snp<snp_amount ; snp++)
	{
	snpname = ids_chip[0]->get_snp_name(snp);
	if(Map.is_snp_in_map(snpname)){outfile<<Map.get_chromosome(snpname.c_str())<<" ";}
	}
outfile<<"\n";
//.........这里部分代码省略.........

//=========================================================================
int Epetra_IntVector::PackAndPrepare(const Epetra_SrcDistObject & Source,
                                     int NumExportIDs,
                                     int * ExportLIDs,
                                     int & LenExports,
                                     char * & Exports,
                                     int & SizeOfPacket,
                                     int * Sizes,
                                     bool & VarSizes,
                                     Epetra_Distributor & Distor)
{
  (void)Sizes;
  (void)VarSizes;
  (void)Distor;
  const Epetra_IntVector & A = dynamic_cast<const Epetra_IntVector &>(Source);

  int j, jj, k;

  int  * From;
  A.ExtractView(&From);
  int MaxElementSize = Map().MaxElementSize();
  bool ConstantElementSize = Map().ConstantElementSize();

  int * FromFirstPointInElementList = 0;
  int * FromElementSizeList = 0;

  if (!ConstantElementSize) {
    FromFirstPointInElementList = A.Map().FirstPointInElementList();
    FromElementSizeList = A.Map().ElementSizeList();
  }

  SizeOfPacket = MaxElementSize * (int)sizeof(int);

  if(NumExportIDs*SizeOfPacket>LenExports) {
    if (LenExports>0) delete [] Exports;
    LenExports = NumExportIDs*SizeOfPacket;
    Exports = new char[LenExports];
  }

  int * ptr;

  if (NumExportIDs>0) {
    ptr = (int *) Exports;

    // Point entry case
    if (MaxElementSize==1) for (j=0; j<NumExportIDs; j++) *ptr++ = From[ExportLIDs[j]];

    // constant element size case
    else if (ConstantElementSize) {

      for (j=0; j<NumExportIDs; j++) {
  jj = MaxElementSize*ExportLIDs[j];
    for (k=0; k<MaxElementSize; k++)
      *ptr++ = From[jj+k];
      }
    }

    // variable element size case
    else {

      int thisSizeOfPacket = MaxElementSize;
      for (j=0; j<NumExportIDs; j++) {
  ptr = (int *) Exports + j*thisSizeOfPacket;
  jj = FromFirstPointInElementList[ExportLIDs[j]];
  int ElementSize = FromElementSizeList[ExportLIDs[j]];
    for (k=0; k<ElementSize; k++)
      *ptr++ = From[jj+k];
      }
    }
  }

  return(0);
}

int main(int argc, char *argv[])
{
  int ierr = 0, forierr = 0;
  bool debug = false;

#ifdef EPETRA_MPI

  // Initialize MPI

  MPI_Init(&argc,&argv);
  int rank; // My process ID

  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  Epetra_MpiComm Comm( MPI_COMM_WORLD );

#else

  int rank = 0;
  Epetra_SerialComm Comm;

#endif

  bool verbose = false;

  // Check if we should print results to standard out
  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;

  int verbose_int = verbose ? 1 : 0;
  Comm.Broadcast(&verbose_int, 1, 0);
  verbose = verbose_int==1 ? true : false;


  //  char tmp;
  //  if (rank==0) cout << "Press any key to continue..."<< std::endl;
  //  if (rank==0) cin >> tmp;
  //  Comm.Barrier();

  Comm.SetTracebackMode(0); // This should shut down any error traceback reporting
  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc();

  if(verbose && MyPID==0)
    cout << Epetra_Version() << std::endl << std::endl;

  if (verbose) cout << "Processor "<<MyPID<<" of "<< NumProc
		    << " is alive."<<endl;

  bool verbose1 = verbose;

  // Redefine verbose to only print on PE 0
  if(verbose && rank!=0) 
		verbose = false;

  int NumMyEquations = 10000;
  int NumGlobalEquations = (NumMyEquations * NumProc) + EPETRA_MIN(NumProc,3);
  if(MyPID < 3) 
    NumMyEquations++;

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

  Epetra_Map Map(NumGlobalEquations, NumMyEquations, 0, Comm);
  
  // Get update list and number of local equations from newly created Map
  int* MyGlobalElements = new int[Map.NumMyElements()];
  Map.MyGlobalElements(MyGlobalElements);

  // Create an integer 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

  int* NumNz = new int[NumMyEquations];

  // 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 (int i = 0; i < NumMyEquations; i++)
    if((MyGlobalElements[i] == 0) || (MyGlobalElements[i] == NumGlobalEquations - 1))
      NumNz[i] = 1;
    else
      NumNz[i] = 2;

  // Create a Epetra_Matrix

  Epetra_CrsMatrix A(Copy, Map, NumNz);
  EPETRA_TEST_ERR(A.IndicesAreGlobal(),ierr);
  EPETRA_TEST_ERR(A.IndicesAreLocal(),ierr);
  
  // Add  rows one-at-a-time
  // Need some vectors to help
  // Off diagonal Values will always be -1


  double* Values = new double[2];
  Values[0] = -1.0; 
	Values[1] = -1.0;
  int* Indices = new int[2];
  double two = 2.0;
  int NumEntries;

  forierr = 0;
  for (int i = 0; i < NumMyEquations; i++) {
//.........这里部分代码省略.........

//=========================================================================
// Allows the source and target (\e this) objects to be compared for compatibility, return nonzero if not.
int EpetraExt_BlockDiagMatrix::CheckSizes(const Epetra_SrcDistObject& Source){
  return &Map() == &Source.Map();
}

void DataDrawer::AddFormat(const char *id, Factory factory)
{
	Mutex::Lock __(sDataDrawer);
	Map().Add(id, (void *)factory);
}

//=========================================================================
int EpetraExt_BlockDiagMatrix::ApplyInverse(const Epetra_MultiVector& X, Epetra_MultiVector& Y) const{
  int info;
  // Sanity Checks
  int NumVectors=X.NumVectors();  
  if(NumVectors!=Y.NumVectors())
    EPETRA_CHK_ERR(-1);
  if(!HasComputed_ && (ApplyMode_==AM_INVERT || ApplyMode_==AM_FACTOR))
    EPETRA_CHK_ERR(-2);
  
  //NTS: MultiVector's MyLength and [] Operators are  "points" level operators
  //not a "block/element" level operators.

  const int *vlist=DataMap_->FirstPointInElementList();
  const int *xlist=Map().FirstPointInElementList();
  const int *blocksize=Map().ElementSizeList();
  
  if(ApplyMode_==AM_MULTIPLY || ApplyMode_==AM_INVERT){
    // Multiply & Invert mode have the same apply
    int NumBlocks=NumMyBlocks();
    for(int i=0;i<NumBlocks;i++){
      int Nb=blocksize[i];
      int vidx0=vlist[i];
      int xidx0=xlist[i];
      for(int j=0;j<NumVectors;j++){	
	if(Nb==1) {
	  // Optimize for size = 1
	  Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
	}
	else if(Nb==2){
	  // Optimize for size = 2
	  Y[j][xidx0  ]=Values_[vidx0  ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
	  Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
	}
	else{
	  // "Large" Block - Use BLAS
	  //void 	GEMV (const char TRANS, const int M, const int N, const double ALPHA, const double *A, const int LDA, const double *X, const double BETA, double *Y, const int INCX=1, const int INCY=1) const 
	  GEMV('N',Nb,Nb,1.0,&Values_[vidx0],Nb,&X[j][xidx0],0.0,&Y[j][xidx0]);
	}
      }   
    }
  }
  else{
    // Factorization mode has a different apply
    int NumBlocks=NumMyBlocks();
    for(int i=0;i<NumBlocks;i++){
      int Nb=blocksize[i];
      int vidx0=vlist[i];
      int xidx0=xlist[i];      
      for(int j=0;j<NumVectors;j++){
	if(Nb==1) {
	  // Optimize for size = 1 - use the inverse
	  Y[j][xidx0]=Values_[vidx0]*X[j][xidx0];
	}
	else if(Nb==2){
	  // Optimize for size = 2 - use the inverse
	  Y[j][xidx0  ]=Values_[vidx0  ]*X[j][xidx0] + Values_[vidx0+2]*X[j][xidx0+1];
	  Y[j][xidx0+1]=Values_[vidx0+1]*X[j][xidx0] + Values_[vidx0+3]*X[j][xidx0+1];
	}
	else{
	  // "Large" Block - use LAPACK
	  //    void 	GETRS (const char TRANS, const int N, const int NRHS, const double *A, const int LDA, const int *IPIV, double *X, const int LDX, int *INFO) const 
	  for(int k=0;k<Nb;k++) Y[j][xidx0+k]=X[j][xidx0+k];
	  LAPACK.GETRS('N',Nb,1,&Values_[vidx0],Nb,&Pivots_[xidx0],&Y[j][xidx0],Nb,&info);
	  if(info) EPETRA_CHK_ERR(info);
	}
      }
    }    
  }  
  return 0;
}

int Drumm1(const Epetra_Map& map, bool verbose)
{
  (void)verbose;
  //Simple 2-element problem (element as in "finite-element") from
  //Clif Drumm. Two triangular elements, one per processor, as shown
  //here:
  //
  //   *----*
  //  3|\  2|
  //   | \  |
  //   | 0\1|
  //   |   \|
  //   *----*
  //  0    1
  //
  //Element 0 on processor 0, element 1 on processor 1.
  //Processor 0 will own nodes 0,1 and processor 1 will own nodes 2,3.
  //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
  //After GlobalAssemble(), the matrix should be as follows:
  //
  //         row 0: 2  1  0  1
  //proc 0   row 1: 1  4  1  2
  //----------------------------------
  //         row 2: 0  1  2  1
  //proc 1   row 3: 1  2  1  4
  //

  int numProcs = map.Comm().NumProc();
  int localProc = map.Comm().MyPID();

  if (numProcs != 2) return(0);

  //so first we'll set up a epetra_test::matrix_data object with
  //contents that match the above-described matrix. (but the
  //matrix_data object will have all 4 rows on each processor)

  int i;
  int rowlengths[4];
  rowlengths[0] = 3;
  rowlengths[1] = 4;
  rowlengths[2] = 3;
  rowlengths[3] = 4;

  epetra_test::matrix_data matdata(4, rowlengths);
  for(i=0; i<4; ++i) {
    for(int j=0; j<matdata.rowlengths()[i]; ++j) {
      matdata.colindices()[i][j] = j;
    }
  }

  matdata.colindices()[0][2] = 3;

  matdata.colindices()[2][0] = 1;
  matdata.colindices()[2][1] = 2;
  matdata.colindices()[2][2] = 3;

  double** coefs = matdata.coefs();
  coefs[0][0] = 2.0; coefs[0][1] = 1.0;                    coefs[0][2] = 1.0;
  coefs[1][0] = 1.0; coefs[1][1] = 4.0; coefs[1][2] = 1.0; coefs[1][3] = 2.0;
                     coefs[2][0] = 1.0; coefs[2][1] = 2.0; coefs[2][2] = 1.0;
  coefs[3][0] = 1.0; coefs[3][1] = 2.0; coefs[3][2] = 1.0; coefs[3][3] = 4.0;

  //now we'll load a Epetra_FECrsMatrix with data that matches the
  //above-described finite-element problem.

  int indexBase = 0, ierr = 0;
  int myNodes[4];
  double values[9];
  values[0] = 2.0;
  values[1] = 1.0;
  values[2] = 1.0;
  values[3] = 1.0;
  values[4] = 2.0;
  values[5] = 1.0;
  values[6] = 1.0;
  values[7] = 1.0;
  values[8] = 2.0;

  int numMyNodes = 2;

  if (localProc == 0) {
    myNodes[0] = 0;
    myNodes[1] = 1;
  }
  else {
    myNodes[0] = 2;
    myNodes[1] = 3;
  }

  Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());

  numMyNodes = 3;

  if (localProc == 0) {
    myNodes[0] = 0;
    myNodes[1] = 1;
    myNodes[2] = 3;
  }
  else {
    myNodes[0] = 1;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	}
	// Close the file.
	status = nc_close(ncid);
	if (status != NC_NOERR) handle_error(status);
      }

      // Convert from size_t to int.
      num_ss = total_rows;
      num_vars = cols0;

      std::cout<<"Number of snapshots: "<< num_ss << std::endl;
      std::cout<<"Length of snapshot : "<< num_vars << std::endl;
    }
    // Broadcast information about size of snapshot matrix.
    comm.Broadcast( &num_ss, 1, 0 );
    comm.Broadcast( &num_vars, 1, 0 );
    //
    // Sync all other processors on the scaling index vector if necessary
    // 
    if (createSSIdx) {
      for (i=0; i<(int)filenames.size(); i++) {
	if ( comm.MyPID() != 0 )
	  scaling_idx.push_back( idx_pair );
	comm.Broadcast( &scaling_idx[i].first, 1, 0 );
	comm.Broadcast( &scaling_idx[i].second, 1, 0 );
      }
      // Set the scaling index vector
      //params_->set("Snapshot Scaling Indices", scaling_idx);   
    }
    //
    // Create maps for new Epetra_MultiVector to hold the snapshots and 
    // temporary Epetra_Vector used by processor 0 to import the information.
    //
    Epetra_Map Map( num_vars, 0, comm );
    Teuchos::RCP<Epetra_MultiVector> newMV = Teuchos::rcp( new Epetra_MultiVector( Map, num_ss ) );
    Epetra_Vector *col_newMV = 0;
    Epetra_Map *Proc0Map = 0;
    int *index = 0;
    float *temp_vec_f = 0;
    double *temp_vec_d = 0;
    //
    if ( comm.MyPID() == 0 ) {
      Proc0Map = new Epetra_Map( num_vars, num_vars, 0, comm );
      temp_vec_f = new float [ num_vars ];
      temp_vec_d = new double [ num_vars ];
      index = new int[ num_vars ];
      for ( i=0; i<num_vars; i++ ) { index[i] = i; }
    } else {
      Proc0Map = new Epetra_Map( num_vars, 0, 0, comm );
    }
    //
    // Create an importer to get this information into the global Epetra_MultiVector
    //
    Epetra_Import importer( Map, *Proc0Map );
    //
    // Processor 0 reads each file and then creates a local Epetra_Vector, which will be
    // imported into the i-th column of the Epetra_MultiVector.
    //
    // Read starting with row "start2[0]" for "count2[0]" rows, as the columns vary from
    // "start2[1]" to "count2[1]", i.e. specifically for this case, read starting with row i 
    // for 1 row, as the columns vary from first column to the last column
    //
    start2[1]=0;
    count2[0]=1;
    count2[1]=num_vars;
    col_ptr = 0;

int Drumm2(const Epetra_Map& map, bool verbose)
{
  //Simple 2-element problem (element as in "finite-element") from
  //Clif Drumm. Two triangular elements, one per processor, as shown
  //here:
  //
  //   *----*
  //  3|\  2|
  //   | \  |
  //   | 0\1|
  //   |   \|
  //   *----*
  //  0    1
  //
  //Element 0 on processor 0, element 1 on processor 1.
  //Processor 0 will own nodes 0,1,3 and processor 1 will own node 2.
  //Each processor will pass a 3x3 element-matrix to Epetra_FECrsMatrix.
  //After GlobalAssemble(), the matrix should be as follows:
  //
  //         row 0: 2  1  0  1
  //proc 0   row 1: 1  4  1  2
  //         row 2: 0  1  2  1
  //----------------------------------
  //proc 1   row 3: 1  2  1  4
  //

  int numProcs = map.Comm().NumProc();
  int localProc = map.Comm().MyPID();

  if (numProcs != 2) return(0);

  int indexBase = 0, ierr = 0;

  double* values = new double[9];
  values[0] = 2.0;
  values[1] = 1.0;
  values[2] = 1.0;
  values[3] = 1.0;
  values[4] = 2.0;
  values[5] = 1.0;
  values[6] = 1.0;
  values[7] = 1.0;
  values[8] = 2.0;

  if (localProc == 0) {
    int numMyNodes = 3;
    int* myNodes = new int[numMyNodes];
    myNodes[0] = 0;
    myNodes[1] = 1;
    myNodes[2] = 3;

    Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());

    int rowLengths = 3;
    Epetra_FECrsMatrix A(Copy, Map, rowLengths);

    EPETRA_TEST_ERR( A.InsertGlobalValues(numMyNodes, myNodes,
			  numMyNodes, myNodes,
			  values, Epetra_FECrsMatrix::ROW_MAJOR),ierr);

    EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );
    EPETRA_TEST_ERR( A.GlobalAssemble(), ierr );

    if (verbose) {
    A.Print(cout);
    }

    //now let's make sure we can do a matvec with this matrix.
    Epetra_Vector x(Map), y(Map);
    x.PutScalar(1.0);
    EPETRA_TEST_ERR( A.Multiply(false, x, y), ierr);

    if (verbose&&localProc==0) {
      cout << "y = A*x, x=1.0's"<<endl;
    }

    if (verbose) {
    y.Print(cout);
    }

    delete [] myNodes;
    delete [] values;
  }
  else {
    int numMyNodes = 1;
    int* myNodes = new int[numMyNodes];
    myNodes[0] = 2;

    Epetra_Map Map(-1, numMyNodes, myNodes, indexBase, map.Comm());

    int rowLengths = 3;
    Epetra_FECrsMatrix A(Copy, Map, rowLengths);

    delete [] myNodes;
    numMyNodes = 3;
    myNodes = new int[numMyNodes];
    myNodes[0] = 1;
    myNodes[1] = 2;
    myNodes[2] = 3;

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

int Drumm3(const Epetra_Map& map, bool verbose)
{
  const Epetra_Comm & Comm = map.Comm();

  /* get number of processors and the name of this processor */

  int Numprocs = Comm.NumProc();
  int MyPID   = Comm.MyPID();

  if (Numprocs != 2) return(0);

  int NumGlobalRows = 4;
  int IndexBase = 0;
  Epetra_Map Map(NumGlobalRows, IndexBase, Comm);

  // Construct FECrsMatrix

  int NumEntriesPerRow = 3;

  Epetra_FECrsMatrix A(Copy, Map, NumEntriesPerRow);

  double ElementArea = 0.5;
  
  int NumCols = 3;
  int* Indices = new int[NumCols];

  if(MyPID==0)  // indices corresponding to element 0 on processor 0
  {
    Indices[0] = 0;
    Indices[1] = 1;
    Indices[2] = 3;
  }
  else if(MyPID==1)  // indices corresponding to element 1 on processor 1
  {
    Indices[0] = 1;
    Indices[1] = 2;
    Indices[2] = 3;
  }

  double* Values = new double[NumCols*NumCols];

// removal term
  Values[0] = 2*ElementArea/12.;
  Values[1] = 1*ElementArea/12.;
  Values[2] = 1*ElementArea/12.;
  Values[3] = 1*ElementArea/12.;
  Values[4] = 2*ElementArea/12.;
  Values[5] = 1*ElementArea/12.;
  Values[6] = 1*ElementArea/12.;
  Values[7] = 1*ElementArea/12.;
  Values[8] = 2*ElementArea/12.;

  A.InsertGlobalValues(NumCols, Indices,
                        Values,
                        Epetra_FECrsMatrix::ROW_MAJOR);

  A.GlobalAssemble();
  A.GlobalAssemble();

//  A.Print(cout);

// Create vectors for CG algorithm

  Epetra_FEVector* bptr = new Epetra_FEVector(A.RowMap(), 1);
  Epetra_FEVector* x0ptr = new Epetra_FEVector(A.RowMap(), 1);

  Epetra_FEVector& b = *bptr;
  Epetra_FEVector& x0 = *x0ptr;

  // source terms
  NumCols = 2;

  if(MyPID==0)  // indices corresponding to element 0 on processor 0
  {
    Indices[0] = 0;
    Indices[1] = 3;

    Values[0] = 1./2.;
    Values[1] = 1./2.;

   }
   else
   {
    Indices[0] = 1;
    Indices[1] = 2;

    Values[0] = 0;
    Values[1] = 0;
   }

  b.SumIntoGlobalValues(NumCols, Indices, Values);

  b.GlobalAssemble();

  if (verbose&&MyPID==0) cout << "b:" << endl;
  if (verbose) {
  b.Print(cout);
  }

  x0 = b;
//.........这里部分代码省略.........

 MultiVec<double>* EpetraMultiVec::Clone ( const int numvecs ) const
 {
   EpetraMultiVec * ptr_apv = new EpetraMultiVec(Map(), numvecs);
   return ptr_apv; // safe upcast.
 }

//=========================================================================
void Epetra_MapColoring::Print(std::ostream& os) const {
  int MyPID = Map().Comm().MyPID();
  int NumProc = Map().Comm().NumProc();
  
  if (MyPID==0) os 
    << std::endl 
    << " *****************************************" << std::endl
    << " Coloring information arranged map element" << std::endl 
    << " *****************************************" << std::endl
    << std::endl;
  for (int iproc=0; iproc < NumProc; iproc++) {
    if (MyPID==iproc) {
      int NumMyElements1 =Map(). NumMyElements();

      if (MyPID==0) {
  os.width(8);
  os <<  "     MyPID"; os << "    ";
  os.width(12);
  os <<  "GID  ";
  os.width(20);
  os <<  "Color  ";
  os << std::endl;
      }
      for (int i=0; i < NumMyElements1; i++) {
  os.width(10);
  os <<  MyPID; os << "    ";
  os.width(10);

    if(Map().GlobalIndicesInt()) {
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
        int * MyGlobalElements1 = Map().MyGlobalElements();
        os << MyGlobalElements1[i] << "    ";
#else
        throw ReportError("Epetra_MapColoring::Print: ERROR, GlobalIndicesInt but no API for it.",-1);
#endif
    }
    else if(Map().GlobalIndicesLongLong())
    {
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
    long long * MyGlobalElements1 = Map().MyGlobalElements64();
    os << MyGlobalElements1[i] << "    ";
#else
        throw ReportError("Epetra_MapColoring::Print: ERROR, GlobalIndicesLongLong but no API for it.",-1);
#endif
    }
    else
    throw ReportError("Epetra_MapColoring::Print: ERROR, Don't know map global index type.",-1);

  os.width(20);
  os <<  ElementColors_[i];
  os << std::endl;
      }
      os << std::flush; 
    }

    // Do a few global ops to give I/O a chance to complete
    Map().Comm().Barrier();
    Map().Comm().Barrier();
    Map().Comm().Barrier();
  }

  if (MyPID==0) os 
    << std::endl 
    << " **************************************" << std::endl
    << " Coloring information arranged by color" << std::endl 
    << " **************************************" << std::endl
    << std::endl;
  {for (int iproc=0; iproc < NumProc; iproc++) {
    if (MyPID==iproc) {
      if (NumColors()==0) os << " No colored elements on processor " << MyPID << std::endl;
      else {
        os << "Number of colors in map = " << NumColors() << std::endl
           << "Default color           = " << DefaultColor() << std::endl << std::endl;
        if (MyPID==0) {
          os.width(8);
          os <<  "     MyPID"; os << "    ";
          os.width(12);
          os <<  "LID  ";
          os.width(20);
          os <<  "Color  ";
          os << std::endl;
        }
        int * ColorValues = ListOfColors();
        for (int ii=0; ii<NumColors(); ii++) {
          int CV = ColorValues[ii];
    int ColorCount = NumElementsWithColor(CV);
    int * LIDList = ColorLIDList(CV);
    
    
    for (int i=0; i < ColorCount; i++) {
      os.width(10);
      os <<  MyPID; os << "    ";
      os.width(10);
      os << LIDList[i] << "    ";
      os.width(20);
      os << CV;
      os << std::endl;
    }
    os << std::flush; 
//.........这里部分代码省略.........

//=========================================================================
Epetra_MapColoring::~Epetra_MapColoring(){


  if (Allocated_ && Map().NumMyElements()>0) delete [] ElementColors_;
  if (ListsAreGenerated_) DeleteLists();
}

void DataDrawer::AddFormat(const char *id, Factory factory)
{
	INTERLOCKED_(sDataDrawer)
		Map().Add(id, (void *)factory);
}

int main(int argc, char *argv[]) {

#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  Epetra_SerialComm Comm;
#endif

  // set global dimension to 5, could be any number
  int NumGlobalElements = 5;
  
  // create a linear map
  Epetra_Map Map(NumGlobalElements,0,Comm);
  
  // local number of rows
  int NumMyElements = Map.NumMyElements();
  // get update list
  int * MyGlobalElements = Map.MyGlobalElements( );

  // dimension of each block
  Epetra_IntSerialDenseVector ElementSizeList(NumMyElements);

  // now construct a funky matrix. Diagonal block of block row i will have
  // dimension i+1 (don't run this code with too many nodes...). The
  // dimension of each block row is recordered in ElementSizeList.
  // Here ElementSizeList is declared as Epetra_IntSerialDenseVector,
  // but an int array is fine as well. 
  // max_blk keeps trace of the max block dimension
  
  int max_blk = 0;
  
  for( int i=0 ; i<NumMyElements ; ++i ) {
    ElementSizeList[i] = 1+MyGlobalElements[i];
    if( ElementSizeList[i] > max_blk ) max_blk =  ElementSizeList[i];
  }
  
  // create a block map based on the already declared point map
  // (used to determine NumMyElements and MyGlobalElements).
  // The same point map can be used for more block maps,
  // just change the input value of ElementSizeList
  Epetra_BlockMap BlockMap(NumGlobalElements,NumMyElements,
			   MyGlobalElements, 
			   ElementSizeList.Values(),0,Comm);

  // create a VBR matrix based on BlockMap
  Epetra_VbrMatrix A(Copy, BlockMap,2);

  int MaxBlockSize = max_blk * max_blk*100;
  
  int Indices[2];
  double* Values; Values = new double[MaxBlockSize];

  // cycle over all the local rows. 
  
  for( int i=0 ; i<NumMyElements ; ++i ) {
    
    // get GID of local row
    int GlobalNode = MyGlobalElements[i];
    // all lines but the last one will have to nonzero block-elements
    Indices[0] = GlobalNode;
    int NumEntries = 1;
    if( GlobalNode != NumGlobalElements-1 ) {
      Indices[1] = GlobalNode+1;
      NumEntries++;
    }

    // with VBR matrices, we have to insert one block at time.
    // This required two more instructions, one to start this
    // process (BeginInsertGlobalValues), and another one to
    // commit the end of submissions (EndSubmitEntries).
    
    A.BeginInsertGlobalValues(GlobalNode, NumEntries, Indices);
    // insert diagonal
    int BlockRows = ElementSizeList[i];
    for( int k=0 ; k<BlockRows * BlockRows ; ++k )
      Values[k] = 1.0*i;
    A.SubmitBlockEntry(Values,BlockRows,BlockRows,BlockRows);

    // insert off diagonal if any
    if( GlobalNode != NumGlobalElements-1 ) {
      int BlockCols = BlockRows+1;
      for( int k=0 ; k<BlockRows * BlockCols ; ++k )
	Values[k] = 1.0*i;
      A.SubmitBlockEntry(Values,BlockRows,BlockRows,BlockCols);
    }
    
    A.EndSubmitEntries();
  }

  A.FillComplete();

  cout << A;

  delete[] Values;

#ifdef HAVE_MPI
  MPI_Finalize();
#endif

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

int main(int argc, char *argv[]) {

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

  bool testFailed;
  bool boolret;
  int MyPID = Comm.MyPID();

  bool verbose = true;
  bool debug = false;
  bool insitu = false;
  std::string which("SM");

  Teuchos::CommandLineProcessor cmdp(false,true);
  cmdp.setOption("verbose","quiet",&verbose,"Print messages and results.");
  cmdp.setOption("debug","nodebug",&debug,"Print debugging information.");
  cmdp.setOption("insitu","exsitu",&insitu,"Perform in situ restarting.");
  cmdp.setOption("sort",&which,"Targetted eigenvalues (SM,LM,SR,LR,SI,or LI).");
  if (cmdp.parse(argc,argv) != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL) {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return -1;
  }

  typedef double ScalarType;
  typedef Teuchos::ScalarTraits<ScalarType>          SCT;
  typedef SCT::magnitudeType               MagnitudeType;
  typedef Epetra_MultiVector                          MV;
  typedef Epetra_Operator                             OP;
  typedef Anasazi::MultiVecTraits<ScalarType,MV>     MVT;
  typedef Anasazi::OperatorTraits<ScalarType,MV,OP>  OPT;

  //  Dimension of the matrix
  int nx = 10;        // Discretization points in any one direction.
  int NumGlobalElements = nx*nx;  // Size of matrix nx*nx

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

  Epetra_Map Map(NumGlobalElements, 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 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 matrix of block structure:

      | T -I          |
      |-I  T -I       |
      |   -I  T       |
      |        ...  -I|
      |           -I T|

   where each block is dimension nx by nx and the matrix is on the order of
   nx*nx.  The block T is a tridiagonal matrix. 
  */

  for (int i=0; i<NumMyElements; i++) {
    if (MyGlobalElements[i] == 0 || MyGlobalElements[i] == NumGlobalElements-1 || 
        MyGlobalElements[i] == nx-1 || MyGlobalElements[i] == nx*(nx-1) ) {
      NumNz[i] = 3;
    }
    else if (MyGlobalElements[i] < nx || MyGlobalElements[i] > nx*(nx-1) || 
             MyGlobalElements[i]%nx == 0 || (MyGlobalElements[i]+1)%nx == 0) {
      NumNz[i] = 4;
    }
    else {
      NumNz[i] = 5;
    }
  }

  // Create an Epetra_Matrix

  Teuchos::RCP<Epetra_CrsMatrix> A = Teuchos::rcp( new Epetra_CrsMatrix(Epetra_DataAccess::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);
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {

  int ierr = 0;

#ifdef EPETRA_MPI

  // Initialize MPI

  MPI_Init(&argc,&argv);
  Epetra_MpiComm Comm(MPI_COMM_WORLD);

#else

  Epetra_SerialComm Comm;

#endif

//  Comm.SetTracebackMode(0); // This should shut down any error tracing
  bool verbose = false;

  // Check if we should print results to standard out
  if (argc>1) if (argv[1][0]=='-' && argv[1][1]=='v') verbose = true;

#ifdef EPETRA_MPI
  int localverbose = verbose ? 1 : 0;
  int globalverbose=0;
  MPI_Allreduce(&localverbose, &globalverbose, 1, MPI_INT, MPI_SUM,
		MPI_COMM_WORLD);
  verbose = (globalverbose>0);
#endif

  int MyPID = Comm.MyPID();
  int NumProc = Comm.NumProc(); 

  if (verbose && MyPID==0)
    cout << Epetra_Version() << endl << endl;

  if (verbose) cout << Comm <<endl;

  // Redefine verbose to only print on PE 0
  //if (verbose && rank!=0) verbose = false;

  int NumMyElements = 4;
  int NumGlobalElements = NumMyElements*NumProc;
  int IndexBase = 0;
  
  Epetra_Map Map(NumGlobalElements, NumMyElements, IndexBase, Comm);

  EPETRA_TEST_ERR( Drumm1(Map, verbose),ierr);

  EPETRA_TEST_ERR( Drumm2(Map, verbose),ierr);

  EPETRA_TEST_ERR( Drumm3(Map, verbose),ierr);


  bool preconstruct_graph = false;

  EPETRA_TEST_ERR( four_quads(Comm, preconstruct_graph, verbose), ierr);

  preconstruct_graph = true;

  EPETRA_TEST_ERR( four_quads(Comm, preconstruct_graph, verbose), ierr);

  EPETRA_TEST_ERR( submatrix_formats(Comm, verbose), ierr);

  EPETRA_TEST_ERR( rectangular(Comm, verbose), ierr);

#ifdef EPETRA_MPI
  MPI_Finalize();
#endif

  return ierr;
}

//=========================================================================
int Epetra_IntVector::CopyAndPermute(const Epetra_SrcDistObject& Source,
                                     int NumSameIDs,
                                     int NumPermuteIDs,
                                     int * PermuteToLIDs,
                                     int *PermuteFromLIDs,
                                     const Epetra_OffsetIndex * Indexor,
                                     Epetra_CombineMode CombineMode)
{
  (void)Indexor;
  const Epetra_IntVector & A = dynamic_cast<const Epetra_IntVector &>(Source);

  if(    CombineMode != Add
      && CombineMode != Zero
      && CombineMode != Insert
      && CombineMode != Average
      && CombineMode != AbsMax )
    EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero

  int * From;
  A.ExtractView(&From);
  int *To = Values_;

  int * ToFirstPointInElementList = 0;
  int * FromFirstPointInElementList = 0;
  int * FromElementSizeList = 0;
  int MaxElementSize = Map().MaxElementSize();
  bool ConstantElementSize = Map().ConstantElementSize();

  if (!ConstantElementSize) {
    ToFirstPointInElementList =   Map().FirstPointInElementList();
    FromFirstPointInElementList = A.Map().FirstPointInElementList();
    FromElementSizeList = A.Map().ElementSizeList();
  }
  int j, jj, jjj, k;

  int NumSameEntries;

  bool Case1 = false;
  bool Case2 = false;
  // bool Case3 = false;

  if (MaxElementSize==1) {
    Case1 = true;
    NumSameEntries = NumSameIDs;
  }
  else if (ConstantElementSize) {
    Case2 = true;
    NumSameEntries = NumSameIDs * MaxElementSize;
  }
  else {
    // Case3 = true;
    NumSameEntries = FromFirstPointInElementList[NumSameIDs];
  }

  // Short circuit for the case where the source and target vector is the same.
  if (To==From) NumSameEntries = 0;

  // Do copy first
  if (NumSameIDs>0)
    if (To!=From) {
      if (CombineMode==Add)
  for (j=0; j<NumSameEntries; j++) To[j] += From[j]; // Add to existing value
      else if(CombineMode==Insert)
  for (j=0; j<NumSameEntries; j++) To[j] = From[j];
      else if(CombineMode==AbsMax)
        for (j=0; j<NumSameEntries; j++) {
    To[j] = EPETRA_MAX( To[j],From[j]);
  }
      // Note:  The following form of averaging is not a true average if more that one value is combined.
      //        This might be an issue in the future, but we leave this way for now.
      else if(CombineMode==Average)
  for (j=0; j<NumSameEntries; j++) {To[j] += From[j]; To[j] /= 2;}
    }
  // Do local permutation next
  if (NumPermuteIDs>0) {

    // Point entry case
    if (Case1) {

      if (CombineMode==Add)
  for (j=0; j<NumPermuteIDs; j++) To[PermuteToLIDs[j]] += From[PermuteFromLIDs[j]]; // Add to existing value
      else if(CombineMode==Insert)
  for (j=0; j<NumPermuteIDs; j++) To[PermuteToLIDs[j]] = From[PermuteFromLIDs[j]];
      else if(CombineMode==AbsMax)
        for (j=0; j<NumPermuteIDs; j++) {
    To[PermuteToLIDs[j]] = EPETRA_MAX( To[PermuteToLIDs[j]],From[PermuteFromLIDs[j]]);
  }
      // Note:  The following form of averaging is not a true average if more that one value is combined.
      //        This might be an issue in the future, but we leave this way for now.
      else if(CombineMode==Average)
  for (j=0; j<NumPermuteIDs; j++) {To[PermuteToLIDs[j]] += From[PermuteFromLIDs[j]]; To[PermuteToLIDs[j]] /= 2;}
    }
    // constant element size case
    else if (Case2) {

      if (CombineMode==Add)
      for (j=0; j<NumPermuteIDs; j++) {
  jj = MaxElementSize*PermuteToLIDs[j];
  jjj = MaxElementSize*PermuteFromLIDs[j];
//.........这里部分代码省略.........

//=========================================================================
// Print method
void EpetraExt_BlockDiagMatrix::Print(std::ostream & os) const{
  int MyPID = DataMap_->Comm().MyPID();
  int NumProc = DataMap_->Comm().NumProc();
  
  for (int iproc=0; iproc < NumProc; iproc++) {
    if (MyPID==iproc) {
      int NumMyElements1 =DataMap_->NumMyElements();
      int MaxElementSize1 = DataMap_->MaxElementSize();
	  int * MyGlobalElements1_int = 0;
	  long long * MyGlobalElements1_LL = 0;
#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES
  if(DataMap_->GlobalIndicesInt()) {
      MyGlobalElements1_int = DataMap_->MyGlobalElements();
  }
  else
#endif
#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES
  if(DataMap_->GlobalIndicesLongLong()) {
      MyGlobalElements1_LL = DataMap_->MyGlobalElements64();
  }
  else
#endif
    throw "EpetraExt_BlockDiagMatrix::Print: GlobalIndices type unknown";

      int * FirstPointInElementList1;
      if (MaxElementSize1!=1) FirstPointInElementList1 = DataMap_->FirstPointInElementList();

      if (MyPID==0) {
	os.width(8);
	os <<  "     MyPID"; os << "    ";
	os.width(12);
	if (MaxElementSize1==1)
	  os <<  "