C++ Epetra_MultiVector::GlobalLength方法代码示例

int Epetra_PETScAIJMatrix::Multiply(bool TransA,
                               const Epetra_MultiVector& X,
                               Epetra_MultiVector& Y) const
{
  (void)TransA;
  int NumVectors = X.NumVectors();
  if (NumVectors!=Y.NumVectors()) EPETRA_CHK_ERR(-1);  // X and Y must have same number of vectors

  double ** xptrs;
  double ** yptrs;
  X.ExtractView(&xptrs);
  Y.ExtractView(&yptrs);
  if (RowMatrixImporter()!=0) {
    if (ImportVector_!=0) {
      if (ImportVector_->NumVectors()!=NumVectors) { delete ImportVector_; ImportVector_= 0;}
    }
    if (ImportVector_==0) ImportVector_ = new Epetra_MultiVector(RowMatrixColMap(),NumVectors);
    ImportVector_->Import(X, *RowMatrixImporter(), Insert);
    ImportVector_->ExtractView(&xptrs);
  }

  double *vals=0;
  int length;
  Vec petscX, petscY;
  int ierr;
  for (int i=0; i<NumVectors; i++) {
#   ifdef HAVE_MPI
    ierr=VecCreateMPIWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
    ierr=VecCreateMPIWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
#   else //FIXME  untested
    ierr=VecCreateSeqWithArray(Comm_->Comm(),X.MyLength(),X.GlobalLength(),xptrs[i],&petscX); CHKERRQ(ierr);
    ierr=VecCreateSeqWithArray(Comm_->Comm(),Y.MyLength(),Y.GlobalLength(),yptrs[i],&petscY); CHKERRQ(ierr);
#   endif

    ierr = MatMult(Amat_,petscX,petscY);CHKERRQ(ierr);

    ierr = VecGetArray(petscY,&vals);CHKERRQ(ierr);
    ierr = VecGetLocalSize(petscY,&length);CHKERRQ(ierr);
    for (int j=0; j<length; j++) yptrs[i][j] = vals[j];
    ierr = VecRestoreArray(petscY,&vals);CHKERRQ(ierr);
  }

  VecDestroy(petscX); VecDestroy(petscY);
  
  double flops = NumGlobalNonzeros();
  flops *= 2.0;
  flops *= (double) NumVectors;
  UpdateFlops(flops);
  return(0);
} //Multiply()

int writeMultiVector(FILE * handle, const Epetra_MultiVector & A, bool mmFormat) {

  int ierr = 0;
  int length = A.GlobalLength();
  int numVectors = A.NumVectors();
  const Epetra_Comm & comm = A.Map().Comm();
  if (comm.MyPID()!=0) {
    if (A.MyLength()!=0) ierr = -1;
  }
  else {
    if (length!=A.MyLength()) ierr = -1;
    for (int j=0; j<numVectors; j++) {
      for (int i=0; i<length; i++) {
	double val = A[j][i];
	if (mmFormat)
	  fprintf(handle, "%22.16e\n", val);
	else
	  fprintf(handle, "%22.16e ", val);
      }
      if (!mmFormat) fprintf(handle, "%s", "\n");
    }
  }
  int ierrGlobal;
  comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
  return(ierrGlobal);
}

// =============================================================================
void
VIO::EpetraMesh::Writer::
setValues( const Epetra_MultiVector          & x,
           const Teuchos::Array<std::string> & scalarsNames
         )
{
  unsigned int numVecs = x.NumVectors();
  unsigned int numVariables = x.GlobalLength();

  unsigned int numNodes = mesh_->getNodesMap()->NumGlobalElements();

  // make sure the sizes match the mesh
  if ( !mesh_.is_null() )
      TEUCHOS_ASSERT_EQUALITY( numVariables, 2*numNodes );

  // cast into a vtkUnstructuredGrid
  vtkSmartPointer<vtkUnstructuredGrid> vtkMesh =
      dynamic_cast<vtkUnstructuredGrid*> ( vtkDataSet_.GetPointer() );
  TEUCHOS_ASSERT_INEQUALITY( 0, !=, vtkMesh );

  // get scalarsNames, and insert default names if empty
  Teuchos::Array<std::string> scNames ( scalarsNames );
  if ( scNames.empty() )
  {
      scNames.resize ( numVecs );
      for ( int vec=0; vec<numVecs; vec++ )
          scNames[vec] = "x" + EpetraExt::toString ( vec );
  }

  // fill the scalar field
  vtkSmartPointer<vtkDoubleArray> scalars =
      vtkSmartPointer<vtkDoubleArray>::New();

  // real and imaginary part
  scalars->SetNumberOfComponents ( 2 );

  for ( int vec=0; vec<numVecs; vec++ )
  {
      scalars->SetName ( scNames[vec].c_str() );
      for ( int k=0; k<numNodes; k++ )
      {
//           const unsigned int dof_id = libmeshMesh_->node(k).dof_number(0,k,0);
          scalars->InsertNextValue ( x[vec][2*k] );
          scalars->InsertNextValue ( x[vec][2*k+1] );
      }
      vtkMesh->GetPointData()->AddArray ( scalars );
  }

  return;
}

