自学教程：C++ EPETRA_CHK_ERR函数代码示例

int AztecOO_StatusTestResNorm::DefineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, 					       Epetra_Vector * Weights, 					       double ScaleValue ){  if (!firstcallDefineScaleForm_) EPETRA_CHK_ERR(-1); // We can only have this routine called once.  firstcallDefineScaleForm_ = false;  scaletype_ = TypeOfScaling;  scalenormtype_ = TypeOfNorm;  scaleweights_ = Weights;  scalevalue_ = ScaleValue;  // These conditions force the residual vector to be computed  if (scaletype_==NormOfInitRes && scalenormtype_!=TwoNorm) resvecrequired_ = true;    return(0);}

int Epetra_SerialDistributor::CreateFromRecvs( const int & NumRemoteIDs,				   const long long * RemoteGIDs,			           const int * RemotePIDs,				   bool Deterministic,			           int & NumExportIDs,				   long long *& ExportGIDs,				   int *& ExportPIDs ){  (void)NumRemoteIDs;  (void)RemoteGIDs;  (void)RemotePIDs;  (void)Deterministic;  (void)NumExportIDs;  (void)ExportGIDs;  (void)ExportPIDs;  EPETRA_CHK_ERR(-1); // This method should never be called  return(-1);}

//=========================================================================int Epetra_MapColoring::UnpackAndCombine(const Epetra_SrcDistObject & Source,           int NumImportIDs,                                         int * ImportLIDs,                                          int LenImports,           char * Imports,                                         int & SizeOfPacket,            Epetra_Distributor & Distor,            Epetra_CombineMode CombineMode,                                         const Epetra_OffsetIndex * Indexor ){  (void)Source;  (void)LenImports;  (void)Imports;  (void)SizeOfPacket;  (void)Distor;  (void)Indexor;  int j;    if(    CombineMode != Add      && CombineMode != Zero      && CombineMode != Insert      && CombineMode != AbsMax )    EPETRA_CHK_ERR(-1); //Unsupported CombinedMode, will default to Zero  if (NumImportIDs<=0) return(0);  int * To = ElementColors_;  int * ptr;  // Unpack it...  ptr = (int *) Imports;      if (CombineMode==Add)    for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] += ptr[j]; // Add to existing value  else if(CombineMode==Insert)    for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = ptr[j];  else if(CombineMode==AbsMax) {    for (j=0; j<NumImportIDs; j++) To[ImportLIDs[j]] = 0;    for (j=0; j<NumImportIDs; j++)  To[ImportLIDs[j]] = EPETRA_MAX( To[ImportLIDs[j]],std::abs(ptr[j]));  }    return(0);}

//=========================================================================int Epetra_LinearProblemRedistor::UpdateRedistProblemValues(Epetra_LinearProblem * ProblemWithNewValues) {		if (!RedistProblemCreated_) EPETRA_CHK_ERR(-1);  // This method can only be called after CreateRedistProblem()	Epetra_RowMatrix * OrigMatrix = ProblemWithNewValues->GetMatrix();	Epetra_MultiVector * OrigLHS = ProblemWithNewValues->GetLHS();	Epetra_MultiVector * OrigRHS = ProblemWithNewValues->GetRHS();			if (OrigMatrix==0) EPETRA_CHK_ERR(-2); // There is no matrix associated with this Problem	Epetra_CrsMatrix * RedistMatrix = dynamic_cast<Epetra_CrsMatrix *>(RedistProblem_->GetMatrix());	// Check if the tranpose should be create or not	if (ConstructTranspose_) {		EPETRA_CHK_ERR(Transposer_->UpdateTransposeValues(OrigMatrix));	}	else {		// If not, then just do the redistribution based on the the RedistMap				EPETRA_CHK_ERR(RedistMatrix->PutScalar(0.0));		// need to do this next step until we generalize the Import/Export ops for CrsMatrix		Epetra_CrsMatrix * OrigCrsMatrix = dynamic_cast<Epetra_CrsMatrix *>(OrigMatrix);		if (OrigCrsMatrix==0) EPETRA_CHK_ERR(-3); // Broken for a RowMatrix at this point		EPETRA_CHK_ERR(RedistMatrix->Export(*OrigCrsMatrix, *RedistExporter_, Add));	}	// Now redistribute the RHS and LHS if non-zero	if (OrigLHS!=0) {		EPETRA_CHK_ERR(RedistProblem_->GetLHS()->Export(*OrigLHS, *RedistExporter_, Add));	}	if (OrigRHS!=0) {		EPETRA_CHK_ERR(RedistProblem_->GetRHS()->Export(*OrigRHS, *RedistExporter_, Add));	}  return(0);}

//=============================================================================int Ifpack_CrsIct::Solve(bool Trans, const Epetra_MultiVector& X, 				Epetra_MultiVector& Y) const {//// This function finds Y such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS//  if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size  bool Upper = true;  bool UnitDiagonal = true;  Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;  Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;  U_->Solve(Upper, true, UnitDiagonal, *X1, *Y1);  Y1->Multiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)  U_->Solve(Upper, false, UnitDiagonal, *Y1, *Y1); // Solve Uy = y  return(0);}

//=========================================================================bool EpetraExt::RowMatrix_Transpose::fwd(){  Epetra_CrsMatrix * TempTransA1 = CreateTransposeLocal(*origObj_);  const Epetra_Export * TransposeExporter=0;  bool DeleteExporter = false;  if(TempTransA1->Exporter()) TransposeExporter = TempTransA1->Exporter();  else {    DeleteExporter=true;    TransposeExporter = new Epetra_Export(TransposeMatrix_->DomainMap(),TransposeMatrix_->RowMap());  }  TransposeMatrix_->PutScalar(0.0);  // Zero out all values of the matrix  EPETRA_CHK_ERR(TransposeMatrix_->Export(*TempTransA1, *TransposeExporter, Add));  if(DeleteExporter) delete TransposeExporter;  delete TempTransA1;  return 0;}

int EpetraVector<T>::GlobalAssemble(Epetra_CombineMode mode){  //In this method we need to gather all the non-local (overlapping) data  //that's been input on each processor, into the (probably) non-overlapping  //distribution defined by the map that 'this' vector was constructed with.  //We don't need to do anything if there's only one processor or if  //ignoreNonLocalEntries_ is true.  if (_vec->Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {    return(0);  }  //First build a map that describes the data in nonlocalIDs_/nonlocalCoefs_.  //We'll use the arbitrary distribution constructor of Map.  Epetra_BlockMap sourceMap(-1, numNonlocalIDs_,                            nonlocalIDs_, nonlocalElementSize_,                            _vec->Map().IndexBase(), _vec->Map().Comm());  //Now build a vector to hold our nonlocalCoefs_, and to act as the source-  //vector for our import operation.  Epetra_MultiVector nonlocalVector(sourceMap, 1);  int i,j;  for(i=0; i<numNonlocalIDs_; ++i) {    for(j=0; j<nonlocalElementSize_[i]; ++j) {      nonlocalVector.ReplaceGlobalValue(nonlocalIDs_[i], j, 0,                                        nonlocalCoefs_[i][j]);    }  }  Epetra_Export exporter(sourceMap, _vec->Map());  EPETRA_CHK_ERR( _vec->Export(nonlocalVector, exporter, mode) );  destroyNonlocalData();  return(0);}

int Epetra_FEVector::GlobalAssemble(Epetra_CombineMode mode,                                    bool reuse_map_and_exporter){  //In this method we need to gather all the non-local (overlapping) data  //that's been input on each processor, into the (probably) non-overlapping  //distribution defined by the map that 'this' vector was constructed with.  //We don't need to do anything if there's only one processor or if  //ignoreNonLocalEntries_ is true.  if (Map().Comm().NumProc() < 2 || ignoreNonLocalEntries_) {    return(0);  }  if (nonlocalMap_ == 0 || !reuse_map_and_exporter) {    createNonlocalMapAndExporter<int_type>();  }  Epetra_MultiVector& nonlocalVector = *nonlocalVector_;  nonlocalVector.PutScalar(0.0);  int elemSize = Map().MaxElementSize();  for(int vi=0; vi<NumVectors(); ++vi) {    for(size_t i=0; i<nonlocalIDs<int_type>().size(); ++i) {      for(int j=0; j<nonlocalElementSize_[i]; ++j) {        nonlocalVector.ReplaceGlobalValue(nonlocalIDs<int_type>()[i], j, vi,                                          nonlocalCoefs_[vi][i*elemSize+j]);      }    }  }  EPETRA_CHK_ERR( Export(nonlocalVector, *exporter_, mode) );  if (reuse_map_and_exporter) {    zeroNonlocalData<int_type>();  }  else {    destroyNonlocalData();  }  return(0);}

int Epetra_BlockMap::RemoteIDList(int NumIDs, const long long * GIDList,          int * PIDList, int * LIDList,          int * SizeList) const{  if(!BlockMapData_->GlobalIndicesLongLong_)    throw ReportError("Epetra_BlockMap::RemoteIDList ERROR, Can't call long long* version for non long long* map.",-1);  if (BlockMapData_->Directory_ == NULL) {    BlockMapData_->Directory_ = Comm().CreateDirectory(*this);  }  Epetra_Directory* directory = BlockMapData_->Directory_;  if (directory == NULL) {    return(-1);  }  EPETRA_CHK_ERR( directory->GetDirectoryEntries(*this, NumIDs, GIDList,             PIDList, LIDList, SizeList) );  return(0);}

int EpetraVector<T>::inputValues(int numIDs,                                 const int * GIDs,                                 const int * numValuesPerID,                                 const double * values,                                 bool accumulate){  if (accumulate) {    libmesh_assert(last_edit == 0 || last_edit == 2);    last_edit = 2;  } else {    libmesh_assert(last_edit == 0 || last_edit == 1);    last_edit = 1;  }  int offset=0;  for(int i=0; i<numIDs; ++i) {    int numValues = numValuesPerID[i];    if (_vec->Map().MyGID(GIDs[i])) {      if (accumulate) {        for(int j=0; j<numValues; ++j) {          _vec->SumIntoGlobalValue(GIDs[i], j, 0, values[offset+j]);        }      }      else {        for(int j=0; j<numValues; ++j) {          _vec->ReplaceGlobalValue(GIDs[i], j, 0, values[offset+j]);        }      }    }    else {      if (!ignoreNonLocalEntries_) {        EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,                                            &(values[offset]), accumulate) );      }    }    offset += numValues;  }  return(0);}

//==============================================================================int LinearProblem_CrsSingletonFilter::ConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,							     Epetra_Export * & RedistributeExporter,							     Epetra_Map * & RedistributeMap) {  int IndexBase = SourceMap->IndexBase();  if (IndexBase!=TargetMap->IndexBase()) EPETRA_CHK_ERR(-1);  const Epetra_Comm & Comm = TargetMap->Comm();  int TargetNumMyElements = TargetMap->NumMyElements();  int SourceNumMyElements = SourceMap->NumMyElements();  // ContiguousTargetMap has same number of elements per PE as TargetMap, but uses contigious indexing   Epetra_Map ContiguousTargetMap(-1, TargetNumMyElements, IndexBase,Comm);  // Same for ContiguousSourceMap  Epetra_Map ContiguousSourceMap(-1, SourceNumMyElements, IndexBase, Comm);  assert(ContiguousSourceMap.NumGlobalElements()==ContiguousTargetMap.NumGlobalElements());  // Now create a vector that contains the global indices of the Source Epetra_MultiVector  Epetra_IntVector SourceIndices(View, ContiguousSourceMap, SourceMap->MyGlobalElements());  // Create an exporter to send the SourceMap global IDs to the target distribution  Epetra_Export Exporter(ContiguousSourceMap, ContiguousTargetMap);    // Create a vector to catch the global IDs in the target distribution  Epetra_IntVector TargetIndices(ContiguousTargetMap);  TargetIndices.Export(SourceIndices, Exporter, Insert);  // Create a new map that describes how the Source MultiVector should be laid out so that it has  // the same number of elements on each processor as the TargetMap  RedistributeMap = new Epetra_Map(-1, TargetNumMyElements, TargetIndices.Values(), IndexBase, Comm);  // This exporter will finally redistribute the Source MultiVector to the same layout as the TargetMap  RedistributeExporter = new Epetra_Export(*SourceMap, *RedistributeMap);  return(0);}

//==============================================================================int Epetra_CrsSingletonFilter::ConstructRedistributeExporter(Epetra_Map * SourceMap, Epetra_Map * TargetMap,                   Epetra_Export * & RedistributeExporter,                   Epetra_Map * & RedistributeMap) {  int IndexBase = SourceMap->IndexBase(); // CJ TODO FIXME long long  if (IndexBase!=TargetMap->IndexBase()) EPETRA_CHK_ERR(-1);  const Epetra_Comm & Comm = TargetMap->Comm();  int TargetNumMyElements = TargetMap->NumMyElements();  int SourceNumMyElements = SourceMap->NumMyElements();  // ContiguousTargetMap has same number of elements per PE as TargetMap, but uses contigious indexing  Epetra_Map ContiguousTargetMap(-1, TargetNumMyElements, IndexBase,Comm);  // Same for ContiguousSourceMap  Epetra_Map ContiguousSourceMap(-1, SourceNumMyElements, IndexBase, Comm);  assert(ContiguousSourceMap.NumGlobalElements64()==ContiguousTargetMap.NumGlobalElements64());  // Now create a vector that contains the global indices of the Source Epetra_MultiVector#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES  Epetra_IntVector *SourceIndices = 0;  Epetra_IntVector *TargetIndices = 0;#endif#ifndef EPETRA_NO_64BIT_GLOBAL_INDICES  Epetra_LongLongVector *SourceIndices_LL = 0;  Epetra_LongLongVector *TargetIndices_LL = 0;#endif  if(SourceMap->GlobalIndicesInt())#ifndef EPETRA_NO_32BIT_GLOBAL_INDICES    SourceIndices = new Epetra_IntVector(View, ContiguousSourceMap, SourceMap->MyGlobalElements());#else    throw "Epetra_CrsSingletonFilter::ConstructRedistributeExporter: GlobalIndicesInt but no int API";#endif  else if(SourceMap->GlobalIndicesLongLong())

//=============================================================================int Ifpack_CrsIct::Multiply(bool Trans, const Epetra_MultiVector& X, 				Epetra_MultiVector& Y) const {//// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS//  if (X.NumVectors()!=Y.NumVectors()) EPETRA_CHK_ERR(-1); // Return error: X and Y not the same size  //bool Upper = true;  //bool Lower = false;  //bool UnitDiagonal = true;  Epetra_MultiVector * X1 = (Epetra_MultiVector *) &X;  Epetra_MultiVector * Y1 = (Epetra_MultiVector *) &Y;  U_->Multiply(false, *X1, *Y1);  Y1->Update(1.0, *X1, 1.0); // Y1 = Y1 + X1 (account for implicit unit diagonal)  Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0); // y = D*y (D_ has inverse of diagonal)  Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1  U_->Multiply(true, Y1temp, *Y1);  Y1->Update(1.0, Y1temp, 1.0); // (account for implicit unit diagonal)  return(0);}

int Epetra_FEVector::inputValues(int numIDs,                                 const int_type* GIDs,                                 const int* numValuesPerID,                                 const double* values,                                 bool suminto,                                 int vectorIndex){  if(!Map().template GlobalIndicesIsType<int_type>())  throw ReportError("Epetra_FEVector::inputValues mismatch between argument types (int/long long) and map type.", -1);  int offset=0;  for(int i=0; i<numIDs; ++i) {    int numValues = numValuesPerID[i];    if (Map().MyGID(GIDs[i])) {      if (suminto) {        for(int j=0; j<numValues; ++j) {          SumIntoGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);        }      }      else {        for(int j=0; j<numValues; ++j) {          ReplaceGlobalValue(GIDs[i], j, vectorIndex, values[offset+j]);        }      }    }    else {      if (!ignoreNonLocalEntries_) {        EPETRA_CHK_ERR( inputNonlocalValues(GIDs[i], numValues,              &(values[offset]), suminto,              vectorIndex) );      }    }    offset += numValues;  }  return(0);}

//=============================================================================int Ifpack_CrsRiluk::Multiply(bool Trans, const Epetra_MultiVector& X, 			      Epetra_MultiVector& Y) const {//// This function finds X such that LDU Y = X or U(trans) D L(trans) Y = X for multiple RHS//      // First generate X and Y as needed for this function  Teuchos::RefCountPtr<Epetra_MultiVector> X1;  Teuchos::RefCountPtr<Epetra_MultiVector> Y1;  EPETRA_CHK_ERR(GenerateXY(Trans, X, Y, &X1, &Y1));#ifdef IFPACK_FLOPCOUNTERS  Epetra_Flops * counter = this->GetFlopCounter();  if (counter!=0) {    L_->SetFlopCounter(*counter);    Y1->SetFlopCounter(*counter);    U_->SetFlopCounter(*counter);  }#endif  if (!Trans) {    EPETRA_CHK_ERR(U_->Multiply(Trans, *X1, *Y1)); //     EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)    EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)    Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1    EPETRA_CHK_ERR(L_->Multiply(Trans, Y1temp, *Y1));    EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)    if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));} // Export computed Y values if needed  }  else {    EPETRA_CHK_ERR(L_->Multiply(Trans, *X1, *Y1));    EPETRA_CHK_ERR(Y1->Update(1.0, *X1, 1.0)); // Y1 = Y1 + X1 (account for implicit unit diagonal)    EPETRA_CHK_ERR(Y1->ReciprocalMultiply(1.0, *D_, *Y1, 0.0)); // y = D*y (D_ has inverse of diagonal)    Epetra_MultiVector Y1temp(*Y1); // Need a temp copy of Y1    EPETRA_CHK_ERR(U_->Multiply(Trans, Y1temp, *Y1));    EPETRA_CHK_ERR(Y1->Update(1.0, Y1temp, 1.0)); // (account for implicit unit diagonal)    if (IsOverlapped_) {EPETRA_CHK_ERR(Y.Export(*Y1,*L_->Exporter(), OverlapMode_));}  }   return(0);}

//==========================================================================int Ifpack_CrsIct::InitValues(const Epetra_CrsMatrix & A) {  int ierr = 0;  int i, j;  int NumIn, NumL, NumU;  bool DiagFound;  int NumNonzeroDiags = 0;  Teuchos::RefCountPtr<Epetra_CrsMatrix> OverlapA = Teuchos::rcp( (Epetra_CrsMatrix *) &A_ , false );  if (LevelOverlap_>0) {    EPETRA_CHK_ERR(-1); // Not implemented yet    //OverlapA = new Epetra_CrsMatrix(Copy, *Graph_.OverlapGraph());    //EPETRA_CHK_ERR(OverlapA->Import(A, *Graph_.OverlapImporter(), Insert));    //EPETRA_CHK_ERR(OverlapA->FillComplete());  }  // Get Maximun Row length  int MaxNumEntries = OverlapA->MaxNumEntries();  vector<int> InI(MaxNumEntries); // Allocate temp space  vector<int> UI(MaxNumEntries);  vector<double> InV(MaxNumEntries);  vector<double> UV(MaxNumEntries);  double *DV;  ierr = D_->ExtractView(&DV); // Get view of diagonal      // First we copy the user's matrix into diagonal vector and U, regardless of fill level  int NumRows = OverlapA->NumMyRows();  for (i=0; i< NumRows; i++) {    OverlapA->ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0]); // Get Values and Indices        // Split into L and U (we don't assume that indices are ordered).        NumL = 0;     NumU = 0;     DiagFound = false;        for (j=0; j< NumIn; j++) {      int k = InI[j];      if (k==i) {	DiagFound = true;	DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_      }      else if (k < 0) return(-1); // Out of range      else if (i<k && k<NumRows) {	UI[NumU] = k;	UV[NumU] = InV[j];	NumU++;      }    }        // Check in things for this row of L and U    if (DiagFound) NumNonzeroDiags++;    if (NumU) U_->InsertMyValues(i, NumU, &UV[0], &UI[0]);      }  U_->FillComplete(A_.OperatorDomainMap(), A_.OperatorRangeMap());  SetValuesInitialized(true);  SetFactored(false);  int ierr1 = 0;  if (NumNonzeroDiags<U_->NumMyRows()) ierr1 = 1;  A_.Comm().MaxAll(&ierr1, &ierr, 1);  EPETRA_CHK_ERR(ierr);  return(0);}

//==========================================================================int Ifpack_CrsIct::Factor() {  // if (!Allocated()) return(-1); // This test is not needed at this time.  All constructors allocate.  if (!ValuesInitialized_) EPETRA_CHK_ERR(-2); // Must have values initialized.  if (Factored()) EPETRA_CHK_ERR(-3); // Can't have already computed factors.  SetValuesInitialized(false);  int i;  int m, n, nz, Nrhs, ldrhs, ldlhs;  int * ptr=0, * ind;  double * val, * rhs, * lhs;  int ierr = Epetra_Util_ExtractHbData(U_.get(), 0, 0, m, n, nz, ptr, ind,			    val, Nrhs, rhs, ldrhs, lhs, ldlhs);  if (ierr<0) EPETRA_CHK_ERR(ierr);  Matrix * Aict;  if (Aict_==0) {    Aict = new Matrix;    Aict_ = (void *) Aict;  }  else Aict = (Matrix *) Aict_;  Matrix * Lict;  if (Lict_==0) {    Lict = new Matrix;    Lict_ = (void *) Lict;  }  else Lict = (Matrix *) Lict_;  Aict->val = val;  Aict->col = ind;  Aict->ptr = ptr;  double *DV;  EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal      crout_ict(m, Aict, DV, Droptol_, Lfil_, Lict, &Ldiag_);  // Get rid of unnecessary data  delete [] ptr;  // Create Epetra View of L from crout_ict  if (LevelOverlap_==0) {    U_ = Teuchos::rcp( new Epetra_CrsMatrix(View, A_.RowMatrixRowMap(), A_.RowMatrixRowMap(),0) );    D_ = Teuchos::rcp( new Epetra_Vector(View, A_.RowMatrixRowMap(), Ldiag_) );  }  else {    EPETRA_CHK_ERR(-1); // LevelOverlap > 0 not implemented yet    //    U_ = new Epetra_CrsMatrix(Copy, OverlapRowMap());    //    D_ = new Epetra_Vector(OverlapRowMap());  }  ptr = Lict->ptr;  ind = Lict->col;  val = Lict->val;      for (i=0; i< m; i++) {    int NumEntries = ptr[i+1]-ptr[i];    int * Indices = ind+ptr[i];    double * Values = val+ptr[i];    U_->InsertMyValues(i, NumEntries, Values, Indices);  }  U_->FillComplete(A_.OperatorDomainMap(), A_.OperatorRangeMap());    D_->Reciprocal(*D_); // Put reciprocal of diagonal in this vector  // Add up flops   double current_flops = 2 * nz; // Just an estimate  double total_flops = 0;      A_.Comm().SumAll(&current_flops, &total_flops, 1); // Get total madds across all PEs  // Now count the rest  total_flops += (double) U_->NumGlobalNonzeros(); // Accounts for multiplier above  total_flops += (double) D_->GlobalLength(); // Accounts for reciprocal of diagonal  UpdateFlops(total_flops); // Update flop count  SetFactored(true);  return(0);}

//=========================================================================int Epetra_Vector::SumIntoMyValues(int NumEntries, int BlockOffset, const double * values, const int * Indices) {  // Use the more general method below  EPETRA_CHK_ERR(ChangeValues(NumEntries, BlockOffset, values, Indices, false, true));  return(0);}

////  Amesos_TestMultiSolver.cpp reads in a matrix in Harwell-Boeing format, //  calls one of the sparse direct solvers, using blocked right hand sides//  and computes the error and residual.  ////  TestSolver ignores the Harwell-Boeing right hand sides, creating//  random right hand sides instead.  ////  Amesos_TestMultiSolver can test either A x = b or A^T x = b.//  This can be a bit confusing because sparse direct solvers //  use compressed column storage - the transpose of Trilinos'//  sparse row storage.////  Matrices://    readA - Serial.  As read from the file.//    transposeA - Serial.  The transpose of readA.//    serialA - if (transpose) then transposeA else readA //    distributedA - readA distributed to all processes//    passA - if ( distributed ) then distributedA else serialA////int Amesos_TestMultiSolver( Epetra_Comm &Comm, char *matrix_file, int numsolves, 		      SparseSolverType SparseSolver, bool transpose,		      int special, AMESOS_MatrixType matrix_type ) {  int iam = Comm.MyPID() ;    //  int hatever;  //  if ( iam == 0 )  std::cin >> hatever ;   Comm.Barrier();  Epetra_Map * readMap;  Epetra_CrsMatrix * readA;   Epetra_Vector * readx;   Epetra_Vector * readb;  Epetra_Vector * readxexact;     std::string FileName = matrix_file ;  int FN_Size = FileName.size() ;   std::string LastFiveBytes = FileName.substr( EPETRA_MAX(0,FN_Size-5), FN_Size );  std::string LastFourBytes = FileName.substr( EPETRA_MAX(0,FN_Size-4), FN_Size );  bool NonContiguousMap = false;   if ( LastFiveBytes == ".triU" ) {     NonContiguousMap = true;     // Call routine to read in unsymmetric Triplet matrix    EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, false, Comm, readMap, readA, readx, 						      readb, readxexact, NonContiguousMap ) );  } else {    if ( LastFiveBytes == ".triS" ) {       NonContiguousMap = true;       // Call routine to read in symmetric Triplet matrix      EPETRA_CHK_ERR( Trilinos_Util_ReadTriples2Epetra( matrix_file, true, Comm, 							readMap, readA, readx, 							readb, readxexact, NonContiguousMap ) );    } else {      if (  LastFourBytes == ".mtx" ) { 	EPETRA_CHK_ERR( Trilinos_Util_ReadMatrixMarket2Epetra( matrix_file, Comm, readMap, 							       readA, readx, readb, readxexact) );      } else {	// Call routine to read in HB problem	Trilinos_Util_ReadHb2Epetra( matrix_file, Comm, readMap, readA, readx, 						     readb, readxexact) ;      }    }  }  Epetra_CrsMatrix transposeA(Copy, *readMap, 0);  Epetra_CrsMatrix *serialA ;   if ( transpose ) {    assert( CrsMatrixTranspose( readA, &transposeA ) == 0 );     serialA = &transposeA ;   } else {    serialA = readA ;   }  // Create uniform distributed map  Epetra_Map map(readMap->NumGlobalElements(), 0, Comm);  Epetra_Map* map_;  if( NonContiguousMap ) {    //    //  map gives us NumMyElements and MyFirstElement;    //    int NumGlobalElements =  readMap->NumGlobalElements();    int NumMyElements = map.NumMyElements();    int MyFirstElement = map.MinMyGID();    std::vector<int> MapMap_( NumGlobalElements );    readMap->MyGlobalElements( &MapMap_[0] ) ;    Comm.Broadcast( &MapMap_[0], NumGlobalElements, 0 ) ;     map_ = new Epetra_Map( NumGlobalElements, NumMyElements, &MapMap_[MyFirstElement], 0, Comm);  } else {    map_ = new Epetra_Map( map ) ;   }//.........这里部分代码省略.........

//=============================================================================int Epetra_DistObject::DoTransfer(const Epetra_SrcDistObject& A,				  Epetra_CombineMode CombineMode,				  int NumSameIDs,				  int NumPermuteIDs,				  int NumRemoteIDs,				  int NumExportIDs,				  int* PermuteToLIDs,				  int* PermuteFromLIDs,				  int* RemoteLIDs,				  int* ExportLIDs,				  int& LenExports,				  char*& Exports,				  int& LenImports,				  char*& Imports,				  Epetra_Distributor& Distor,				  bool DoReverse,                                  const Epetra_OffsetIndex * Indexor){  EPETRA_CHK_ERR(CheckSizes(A));  if (NumSameIDs + NumPermuteIDs > 0) {    EPETRA_CHK_ERR(CopyAndPermute(A, NumSameIDs, NumPermuteIDs, PermuteToLIDs, PermuteFromLIDs,Indexor, CombineMode));  }  // Once CopyAndPermute is done, switch to Add so rest works as before.  if(CombineMode == Epetra_AddLocalAlso) {    CombineMode = Add;  }  if (CombineMode==Zero)    return(0); // All done if CombineMode only involves copying and permuting  int SizeOfPacket;  bool VarSizes = false;  if( NumExportIDs > 0) {    delete [] Sizes_;    Sizes_ = new int[NumExportIDs];  }  EPETRA_CHK_ERR(PackAndPrepare(A, NumExportIDs, ExportLIDs,                 LenExports, Exports, SizeOfPacket, Sizes_, VarSizes, Distor));  if ((DistributedGlobal_ && DoReverse) || (A.Map().DistributedGlobal() && !DoReverse)) {    if (DoReverse) {      // Do the exchange of remote data      if( VarSizes ) {        EPETRA_CHK_ERR(Distor.DoReverse(Exports, SizeOfPacket, Sizes_, LenImports, Imports));      }      else {        EPETRA_CHK_ERR(Distor.DoReverse(Exports, SizeOfPacket, LenImports, Imports));      }    }    else {      // Do the exchange of remote data      if( VarSizes ) {        EPETRA_CHK_ERR(Distor.Do(Exports, SizeOfPacket, Sizes_, LenImports, Imports));      }      else {        EPETRA_CHK_ERR(Distor.Do(Exports, SizeOfPacket, LenImports, Imports));      }    }    EPETRA_CHK_ERR(UnpackAndCombine(A, NumRemoteIDs, RemoteLIDs, LenImports, Imports, SizeOfPacket, Distor, CombineMode, Indexor));  }  return(0);}

//.........这里部分代码省略.........      starts_to_ = new int[nsends_];      lengths_to_ = new int[nsends_];    }    int index = numDeadIndices;  // Leave off the dead indices (PID = -1)    int proc;    for( i = 0; i < nsends_; ++i )    {      starts_to_[i] = index;      proc = ExportPIDs[index];      procs_to_[i] = proc;      index += starts[proc];    }    if( nsends_ )      Sort_ints_( procs_to_, starts_to_, nsends_ );    max_send_length_ = 0;    for( i = 0; i < nsends_; ++i )    {      proc = procs_to_[i];      lengths_to_[i] = starts[proc];      if( (proc != my_proc) && (lengths_to_[i] > max_send_length_) )        max_send_length_ = lengths_to_[i];    }  }  else //not grouped by processor, need send buffer and indices_to_  {    if( starts[0] != 0 ) nsends_ = 1;    for( i = 1; i < nprocs; i++ )    {      if( starts[i] != 0 ) ++nsends_;      starts[i] += starts[i-1];    }    for( i = nprocs-1; i != 0; i-- )      starts[i] = starts[i-1];    starts[0] = 0;    if (nactive>0) {      indices_to_ = new int[ nactive ];      size_indices_to_ = nactive;    }    for( i = 0; i < NumExportIDs; i++ )    if( ExportPIDs[i] >= 0 )    {      indices_to_[ starts[ ExportPIDs[i] ] ] = i;      ++starts[ ExportPIDs[i] ];    }    //Reconstuct starts array to index into indices_to.    for( i = nprocs-1; i != 0; i-- )      starts[i] = starts[i-1];    starts[0] = 0;    starts[nprocs] = nactive;    if (nsends_>0) {      lengths_to_ = new int[ nsends_ ];      procs_to_ = new int[ nsends_ ];      starts_to_ = new int[ nsends_ ];    }    int j = 0;    max_send_length_ = 0;    for( i = 0; i < nprocs; i++ )      if( starts[i+1] != starts[i] )      {        lengths_to_[j] = starts[i+1] - starts[i];        starts_to_[j] = starts[i];        if( ( i != my_proc ) && ( lengths_to_[j] > max_send_length_ ) )          max_send_length_ = lengths_to_[j];        procs_to_[j] = i;        j++;      }  }      delete [] starts;  nsends_ -= self_msg_;  //Invert map to see what msgs are received and what length  EPETRA_CHK_ERR( ComputeRecvs_( my_proc, nprocs ) );  if (nrecvs_>0) {    if( !request_ ) {      request_ = new MPI_Request[ nrecvs_ ];      status_ = new MPI_Status[ nrecvs_ ];    }  }  NumRemoteIDs = total_recv_length_;  return 0;}

//==============================================================================//---------------------------------------------------------------------------//ComputeSends Method//---------------------------------------------------------------------------int Epetra_MpiDistributor::ComputeSends_( int num_imports,				const int *& import_ids,				const int *& import_procs,				int & num_exports,				int *& export_ids,				int *& export_procs,				int my_proc ) {   Epetra_MpiDistributor tmp_plan(*epComm_);  int i;  int * proc_list = 0;  int * import_objs = 0;  char * c_export_objs = 0;  if( num_imports > 0 )  {    proc_list = new int[ num_imports ];    import_objs = new int[ 2 * num_imports ];    for( i = 0; i < num_imports; i++ )    {      proc_list[i] = import_procs[i];      import_objs[2*i] = import_ids[i];      import_objs[2*i+1] = my_proc;    }  }  EPETRA_CHK_ERR(tmp_plan.CreateFromSends( num_imports, proc_list,					   true, num_exports) );  if( num_exports > 0 )  {    //export_objs = new int[ 2 * num_exports ];    export_ids = new int[ num_exports ];    export_procs = new int[ num_exports ];  }  else  {    export_ids = 0;    export_procs = 0;  }  int len_c_export_objs = 0;  EPETRA_CHK_ERR( tmp_plan.Do(reinterpret_cast<char *> (import_objs),			      2 * (int)sizeof( int ), 			      len_c_export_objs,			      c_export_objs) );  int * export_objs = reinterpret_cast<int *>(c_export_objs);  for( i = 0; i < num_exports; i++ ) {    export_ids[i] = export_objs[2*i];    export_procs[i] = export_objs[2*i+1];  }  if( proc_list != 0 ) delete [] proc_list;  if( import_objs != 0 ) delete [] import_objs;  if( len_c_export_objs != 0 ) delete [] c_export_objs;  return(0);}

//==============================================================================int LinearProblem_CrsSingletonFilter::ConstructReducedProblem(Epetra_LinearProblem * Problem) {  int i, j;  if (HaveReducedProblem_) EPETRA_CHK_ERR(-1); // Setup already done once.  Cannot do it again  if (Problem==0) EPETRA_CHK_ERR(-2); // Null problem pointer  FullProblem_ = Problem;  FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());  if (FullMatrix_==0) EPETRA_CHK_ERR(-3); // Need a RowMatrix  if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-4); // Need a RHS  if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-5); // Need a LHS  // Generate reduced row and column maps  Epetra_MapColoring & RowMapColors = *RowMapColors_;  Epetra_MapColoring & ColMapColors = *ColMapColors_;  ReducedMatrixRowMap_ = RowMapColors.GenerateMap(0);  ReducedMatrixColMap_ = ColMapColors.GenerateMap(0);  // Create domain and range map colorings by exporting map coloring of column and row maps  if (FullMatrix()->RowMatrixImporter()!=0) {    Epetra_MapColoring DomainMapColors(FullMatrixDomainMap());    EPETRA_CHK_ERR(DomainMapColors.Export(*ColMapColors_, *FullMatrix()->RowMatrixImporter(), AbsMax));    OrigReducedMatrixDomainMap_ = DomainMapColors.GenerateMap(0);  }  else    OrigReducedMatrixDomainMap_ = ReducedMatrixColMap_;  if (FullMatrixIsCrsMatrix_) {    if (FullCrsMatrix()->Exporter()!=0) { // Non-trivial exporter      Epetra_MapColoring RangeMapColors(FullMatrixRangeMap());      EPETRA_CHK_ERR(RangeMapColors.Export(*RowMapColors_, *FullCrsMatrix()->Exporter(), 					   AbsMax));      ReducedMatrixRangeMap_ = RangeMapColors.GenerateMap(0);    }    else      ReducedMatrixRangeMap_ = ReducedMatrixRowMap_;  }  else    ReducedMatrixRangeMap_ = ReducedMatrixRowMap_;  // Check to see if the reduced system domain and range maps are the same.  // If not, we need to remap entries of the LHS multivector so that they are distributed  // conformally with the rows of the reduced matrix and the RHS multivector  SymmetricElimination_ = ReducedMatrixRangeMap_->SameAs(*OrigReducedMatrixDomainMap_);  if (!SymmetricElimination_)     ConstructRedistributeExporter(OrigReducedMatrixDomainMap_, ReducedMatrixRangeMap_, 				  RedistributeDomainExporter_, ReducedMatrixDomainMap_);  else {    ReducedMatrixDomainMap_ = OrigReducedMatrixDomainMap_;    OrigReducedMatrixDomainMap_ = 0;    RedistributeDomainExporter_ = 0;  }    // Create pointer to Full RHS, LHS  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();  int NumVectors = FullLHS->NumVectors();  // Create importers//  cout << "RedDomainMap/n";//  cout << *ReducedMatrixDomainMap();//  cout << "FullDomainMap/n";//  cout << FullMatrixDomainMap();  Full2ReducedLHSImporter_ = new Epetra_Import(*ReducedMatrixDomainMap(), FullMatrixDomainMap());//  cout << "RedRowMap/n";//  cout << *ReducedMatrixRowMap();//  cout << "FullRHSMap/n";//  cout << FullRHS->Map();  Full2ReducedRHSImporter_ = new Epetra_Import(*ReducedMatrixRowMap(), FullRHS->Map());  // Construct Reduced Matrix  ReducedMatrix_ = new Epetra_CrsMatrix(Copy, *ReducedMatrixRowMap(), *ReducedMatrixColMap(), 0);  // Create storage for temporary X values due to explicit elimination of rows  tempExportX_ = new Epetra_MultiVector(FullMatrixColMap(), NumVectors);  int NumEntries;  int * Indices;  double * Values;  int NumMyRows = FullMatrix()->NumMyRows();  int ColSingletonCounter = 0;  for (i=0; i<NumMyRows; i++) {    int curGRID = FullMatrixRowMap().GID(i);    if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix      EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (Indices are global)            int ierr = ReducedMatrix()->InsertGlobalValues(curGRID, NumEntries, 						     Values, Indices); // Insert into reduce matrix      // Positive errors will occur because we are submitting col entries that are not part of      // reduced system.  However, because we specified a column map to the ReducedMatrix constructor      // these extra column entries will be ignored and we will be politely reminded by a positive      // error code      if (ierr<0) EPETRA_CHK_ERR(ierr);     }    else {      EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row//.........这里部分代码省略.........

//==============================================================================int LinearProblem_CrsSingletonFilter::UpdateReducedProblem(Epetra_LinearProblem * Problem) {  int i, j;  if (Problem==0) EPETRA_CHK_ERR(-1); // Null problem pointer  FullProblem_ = Problem;  FullMatrix_ = dynamic_cast<Epetra_RowMatrix *>(Problem->GetMatrix());  if (FullMatrix_==0) EPETRA_CHK_ERR(-2); // Need a RowMatrix  if (Problem->GetRHS()==0) EPETRA_CHK_ERR(-3); // Need a RHS  if (Problem->GetLHS()==0) EPETRA_CHK_ERR(-4); // Need a LHS  if (!HaveReducedProblem_) EPETRA_CHK_ERR(-5); // Must have set up reduced problem  // Create pointer to Full RHS, LHS  Epetra_MultiVector * FullRHS = FullProblem()->GetRHS();  Epetra_MultiVector * FullLHS = FullProblem()->GetLHS();  int NumVectors = FullLHS->NumVectors();  int NumEntries;  int * Indices;  double * Values;  int NumMyRows = FullMatrix()->NumMyRows();  int ColSingletonCounter = 0;  for (i=0; i<NumMyRows; i++) {    int curGRID = FullMatrixRowMap().GID(i);    if (ReducedMatrixRowMap()->MyGID(curGRID)) { // Check if this row should go into reduced matrix      EPETRA_CHK_ERR(GetRowGCIDs(i, NumEntries, Values, Indices)); // Get current row (indices global)      int ierr = ReducedMatrix()->ReplaceGlobalValues(curGRID, NumEntries, 						      Values, Indices);      // Positive errors will occur because we are submitting col entries that are not part of      // reduced system.  However, because we specified a column map to the ReducedMatrix constructor      // these extra column entries will be ignored and we will be politely reminded by a positive      // error code      if (ierr<0) EPETRA_CHK_ERR(ierr);     }    // Otherwise if singleton row we explicitly eliminate this row and solve for corresponding X value    else {      EPETRA_CHK_ERR(GetRow(i, NumEntries, Values, Indices)); // Get current row      if (NumEntries==1) {	double pivot = Values[0];	if (pivot==0.0) EPETRA_CHK_ERR(-1); // Encountered zero row, unable to continue	int indX = Indices[0];	for (j=0; j<NumVectors; j++)	  (*tempExportX_)[j][indX] = (*FullRHS)[j][i]/pivot;      }      // Otherwise, this is a singleton column and we will scan for the pivot element needed       // for post-solve equations      else {	j = ColSingletonPivotLIDs_[ColSingletonCounter];	double pivot = Values[j];	if (pivot==0.0) EPETRA_CHK_ERR(-2); // Encountered zero column, unable to continue	ColSingletonPivots_[ColSingletonCounter] = pivot;	ColSingletonCounter++;      }    }  }  assert(ColSingletonCounter==NumMyColSingletons_); // Sanity test  // Update Reduced LHS (Puts any initial guess values into reduced system)  ReducedLHS_->PutScalar(0.0); // zero out Reduced LHS  EPETRA_CHK_ERR(ReducedLHS_->Import(*FullLHS, *Full2ReducedLHSImporter_, Insert));  FullLHS->PutScalar(0.0); // zero out Full LHS since we will inject values as we get them  // Construct Reduced RHS  // Zero out temp space  tempX_->PutScalar(0.0);  tempB_->PutScalar(0.0);    //Inject known X values into tempX for purpose of computing tempB = FullMatrix*tempX  // Also inject into full X since we already know the solution  if (FullMatrix()->RowMatrixImporter()!=0) {    EPETRA_CHK_ERR(tempX_->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));    EPETRA_CHK_ERR(FullLHS->Export(*tempExportX_, *FullMatrix()->RowMatrixImporter(), Add));  }  else {    tempX_->Update(1.0, *tempExportX_, 0.0);    FullLHS->Update(1.0, *tempExportX_, 0.0);  }  EPETRA_CHK_ERR(FullMatrix()->Multiply(false, *tempX_, *tempB_));  EPETRA_CHK_ERR(tempB_->Update(1.0, *FullRHS, -1.0)); // tempB now has influence of already-known X values  ReducedRHS_->PutScalar(0.0);  EPETRA_CHK_ERR(ReducedRHS_->Import(*tempB_, *Full2ReducedRHSImporter_, Insert));    return(0);}

int Ifpack_CrsRiluk::InitAllValues(const Epetra_RowMatrix & OverlapA, int MaxNumEntries) {  int ierr = 0;  int i, j;  int NumIn, NumL, NumU;  bool DiagFound;  int NumNonzeroDiags = 0;  vector<int> InI(MaxNumEntries); // Allocate temp space  vector<int> LI(MaxNumEntries);  vector<int> UI(MaxNumEntries);  vector<double> InV(MaxNumEntries);  vector<double> LV(MaxNumEntries);  vector<double> UV(MaxNumEntries);  bool ReplaceValues = (L_->StaticGraph() || L_->IndicesAreLocal()); // Check if values should be inserted or replaced  if (ReplaceValues) {    L_->PutScalar(0.0); // Zero out L and U matrices    U_->PutScalar(0.0);  }  D_->PutScalar(0.0); // Set diagonal values to zero  double *DV;  EPETRA_CHK_ERR(D_->ExtractView(&DV)); // Get view of diagonal      // First we copy the user's matrix into L and U, regardless of fill level  for (i=0; i< NumMyRows(); i++) {    EPETRA_CHK_ERR(OverlapA.ExtractMyRowCopy(i, MaxNumEntries, NumIn, &InV[0], &InI[0])); // Get Values and Indices        // Split into L and U (we don't assume that indices are ordered).        NumL = 0;     NumU = 0;     DiagFound = false;        for (j=0; j< NumIn; j++) {      int k = InI[j];      if (k==i) {	DiagFound = true;	DV[i] += Rthresh_ * InV[j] + EPETRA_SGN(InV[j]) * Athresh_; // Store perturbed diagonal in Epetra_Vector D_      }      else if (k < 0) {EPETRA_CHK_ERR(-1);} // Out of range      else if (k < i) {	LI[NumL] = k;	LV[NumL] = InV[j];	NumL++;      }      else if (k<NumMyRows()) {	UI[NumU] = k;	UV[NumU] = InV[j];	NumU++;      }    }        // Check in things for this row of L and U    if (DiagFound) NumNonzeroDiags++;    else DV[i] = Athresh_;    if (NumL) {      if (ReplaceValues) {	EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, &LV[0], &LI[0]));      }      else {	EPETRA_CHK_ERR(L_->InsertMyValues(i, NumL, &LV[0], &LI[0]));      }    }    if (NumU) {      if (ReplaceValues) {	EPETRA_CHK_ERR(U_->ReplaceMyValues(i, NumU, &UV[0], &UI[0]));      }      else {	EPETRA_CHK_ERR(U_->InsertMyValues(i, NumU, &UV[0], &UI[0]));      }    }      }  if (!ReplaceValues) {    // The domain of L and the range of U are exactly their own row maps (there is no communication).    // The domain of U and the range of L must be the same as those of the original matrix,    // However if the original matrix is a VbrMatrix, these two latter maps are translation from    // a block map to a point map.    EPETRA_CHK_ERR(L_->FillComplete(L_->RowMatrixColMap(), *L_RangeMap_));    EPETRA_CHK_ERR(U_->FillComplete(*U_DomainMap_, U_->RowMatrixRowMap()));  }  // At this point L and U have the values of A in the structure of L and U, and diagonal vector D  SetValuesInitialized(true);  SetFactored(false);//.........这里部分代码省略.........

int runTests(Epetra_Map & map, Epetra_CrsMatrix & A, Epetra_Vector & x, Epetra_Vector & b, Epetra_Vector & xexact, bool verbose) {  int ierr = 0;  // Create MultiVectors and put x, b, xexact in both columns of X, B, and Xexact, respectively.  Epetra_MultiVector X( map, 2, false );  Epetra_MultiVector B( map, 2, false );  Epetra_MultiVector Xexact( map, 2, false );  for (int i=0; i<X.NumVectors(); ++i) {    *X(i) = x;    *B(i) = b;    *Xexact(i) = xexact;  }  double residual;  std::vector<double> residualmv(2);  residual = A.NormInf(); double rAInf = residual;  if (verbose) std::cout << "Inf Norm of A                                                     = " << residual << std::endl;  residual = A.NormOne(); double rAOne = residual;  if (verbose) std::cout << "One Norm of A                                                     = " << residual << std::endl;  xexact.Norm2(&residual); double rxx = residual;	  Xexact.Norm2(&residualmv[0]); std::vector<double> rXX( residualmv );	  if (verbose) std::cout << "Norm of xexact                                                    = " << residual << std::endl;  if (verbose) std::cout << "Norm of Xexact                                                    = (" << residualmv[0] << ", " <<residualmv[1] <<")"<< std::endl;  Epetra_Vector tmp1(map);  Epetra_MultiVector tmp1mv(map,2,false);  A.Multiply(false, xexact, tmp1);  A.Multiply(false, Xexact, tmp1mv);  tmp1.Norm2(&residual); double rAx = residual;  tmp1mv.Norm2(&residualmv[0]); std::vector<double> rAX( residualmv );  if (verbose) std::cout << "Norm of Ax                                                        = " << residual << std::endl;  if (verbose) std::cout << "Norm of AX                                                        = (" << residualmv[0] << ", " << residualmv[1] <<")"<< std::endl;  b.Norm2(&residual); double rb = residual;  B.Norm2(&residualmv[0]); std::vector<double> rB( residualmv );  if (verbose) std::cout << "Norm of b (should equal norm of Ax)                               = " << residual << std::endl;  if (verbose) std::cout << "Norm of B (should equal norm of AX)                               = (" << residualmv[0] << ", " << residualmv[1] <<")"<< std::endl;  tmp1.Update(1.0, b, -1.0);  tmp1mv.Update(1.0, B, -1.0);  tmp1.Norm2(&residual);  tmp1mv.Norm2(&residualmv[0]);  if (verbose) std::cout << "Norm of difference between compute Ax and Ax from file            = " << residual << std::endl;  if (verbose) std::cout << "Norm of difference between compute AX and AX from file            = (" << residualmv[0] << ", " << residualmv[1] <<")"<< std::endl;  map.Comm().Barrier();  EPETRA_CHK_ERR(EpetraExt::BlockMapToMatrixMarketFile("Test_map.mm", map, "Official EpetraExt test map", 						       "This is the official EpetraExt test map generated by the EpetraExt regression tests"));  EPETRA_CHK_ERR(EpetraExt::RowMatrixToMatrixMarketFile("Test_A.mm", A, "Official EpetraExt test matrix", 							"This is the official EpetraExt test matrix generated by the EpetraExt regression tests"));  EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_x.mm", x, "Official EpetraExt test initial guess", 						     "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));  EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvX.mm", X, "Official EpetraExt test initial guess", 					      	          "This is the official EpetraExt test initial guess generated by the EpetraExt regression tests"));				         EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_xexact.mm", xexact, "Official EpetraExt test exact solution", 						     "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));  EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvXexact.mm", Xexact, "Official EpetraExt test exact solution", 						          "This is the official EpetraExt test exact solution generated by the EpetraExt regression tests"));				         EPETRA_CHK_ERR(EpetraExt::VectorToMatrixMarketFile("Test_b.mm", b, "Official EpetraExt test right hand side", 						     "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));  EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatrixMarketFile("Test_mvB.mm", B, "Official EpetraExt test right hand side", 						          "This is the official EpetraExt test right hand side generated by the EpetraExt regression tests"));				         EPETRA_CHK_ERR(EpetraExt::MultiVectorToMatlabFile("Test_mvB.mat", B));   EPETRA_CHK_ERR(EpetraExt::RowMatrixToMatlabFile("Test_A.dat", A));  Epetra_Map * map1;  Epetra_CrsMatrix * A1;   Epetra_CrsMatrix * A2;   Epetra_CrsMatrix * A3;   Epetra_Vector * x1;   Epetra_Vector * b1;  Epetra_Vector * xexact1;  Epetra_MultiVector * X1;   Epetra_MultiVector * B1;  Epetra_MultiVector * Xexact1;  EpetraExt::MatrixMarketFileToMap("Test_map.mm", map.Comm(), map1);  if (map.SameAs(*map1)) {    if (verbose) std::cout << "Maps are equal.  In/Out works." << std::endl;  }  else {    if (verbose) std::cout << "Maps are not equal.  In/Out fails." << std::endl;    ierr += 1;  }  EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", *map1, A1));  // If map is zero-based, then we can compare to the convenient reading versions  if (map1->IndexBase()==0) EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToCrsMatrix("Test_A.mm", map1->Comm(), A2));  if (map1->IndexBase()==0) EPETRA_CHK_ERR(EpetraExt::MatlabFileToCrsMatrix("Test_A.dat", map1->Comm(), A3));  EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_x.mm", *map1, x1));  EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_xexact.mm", *map1, xexact1));  EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToVector("Test_b.mm", *map1, b1));  EPETRA_CHK_ERR(EpetraExt::MatrixMarketFileToMultiVector("Test_mvX.mm", *map1, X1));//.........这里部分代码省略.........

//==========================================================================int Ifpack_CrsRiluk::Factor() {  // if (!Allocated()) return(-1); // This test is not needed at this time.  All constructors allocate.  if (!ValuesInitialized()) return(-2); // Must have values initialized.  if (Factored()) return(-3); // Can't have already computed factors.  SetValuesInitialized(false);  // MinMachNum should be officially defined, for now pick something a little   // bigger than IEEE underflow value  double MinDiagonalValue = Epetra_MinDouble;  double MaxDiagonalValue = 1.0/MinDiagonalValue;  int ierr = 0;  int i, j, k;  int * LI=0, * UI = 0;  double * LV=0, * UV = 0;  int NumIn, NumL, NumU;  // Get Maximun Row length  int MaxNumEntries = L_->MaxNumEntries() + U_->MaxNumEntries() + 1;  vector<int> InI(MaxNumEntries); // Allocate temp space  vector<double> InV(MaxNumEntries);  vector<int> colflag(NumMyCols());  double *DV;  ierr = D_->ExtractView(&DV); // Get view of diagonal#ifdef IFPACK_FLOPCOUNTERS  int current_madds = 0; // We will count multiply-add as they happen#endif  // Now start the factorization.  // Need some integer workspace and pointers  int NumUU;   int * UUI;  double * UUV;  for (j=0; j<NumMyCols(); j++) colflag[j] = - 1;  for(i=0; i<NumMyRows(); i++) { // Fill InV, InI with current row of L, D and U combined    NumIn = MaxNumEntries;    EPETRA_CHK_ERR(L_->ExtractMyRowCopy(i, NumIn, NumL, &InV[0], &InI[0]));    LV = &InV[0];    LI = &InI[0];    InV[NumL] = DV[i]; // Put in diagonal    InI[NumL] = i;        EPETRA_CHK_ERR(U_->ExtractMyRowCopy(i, NumIn-NumL-1, NumU, &InV[NumL+1], &InI[NumL+1]));    NumIn = NumL+NumU+1;    UV = &InV[NumL+1];    UI = &InI[NumL+1];    // Set column flags    for (j=0; j<NumIn; j++) colflag[InI[j]] = j;    double diagmod = 0.0; // Off-diagonal accumulator    for (int jj=0; jj<NumL; jj++) {      j = InI[jj];      double multiplier = InV[jj]; // current_mults++;      InV[jj] *= DV[j];            EPETRA_CHK_ERR(U_->ExtractMyRowView(j, NumUU, UUV, UUI)); // View of row above      if (RelaxValue_==0.0) {	for (k=0; k<NumUU; k++) {	  int kk = colflag[UUI[k]];	  if (kk>-1) {	    InV[kk] -= multiplier*UUV[k];#ifdef IFPACK_FLOPCOUNTERS	    current_madds++;#endif	  }	}      }      else {	for (k=0; k<NumUU; k++) {	  int kk = colflag[UUI[k]];	  if (kk>-1) InV[kk] -= multiplier*UUV[k];	  else diagmod -= multiplier*UUV[k];#ifdef IFPACK_FLOPCOUNTERS	  current_madds++;#endif	}      }     }    if (NumL) {      EPETRA_CHK_ERR(L_->ReplaceMyValues(i, NumL, LV, LI));  // Replace current row of L    }    DV[i] = InV[NumL]; // Extract Diagonal value//.........这里部分代码省略.........

int main(int argc, char *argv[]) {#ifdef EPETRA_MPI  MPI_Init(&argc,&argv);  Epetra_MpiComm comm (MPI_COMM_WORLD);#else  Epetra_SerialComm comm;#endif  int MyPID = comm.MyPID();  bool verbose = false;  bool verbose1 = false;   // Check if we should print results to standard out  if (argc > 1) {    if ((argv[1][0] == '-') && (argv[1][1] == 'v')) {      verbose1 = true;      if (MyPID==0) verbose = true;    }  }  if (verbose)    std::cout << EpetraExt::EpetraExt_Version() << std::endl << std::endl;  if (verbose1) std::cout << comm << std::endl;  // Uncomment the next three lines to debug in mpi mode  //int tmp;  //if (MyPID==0) cin >> tmp;  //comm.Barrier();  Epetra_Map * map;  Epetra_CrsMatrix * A;   Epetra_Vector * x;   Epetra_Vector * b;  Epetra_Vector * xexact;  int nx = 20*comm.NumProc();  int ny = 30;  int npoints = 7;  int xoff[] = {-1,  0,  1, -1,  0,  1,  0};  int yoff[] = {-1, -1, -1,  0,  0,  0,  1};     int ierr = 0;  // Call routine to read in HB problem 0-base  Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact);  ierr += runTests(*map, *A, *x, *b, *xexact, verbose);  delete A;  delete x;  delete b;  delete xexact;  delete map;  // Call routine to read in HB problem 1-base  Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, 1);  ierr += runTests(*map, *A, *x, *b, *xexact, verbose);  delete A;  delete x;  delete b;  delete xexact;  delete map;  // Call routine to read in HB problem -1-base  Trilinos_Util_GenerateCrsProblem(nx, ny, npoints, xoff, yoff, comm, map, A, x, b, xexact, -1);  ierr += runTests(*map, *A, *x, *b, *xexact, verbose);  delete A;  delete x;  delete b;  delete xexact;  delete map;  int nx1 = 5;  int ny1 = 4;  Poisson2dOperator Op(nx1, ny1, comm);  ierr += runOperatorTests(Op, verbose);  generateHyprePrintOut("MyMatrixFile", comm);  EPETRA_CHK_ERR(EpetraExt::HypreFileToCrsMatrix("MyMatrixFile", comm, A));    runHypreTest(*A);  delete A;  #ifdef EPETRA_MPI  MPI_Finalize() ;#endif  return(ierr);}