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);
}