// Convert a Epetra_MultiVector with assumed block structure dictated by the
// vector space into a Thyra::MultiVectorBase object.
// const Teuchos::RCP<const Thyra::MultiVectorBase<double> > blockEpetraToThyra(const Epetra_MultiVector & e,const Teuchos::RCP<const Thyra::VectorSpaceBase<double> > & vs)
void blockEpetraToThyra(const Epetra_MultiVector & epetraX,const Teuchos::Ptr<Thyra::MultiVectorBase<double> > & thyraX) 
{
   TEUCHOS_ASSERT(thyraX->range()->dim()==epetraX.GlobalLength());

   // extract local information from the Epetra_MultiVector
   int leadingDim=0,numVectors=0,localDim=0;
   double * epetraData=0;
   epetraX.ExtractView(&epetraData,&leadingDim);

   numVectors = epetraX.NumVectors();

   blockEpetraToThyra(numVectors,epetraData,leadingDim,thyraX.ptr(),localDim);   

   TEUCHOS_ASSERT(localDim==epetraX.MyLength());
}

// ============================================================================
void EpetraExt::XMLWriter::
Write(const std::string& Label, const Epetra_MultiVector& MultiVector)
{
  TEUCHOS_TEST_FOR_EXCEPTION(IsOpen_ == false, std::logic_error,
                     "No file has been opened");

  int Length = MultiVector.GlobalLength();
  int NumVectors = MultiVector.NumVectors();

  if (Comm_.MyPID() == 0)
  {
    std::ofstream of(FileName_.c_str(), std::ios::app);

    of << "<MultiVector Label=\"" << Label 
      << "\" Length=\"" << Length << '"'
      << " NumVectors=\"" << NumVectors << '"'
      << " Type=\"double\">" << std::endl;
  }


  for (int iproc = 0; iproc < Comm_.NumProc(); iproc++)
  {
    if (iproc == Comm_.MyPID())
    {
      std::ofstream of(FileName_.c_str(), std::ios::app);

      of.precision(15);
      for (int i = 0; i < MultiVector.MyLength(); ++i)
      {
        for (int j = 0; j < NumVectors; ++j)
          of << std::setiosflags(std::ios::scientific) << MultiVector[j][i] << " ";
        of << std::endl;
      }
      of.close();
    }
    Comm_.Barrier();
  }

  if (Comm_.MyPID() == 0)
  {
    std::ofstream of(FileName_.c_str(), std::ios::app);
    of << "</MultiVector>" << std::endl;
    of.close();
  }
}

int DoCopyMultiVector(double** matlabApr, const Epetra_MultiVector& A) {

  int ierr = 0;
  int length = A.GlobalLength();
  int numVectors = A.NumVectors();
  const Epetra_Comm & comm = A.Map().Comm();
  if (comm.MyPID()!=0) {
    if (A.MyLength()!=0) ierr = -1;
  }
  else {
    if (length!=A.MyLength()) ierr = -1;
    double* matlabAvalues = *matlabApr;
    double* Aptr = A.Values();
    memcpy((void *)matlabAvalues, (void *)Aptr, sizeof(*Aptr) * length * numVectors);
    *matlabApr += length;   
  }
  int ierrGlobal;
  comm.MinAll(&ierr, &ierrGlobal, 1); // If any processor has -1, all return -1
  return(ierrGlobal);
}

// Convert a Thyra::MultiVectorBase object to a Epetra_MultiVector object with
// the map defined by the Epetra_Map.
// const Teuchos::RCP<const Epetra_MultiVector> 
// blockThyraToEpetra(const Teuchos::RCP<const Thyra::MultiVectorBase<double> > & tX,const RCP<const Epetra_Map> & map)
void blockThyraToEpetra(const Teuchos::RCP<const Thyra::MultiVectorBase<double> > & thyraX,Epetra_MultiVector & epetraX)
{
   // build an Epetra_MultiVector object
   int numVectors = thyraX->domain()->dim();

   // make sure the number of vectors are the same
   TEUCHOS_ASSERT(numVectors==epetraX.NumVectors());
   TEUCHOS_ASSERT(thyraX->range()->dim()==epetraX.GlobalLength());

   // extract local information from the Epetra_MultiVector
   int leadingDim=0,localDim=0;
   double * epetraData=0;
   epetraX.ExtractView(&epetraData,&leadingDim);

   // perform recursive copy
   blockThyraToEpetra(numVectors,epetraData,leadingDim,thyraX,localDim);

   // sanity check
   TEUCHOS_ASSERT(localDim==epetraX.Map().NumMyElements());
}

int MultiVectorToMatrixMarketFile( const char *filename, const Epetra_MultiVector & A, 
				 const char * matrixName,
				 const char *matrixDescription, 
				 bool writeHeader) {
  int M = A.GlobalLength();
  int N = A.NumVectors();

  FILE * handle = 0;

  if (A.Map().Comm().MyPID()==0) { // Only PE 0 does this section

    handle = fopen(filename,"w");
    if (!handle) return(-1);
    MM_typecode matcode;
    mm_initialize_typecode(&matcode);
    mm_set_matrix(&matcode);
    mm_set_array(&matcode);
    mm_set_real(&matcode);

    if (writeHeader==true) { // Only write header if requested (true by default)
    
      if (mm_write_banner(handle, matcode)) return(-1);
      
      if (matrixName!=0) fprintf(handle, "%% \n%% %s\n", matrixName);
      if (matrixDescription!=0) fprintf(handle, "%% %s\n%% \n", matrixDescription);
      
      if (mm_write_mtx_array_size(handle, M, N)) return(-1);
    }
  }
    
  if (MultiVectorToMatrixMarketHandle(handle, A)) return(-1); // Everybody calls this routine

  if (A.Map().Comm().MyPID()==0) // Only PE 0 opened a file
    if (fclose(handle)) return(-1);
  return(0);
}

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

  int fail = 0, dim=0;

#ifdef HAVE_MPI
  MPI_Init(&argc, &argv);
  MPI_Comm_rank(MPI_COMM_WORLD, &localProc);
  MPI_Comm_size(MPI_COMM_WORLD, &numProcs);
  const Epetra_MpiComm Comm(MPI_COMM_WORLD);
#else
  const Epetra_SerialComm Comm;
#endif

  // =============================================================
  // get command line options
  // =============================================================

  Teuchos::CommandLineProcessor clp(false,true);

  std::string *inputFile = new std::string("simple.coords");
  bool verbose = false;

  clp.setOption( "f", inputFile, "Name of coordinate input file");

  clp.setOption( "v", "q", &verbose,
		"Display coordinates and weights before and after partitioning.");

  Teuchos::CommandLineProcessor::EParseCommandLineReturn parse_return =
    clp.parse(argc,argv);

  if( parse_return == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED){
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 0;
  }
  if( parse_return != Teuchos::CommandLineProcessor::PARSE_SUCCESSFUL ) {
#ifdef HAVE_MPI
    MPI_Finalize();
#endif
    return 1;
  }

  // =============================================================
  // Open file of coordinates and distribute them across processes
  // so they are unbalanced.
  // =============================================================

  Epetra_MultiVector *mv = ispatest::file2multivector(Comm, *inputFile);

  if (!mv || ((dim = mv->NumVectors()) < 1)){
    if (localProc == 0)
      std::cerr << "Invalid input file " << *inputFile << std::endl;
    exit(1);
  }

  if (localProc == 0){
    std::cerr << "Found input file " << *inputFile << ", " ;
    std::cerr << dim << " dimensional coordinates" << std::endl;
  }

  delete inputFile;

  int base = mv->Map().IndexBase();
  int globalSize = mv->GlobalLength();
  int myShare = 0;
  int n = numProcs - 1;

  if (n){
    if (localProc < n){
      int oneShare = globalSize / n;
      int leftOver = globalSize - (n * oneShare);
      myShare = oneShare + ((localProc < leftOver) ? 1 : 0);
    }
  }
  else{
    myShare = globalSize;
  }

  Epetra_BlockMap unbalancedMap(globalSize, myShare, 1, base, mv->Map().Comm());
  Epetra_Import importer(unbalancedMap, mv->Map());
  Epetra_MultiVector umv(unbalancedMap, dim);
  umv.Import(*mv, importer, Insert);

  delete mv;

  Teuchos::RCP<const Epetra_MultiVector> coords =
    Teuchos::rcp(new const Epetra_MultiVector(umv));

  // =============================================================
  // Create some different coordinate weight vectors
  // =============================================================

  Epetra_MultiVector *unitWgts = ispatest::makeWeights(coords->Map(), &ispatest::unitWeights);
  Epetra_MultiVector *veeWgts = ispatest::makeWeights(coords->Map(), &ispatest::veeWeights);
  Epetra_MultiVector *altWgts = ispatest::makeWeights(coords->Map(), &ispatest::alternateWeights);

  Teuchos::RCP<const Epetra_MultiVector> unit_weights_rcp = Teuchos::rcp(unitWgts);
  Teuchos::RCP<const Epetra_MultiVector> vee_weights_rcp = Teuchos::rcp(veeWgts);
  Teuchos::RCP<const Epetra_MultiVector> alt_weights_rcp = Teuchos::rcp(altWgts);
//.........这里部分代码省略.........

int ARPACKm3::reSolve(int numEigen, Epetra_MultiVector &Q, double *lambda, int startingEV) {

  // Computes eigenvalues and the corresponding eigenvectors
  // of the generalized eigenvalue problem
  // 
  //      K X = M X Lambda
  // 
  // using ARPACK (mode 3).
  //
  // The convergence test is provided by ARPACK.
  //
  // Note that if M is not specified, then  K X = X Lambda is solved.
  // (using the mode for generalized eigenvalue problem).
  // 
  // Input variables:
  // 
  // numEigen  (integer) = Number of eigenmodes requested
  // 
  // Q (Epetra_MultiVector) = Initial search space
  //                   The number of columns of Q defines the size of search space (=NCV).
  //                   The rows of X are distributed across processors.
  //                   As a rule of thumb in ARPACK User's guide, NCV >= 2*numEigen.
  //                   At exit, the first numEigen locations contain the eigenvectors requested.
  // 
  // lambda (array of doubles) = Converged eigenvalues
  //                   The length of this array is equal to the number of columns in Q.
  //                   At exit, the first numEigen locations contain the eigenvalues requested.
  // 
  // startingEV (integer) = Number of eigenmodes already stored in Q
  //                   A linear combination of these vectors is made to define the starting
  //                   vector, placed in resid.
  //
  // Return information on status of computation
  // 
  // info >=   0 >> Number of converged eigenpairs at the end of computation
  // 
  // // Failure due to input arguments
  // 
  // info = -  1 >> The stiffness matrix K has not been specified.
  // info = -  2 >> The maps for the matrix K and the matrix M differ.
  // info = -  3 >> The maps for the matrix K and the preconditioner P differ.
  // info = -  4 >> The maps for the vectors and the matrix K differ.
  // info = -  5 >> Q is too small for the number of eigenvalues requested.
  // info = -  6 >> Q is too small for the computation parameters.
  // 
  // info = -  8 >> numEigen must be smaller than the dimension of the matrix.
  //
  // info = - 30 >> MEMORY
  //
  // See ARPACK documentation for the meaning of INFO

  if (numEigen <= startingEV) {
    return numEigen;
  }

  int info = myVerify.inputArguments(numEigen, K, M, 0, Q, minimumSpaceDimension(numEigen));
  if (info < 0)
    return info;

  int myPid = MyComm.MyPID();

  int localSize = Q.MyLength();
  int NCV = Q.NumVectors();
  int knownEV = 0;

  if (NCV > Q.GlobalLength()) {
    if (numEigen >= Q.GlobalLength()) {
      cerr << endl;
      cerr << " !! The number of requested eigenvalues must be smaller than the dimension";
      cerr << " of the matrix !!\n";
      cerr << endl;
      return -8;
    }
    NCV = Q.GlobalLength();
  }

  int localVerbose = verbose*(myPid == 0);

  // Define data for ARPACK
  highMem = (highMem > currentSize()) ? highMem : currentSize();

  int ido = 0;

  int lwI = 22 + NCV;
  int *wI = new (nothrow) int[lwI];
  if (wI == 0) {
    return -30;
  }
  memRequested += sizeof(int)*lwI/(1024.0*1024.0);

  int *iparam = wI;
  int *ipntr = wI + 11;
  int *select = wI + 22;

  int lworkl = NCV*(NCV+8);
  int lwD = lworkl + 4*localSize;
  double *wD = new (nothrow) double[lwD];
  if (wD == 0) {
    delete[] wI;
    return -30;
//.........这里部分代码省略.........

int ModifiedARPACKm3::reSolve(int numEigen, Epetra_MultiVector &Q, double *lambda, 
                              int startingEV, const Epetra_MultiVector *orthoVec) {

  // Computes the smallest eigenvalues and the corresponding eigenvectors
  // of the generalized eigenvalue problem
  // 
  //      K X = M X Lambda
  // 
  // using ModifiedARPACK (mode 3).
  //
  // The convergence test is performed outisde of ARPACK
  //
  //                      || Kx - Mx lambda || < tol*lambda
  //
  // The norm ||.|| can be specified by the user through the array normWeight.
  // By default, the L2 Euclidean norm is used.
  //
  // Note that if M is not specified, then  K X = X Lambda is solved.
  // (using the mode for generalized eigenvalue problem).
  // 
  // Input variables:
  // 
  // numEigen  (integer) = Number of eigenmodes requested
  // 
  // Q (Epetra_MultiVector) = Initial search space
  //                   The number of columns of Q defines the size of search space (=NCV).
  //                   The rows of X are distributed across processors.
  //                   As a rule of thumb in ARPACK User's guide, NCV >= 2*numEigen.
  //                   At exit, the first numEigen locations contain the eigenvectors requested.
  // 
  // lambda (array of doubles) = Converged eigenvalues
  //                   The length of this array is equal to the number of columns in Q.
  //                   At exit, the first numEigen locations contain the eigenvalues requested.
  // 
  // startingEV (integer) = Number of eigenmodes already stored in Q
  //                   A linear combination of these vectors is made to define the starting
  //                   vector, placed in resid.
  //
  // orthoVec (Pointer to Epetra_MultiVector) = Space to be orthogonal to
  //                   The computation is performed in the orthogonal of the space spanned
  //                   by the columns vectors in orthoVec.
  //
  // Return information on status of computation
  // 
  // info >=   0 >> Number of converged eigenpairs at the end of computation
  // 
  // // Failure due to input arguments
  // 
  // info = -  1 >> The stiffness matrix K has not been specified.
  // info = -  2 >> The maps for the matrix K and the matrix M differ.
  // info = -  3 >> The maps for the matrix K and the preconditioner P differ.
  // info = -  4 >> The maps for the vectors and the matrix K differ.
  // info = -  5 >> Q is too small for the number of eigenvalues requested.
  // info = -  6 >> Q is too small for the computation parameters.
  // 
  // info = -  8 >> numEigen must be smaller than the dimension of the matrix.
  //
  // info = - 30 >> MEMORY
  //
  // See ARPACK documentation for the meaning of INFO

  if (numEigen <= startingEV) {
    return numEigen;
  }

  int info = myVerify.inputArguments(numEigen, K, M, 0, Q, minimumSpaceDimension(numEigen));
  if (info < 0)
    return info;

  int myPid = MyComm.MyPID();

  int localSize = Q.MyLength();
  int NCV = Q.NumVectors();
  int knownEV = 0;

  if (NCV > Q.GlobalLength()) {
    if (numEigen >= Q.GlobalLength()) {
      cerr << endl;
      cerr << " !! The number of requested eigenvalues must be smaller than the dimension";
      cerr << " of the matrix !!\n";
      cerr << endl;
      return -8;
    }
    NCV = Q.GlobalLength();
  }

  // Get the weight for approximating the M-inverse norm
  Epetra_Vector *vectWeight = 0;
  if (normWeight) {
    vectWeight = new Epetra_Vector(View, Q.Map(), normWeight);
  }

  int localVerbose = verbose*(myPid == 0);

  // Define data for ARPACK
  //
  // UH (10/17/03) Note that workl is also used 
  //               * to store the eigenvectors of the tridiagonal matrix
  //               * as a workspace for DSTEQR
  //               * as a workspace for recovering the global eigenvectors
//.........这里部分代码省略.........

  int  BuildMatrixTests (Epetra_MultiVector & C,
		       const char TransA, const char TransB, 
		       const double alpha, 
		       Epetra_MultiVector& A, 
		       Epetra_MultiVector& B,
		       const double beta,
		       Epetra_MultiVector& C_GEMM ) {

    // For given values of TransA, TransB, alpha and beta, a (possibly
    // zero) filled Epetra_MultiVector C, and allocated
    // Epetra_MultiVectors A, B and C_GEMM this routine will generate values for 
    // Epetra_MultiVectors A, B and C_GEMM such that, if A, B and (this) are 
    // used with GEMM in this class, the results should match the results 
    // generated by this routine.

    // Test for Strided multivectors (required for GEMM ops)

    if (!A.ConstantStride()   ||
	!B.ConstantStride()      ||
	!C_GEMM.ConstantStride() ||
	!C.ConstantStride()) return(-1); // Error 

    int i, j;
    double fi, fj;  // Used for casting loop variables to floats

    // Get a view of the MultiVectors

    double *Ap      = 0;
    double *Bp      = 0;
    double *Cp      = 0;
    double *C_GEMMp = 0;

    int A_nrows = A.MyLength();
    int A_ncols = A.NumVectors();
    int B_nrows = B.MyLength();
    int B_ncols = B.NumVectors();
    int C_nrows = C.MyLength();
    int C_ncols = C.NumVectors();
    int A_Stride         = 0;
    int B_Stride         = 0;
    int C_Stride         = 0;
    int C_GEMM_Stride    = 0;

    A.ExtractView(&Ap, &A_Stride);
    B.ExtractView(&Bp, &B_Stride);
    C.ExtractView(&Cp, &C_Stride);
    C_GEMM.ExtractView(&C_GEMMp, &C_GEMM_Stride);

      // Define some useful constants


    int opA_ncols = (TransA=='N') ? A.NumVectors() : A.MyLength();
    int opB_nrows = (TransB=='N') ? B.MyLength() : B.NumVectors();
  
    int C_global_inner_dim  = (TransA=='N') ? A.NumVectors() : A.GlobalLength();
  

    bool A_is_local = (!A.DistributedGlobal());
    bool B_is_local = (!B.DistributedGlobal());
    bool C_is_local = (!C.DistributedGlobal());

    int A_IndexBase = A.Map().IndexBase();
    int B_IndexBase = B.Map().IndexBase();
  
    // Build two new maps that we can use for defining global equation indices below
    Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
    Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());

    int* A_MyGlobalElements = new int[A_nrows];
    A_Map->MyGlobalElements(A_MyGlobalElements);
    int* B_MyGlobalElements = new int[B_nrows];
    B_Map->MyGlobalElements(B_MyGlobalElements);

  // Check for compatible dimensions

    if (C.MyLength()        != C_nrows     ||
	opA_ncols      != opB_nrows   ||
	C.NumVectors()    != C_ncols     ||
	C.MyLength()        != C_GEMM.MyLength()        ||
	C.NumVectors()    != C_GEMM.NumVectors()      ) {
      delete A_Map;
      delete B_Map;
      delete [] A_MyGlobalElements;
      delete [] B_MyGlobalElements;
      return(-2); // Return error
    }

    bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1 above
    bool Case2 = (!A_is_local && !B_is_local &&  C_is_local && TransA=='T' );// Case 2
    bool Case3 = (!A_is_local &&  B_is_local && !C_is_local && TransA=='N');// Case 3
  
    // Test for meaningful cases

    if (!(Case1 || Case2 || Case3)) {
      delete A_Map;
      delete B_Map;
      delete [] A_MyGlobalElements;
      delete [] B_MyGlobalElements;
      return(-3); // Meaningless case
    }
//.........这里部分代码省略.........

int  BuildMultiVectorTests (Epetra_MultiVector & C, const double alpha, 
				Epetra_MultiVector& A, 
				Epetra_MultiVector& sqrtA,
				Epetra_MultiVector& B,
				Epetra_MultiVector& C_alphaA,
				Epetra_MultiVector& C_alphaAplusB,
				Epetra_MultiVector& C_plusB,
				double* const dotvec_AB,
				double* const norm1_A,
				double* const norm2_sqrtA,
				double* const norminf_A,
				double* const normw_A,
				Epetra_MultiVector& Weights,
				double* const minval_A,
				double* const maxval_A,
				double* const meanval_A ) {

  // For given values alpha and a (possibly zero) filled 
  // Epetra_MultiVector (the this object), allocated double * arguments dotvec_AB, 
  // norm1_A, and norm2_A, and allocated Epetra_MultiVectors A, sqrtA,
  // B, C_alpha, C_alphaAplusB and C_plusB, this method will generate values for 
  // Epetra_MultiVectors A, B and all of the additional arguments on
  // the list above such that, if A, B and (this) are used with the methods in 
  // this class, the results should match the results generated by this routine.
  // Specifically, the results in dotvec_AB should match those from a call to 
  // A.dotProd (B,dotvec).  Similarly for other routines.
  
  int i,j;
  double fi, fj;  // Used for casting loop variables to floats
  // Define some useful constants
  
  int A_nrows = A.MyLength();
  int A_ncols = A.NumVectors();
  int sqrtA_nrows = sqrtA.MyLength();
  int sqrtA_ncols = sqrtA.NumVectors();
  int B_nrows = B.MyLength();
  int B_ncols = B.NumVectors();
  
  double **Ap = 0;
  double **sqrtAp = 0;
  double **Bp = 0;
  double **Cp = 0;
  double **C_alphaAp = 0;
  double **C_alphaAplusBp = 0;
  double **C_plusBp = 0;
  double **Weightsp = 0;

  A.ExtractView(&Ap);
  sqrtA.ExtractView(&sqrtAp);
  B.ExtractView(&Bp);
  C.ExtractView(&Cp);
  C_alphaA.ExtractView(&C_alphaAp);
  C_alphaAplusB.ExtractView(&C_alphaAplusBp);
  C_plusB.ExtractView(&C_plusBp);
  Weights.ExtractView(&Weightsp);
  
  bool A_is_local = (A.MyLength() == A.GlobalLength());
  bool B_is_local = (B.MyLength() == B.GlobalLength());
  bool C_is_local = (C.MyLength()    == C.GlobalLength());

  int A_IndexBase = A.Map().IndexBase();
  int B_IndexBase = B.Map().IndexBase();
  
    // Build two new maps that we can use for defining global equation indices below
    Epetra_Map * A_Map = new Epetra_Map(-1, A_nrows, A_IndexBase, A.Map().Comm());
    Epetra_Map * B_Map = new Epetra_Map(-1, B_nrows, B_IndexBase, B.Map().Comm());

    int* A_MyGlobalElements = new int[A_nrows];
    A_Map->MyGlobalElements(A_MyGlobalElements);
    int* B_MyGlobalElements = new int[B_nrows];
    B_Map->MyGlobalElements(B_MyGlobalElements);

  // Check for compatible dimensions
  
  if (C.MyLength()        != A_nrows     ||
      A_nrows        != B_nrows     ||
      C.NumVectors()    != A_ncols     ||
      A_ncols        != B_ncols     ||
      sqrtA_nrows    != A_nrows     ||
      sqrtA_ncols    != A_ncols     ||
      C.MyLength()        != C_alphaA.MyLength()     ||
      C.NumVectors()    != C_alphaA.NumVectors() ||
      C.MyLength()        != C_alphaAplusB.MyLength()     ||
      C.NumVectors()    != C_alphaAplusB.NumVectors() ||
      C.MyLength()        != C_plusB.MyLength()      ||
      C.NumVectors()    != C_plusB.NumVectors()     ) return(-2); // Return error
  
 
  bool Case1 = ( A_is_local &&  B_is_local &&  C_is_local);  // Case 1
  bool Case2 = (!A_is_local && !B_is_local && !C_is_local);// Case 2
  
  // Test for meaningful cases
  
  if (!(Case1 || Case2)) return(-3); // Meaningless case
  
  /* Fill A and B with values as follows:
     
     If A_is_local is false:
     A(i,j) = A_MyGlobalElements[i]*j,     i=1,...,numLocalEquations, j=1,...,NumVectors
     
//.........这里部分代码省略.........