C++ Epetra_RowMatrix::NumGlobalNonzeros方法代码示例

void ModeLaplace1DQ1::memoryInfo() const {

  int myPid = MyComm.MyPID();

  Epetra_RowMatrix *Mat = dynamic_cast<Epetra_RowMatrix *>(M);
  if ((myPid == 0) && (Mat)) {
    cout << " Total number of nonzero entries in mass matrix      = ";
    cout.width(15);
    cout << Mat->NumGlobalNonzeros() << endl;
    double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
    memSize += 2*Mat->NumGlobalRows()*sizeof(int);
    cout << " Memory requested for mass matrix per processor      = (EST) ";
    cout.precision(2);
    cout.width(6);
    cout.setf(ios::fixed, ios::floatfield);
    cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
    cout << endl;
  }

  Mat = dynamic_cast<Epetra_RowMatrix *>(K);
  if ((myPid == 0) && (Mat)) {
    cout << " Total number of nonzero entries in stiffness matrix = ";
    cout.width(15);
    cout << Mat->NumGlobalNonzeros() << endl;
    double memSize = Mat->NumGlobalNonzeros()*(sizeof(double) + sizeof(int));
    memSize += 2*Mat->NumGlobalRows()*sizeof(int);
    cout << " Memory requested for stiffness matrix per processor = (EST) ";
    cout.precision(2);
    cout.width(6);
    cout.setf(ios::fixed, ios::floatfield);
    cout << memSize/1024.0/1024.0/MyComm.NumProc() << " MB " << endl;
    cout << endl;
  }

}

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

  int Rows = Matrix.NumGlobalRows();
  int Cols = Matrix.NumGlobalRows();
  int Nonzeros = Matrix.NumGlobalNonzeros();

  if (Comm_.MyPID() == 0)
  {
    std::ofstream of(FileName_.c_str(), std::ios::app);
    of << "<PointMatrix Label=\"" << Label << '"'
      << " Rows=\"" << Rows << '"'
      << " Columns=\"" << Cols<< '"'
      << " Nonzeros=\"" << Nonzeros << '"'
      << " Type=\"double\" StartingIndex=\"0\">" << std::endl;
  }

  int Length = Matrix.MaxNumEntries();
  std::vector<int> Indices(Length);
  std::vector<double> Values(Length);

  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 < Matrix.NumMyRows(); ++i)
      {
        int NumMyEntries;
        Matrix.ExtractMyRowCopy(i, Length, NumMyEntries, &Values[0], &Indices[0]);

        int GRID = Matrix.RowMatrixRowMap().GID(i);

        for (int j = 0; j < NumMyEntries; ++j)
          of << GRID << " " << Matrix.RowMatrixColMap().GID(Indices[j])
             << " " << std::setiosflags(std::ios::scientific) << Values[j] << std::endl;
      }
      of.close();
    }
    Comm_.Barrier();
  }

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

int Amesos_Scalapack::NumericFactorization() {
  
  if( debug_ == 1 ) std::cout << "Entering `NumericFactorization()'" << std::endl;
  
  NumNumericFact_++;
  
  iam_ = Comm().MyPID();
  
  Epetra_RowMatrix *RowMatrixA = dynamic_cast<Epetra_RowMatrix *>(Problem_->GetOperator());
  const Epetra_Map &OriginalMap = RowMatrixA->RowMatrixRowMap() ; 
  NumGlobalElements_ = OriginalMap.NumGlobalElements();

  NumGlobalNonzeros_ = RowMatrixA->NumGlobalNonzeros();
  
  RedistributeA();
  ConvertToScalapack();
  
  return PerformNumericFactorization( );
}

//=============================================================================
int Amesos_Dscpack::PerformSymbolicFactorization()
{
  ResetTimer(0);
  ResetTimer(1);

  MyPID_    = Comm().MyPID();
  NumProcs_ = Comm().NumProc();
  
  Epetra_RowMatrix *RowMatrixA = Problem_->GetMatrix();
  if (RowMatrixA == 0)
    AMESOS_CHK_ERR(-1);

  const Epetra_Map& OriginalMap = RowMatrixA->RowMatrixRowMap() ;
  const Epetra_MpiComm& comm1   = dynamic_cast<const Epetra_MpiComm &> (Comm());
  int numrows                   = RowMatrixA->NumGlobalRows();
  int numentries                = RowMatrixA->NumGlobalNonzeros();

  Teuchos::RCP<Epetra_CrsGraph> Graph;

  Epetra_CrsMatrix* CastCrsMatrixA = 
    dynamic_cast<Epetra_CrsMatrix*>(RowMatrixA); 

  if (CastCrsMatrixA)
  {
    Graph = Teuchos::rcp(const_cast<Epetra_CrsGraph*>(&(CastCrsMatrixA->Graph())), false);
  }
  else
  {
    int MaxNumEntries = RowMatrixA->MaxNumEntries();
    Graph = Teuchos::rcp(new Epetra_CrsGraph(Copy, OriginalMap, MaxNumEntries));

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

    for (int i = 0 ; i < RowMatrixA->NumMyRows() ; ++i)
    {
      int NumEntries;
      RowMatrixA->ExtractMyRowCopy(i, MaxNumEntries, NumEntries,
                                   &Values[0], &Indices[0]);

      for (int j = 0 ; j < NumEntries ; ++j)
        Indices[j] = RowMatrixA->RowMatrixColMap().GID(Indices[j]);

      int GlobalRow = RowMatrixA->RowMatrixRowMap().GID(i);
      Graph->InsertGlobalIndices(GlobalRow, NumEntries, &Indices[0]);
    }

    Graph->FillComplete();
  }

  //
  //  Create a replicated map and graph 
  //
  std::vector<int> AllIDs( numrows ) ; 
  for ( int i = 0; i < numrows ; i++ ) AllIDs[i] = i ; 

  Epetra_Map      ReplicatedMap( -1, numrows, &AllIDs[0], 0, Comm());
  Epetra_Import   ReplicatedImporter(ReplicatedMap, OriginalMap);
  Epetra_CrsGraph ReplicatedGraph( Copy, ReplicatedMap, 0 ); 

  AMESOS_CHK_ERR(ReplicatedGraph.Import(*Graph, ReplicatedImporter, Insert));
  AMESOS_CHK_ERR(ReplicatedGraph.FillComplete());

  //
  //  Convert the matrix to Ap, Ai
  //
  std::vector <int> Replicates(numrows);
  std::vector <int> Ap(numrows + 1);
  std::vector <int> Ai(EPETRA_MAX(numrows, numentries));

  for( int i = 0 ; i < numrows; i++ ) Replicates[i] = 1; 
  
  int NumEntriesPerRow ;
  int *ColIndices = 0 ;
  int Ai_index = 0 ; 
  for ( int MyRow = 0; MyRow <numrows; MyRow++ ) {
    AMESOS_CHK_ERR( ReplicatedGraph.ExtractMyRowView( MyRow, NumEntriesPerRow, ColIndices ) );
    Ap[MyRow] = Ai_index ; 
    for ( int j = 0; j < NumEntriesPerRow; j++ ) { 
      Ai[Ai_index] = ColIndices[j] ; 
      Ai_index++;
    }
  }
  assert( Ai_index == numentries ) ; 
  Ap[ numrows ] = Ai_index ; 
  
  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_, 0);

  ResetTimer(0);

  //
  //  Call Dscpack Symbolic Factorization
  //  
  int OrderCode = 2;
  std::vector<double> MyANonZ;
  
  NumLocalNonz = 0 ; 
  GlobalStructNewColNum = 0 ; 
  GlobalStructNewNum = 0 ;  
//.........这里部分代码省略.........

int check(Epetra_RowMatrix& A, Epetra_RowMatrix & B, bool verbose)  {

    int ierr = 0;
    EPETRA_TEST_ERR(!A.Comm().NumProc()==B.Comm().NumProc(),ierr);
    EPETRA_TEST_ERR(!A.Comm().MyPID()==B.Comm().MyPID(),ierr);
    EPETRA_TEST_ERR(!A.Filled()==B.Filled(),ierr);
    EPETRA_TEST_ERR(!A.HasNormInf()==B.HasNormInf(),ierr);
    EPETRA_TEST_ERR(!A.LowerTriangular()==B.LowerTriangular(),ierr);
    EPETRA_TEST_ERR(!A.Map().SameAs(B.Map()),ierr);
    EPETRA_TEST_ERR(!A.MaxNumEntries()==B.MaxNumEntries(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalCols()==B.NumGlobalCols(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalDiagonals()==B.NumGlobalDiagonals(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalNonzeros()==B.NumGlobalNonzeros(),ierr);
    EPETRA_TEST_ERR(!A.NumGlobalRows()==B.NumGlobalRows(),ierr);
    EPETRA_TEST_ERR(!A.NumMyCols()==B.NumMyCols(),ierr);
    EPETRA_TEST_ERR(!A.NumMyDiagonals()==B.NumMyDiagonals(),ierr);
    EPETRA_TEST_ERR(!A.NumMyNonzeros()==B.NumMyNonzeros(),ierr);
    for (int i=0; i<A.NumMyRows(); i++) {
        int nA, nB;
        A.NumMyRowEntries(i,nA);
        B.NumMyRowEntries(i,nB);
        EPETRA_TEST_ERR(!nA==nB,ierr);
    }
    EPETRA_TEST_ERR(!A.NumMyRows()==B.NumMyRows(),ierr);
    EPETRA_TEST_ERR(!A.OperatorDomainMap().SameAs(B.OperatorDomainMap()),ierr);
    EPETRA_TEST_ERR(!A.OperatorRangeMap().SameAs(B.OperatorRangeMap()),ierr);
    EPETRA_TEST_ERR(!A.RowMatrixColMap().SameAs(B.RowMatrixColMap()),ierr);
    EPETRA_TEST_ERR(!A.RowMatrixRowMap().SameAs(B.RowMatrixRowMap()),ierr);
    EPETRA_TEST_ERR(!A.UpperTriangular()==B.UpperTriangular(),ierr);
    EPETRA_TEST_ERR(!A.UseTranspose()==B.UseTranspose(),ierr);

    int NumVectors = 5;
    {   // No transpose case
        Epetra_MultiVector X(A.OperatorDomainMap(), NumVectors);
        Epetra_MultiVector YA1(A.OperatorRangeMap(), NumVectors);
        Epetra_MultiVector YA2(YA1);
        Epetra_MultiVector YB1(YA1);
        Epetra_MultiVector YB2(YA1);
        X.Random();

        bool transA = false;
        A.SetUseTranspose(transA);
        B.SetUseTranspose(transA);
        A.Apply(X,YA1);
        A.Multiply(transA, X, YA2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2,"A Multiply and A Apply", verbose),ierr);
        B.Apply(X,YB1);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1,"A Multiply and B Multiply", verbose),ierr);
        B.Multiply(transA, X, YB2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2,"A Multiply and B Apply", verbose), ierr);

    }
    {   // transpose case
        Epetra_MultiVector X(A.OperatorRangeMap(), NumVectors);
        Epetra_MultiVector YA1(A.OperatorDomainMap(), NumVectors);
        Epetra_MultiVector YA2(YA1);
        Epetra_MultiVector YB1(YA1);
        Epetra_MultiVector YB2(YA1);
        X.Random();

        bool transA = true;
        A.SetUseTranspose(transA);
        B.SetUseTranspose(transA);
        A.Apply(X,YA1);
        A.Multiply(transA, X, YA2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YA2, "A Multiply and A Apply (transpose)", verbose),ierr);
        B.Apply(X,YB1);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB1, "A Multiply and B Multiply (transpose)", verbose),ierr);
        B.Multiply(transA, X,YB2);
        EPETRA_TEST_ERR(checkMultiVectors(YA1,YB2, "A Multiply and B Apply (transpose)", verbose),ierr);

    }

    Epetra_Vector diagA(A.RowMatrixRowMap());
    EPETRA_TEST_ERR(A.ExtractDiagonalCopy(diagA),ierr);
    Epetra_Vector diagB(B.RowMatrixRowMap());
    EPETRA_TEST_ERR(B.ExtractDiagonalCopy(diagB),ierr);
    EPETRA_TEST_ERR(checkMultiVectors(diagA,diagB, "ExtractDiagonalCopy", verbose),ierr);

    Epetra_Vector rowA(A.RowMatrixRowMap());
    EPETRA_TEST_ERR(A.InvRowSums(rowA),ierr);
    Epetra_Vector rowB(B.RowMatrixRowMap());
    EPETRA_TEST_ERR(B.InvRowSums(rowB),ierr)
    EPETRA_TEST_ERR(checkMultiVectors(rowA,rowB, "InvRowSums", verbose),ierr);

    Epetra_Vector colA(A.RowMatrixColMap());
    EPETRA_TEST_ERR(A.InvColSums(colA),ierr);
    Epetra_Vector colB(B.RowMatrixColMap());
    EPETRA_TEST_ERR(B.InvColSums(colB),ierr);
    EPETRA_TEST_ERR(checkMultiVectors(colA,colB, "InvColSums", verbose),ierr);

    EPETRA_TEST_ERR(checkValues(A.NormInf(), B.NormInf(), "NormInf before scaling", verbose), ierr);
    EPETRA_TEST_ERR(checkValues(A.NormOne(), B.NormOne(), "NormOne before scaling", verbose),ierr);

    EPETRA_TEST_ERR(A.RightScale(colA),ierr);
    EPETRA_TEST_ERR(B.RightScale(colB),ierr);


    EPETRA_TEST_ERR(A.LeftScale(rowA),ierr);
    EPETRA_TEST_ERR(B.LeftScale(rowB),ierr);
//.........这里部分代码省略.........

//=============================================================================
int Amesos_Mumps::ConvertToTriplet(const bool OnlyValues)
{

  Epetra_RowMatrix* ptr;
  if (Comm().NumProc() == 1)
    ptr = &Matrix();
  else {
    ptr = &RedistrMatrix(true);
  }

  ResetTimer();
  
#ifdef EXTRA_DEBUG_INFO
  Epetra_CrsMatrix* Eptr = dynamic_cast<Epetra_CrsMatrix*>( ptr );
  if ( ptr->NumGlobalNonzeros() < 300 ) SetICNTL(4,3 );  // Enable more debug info for small matrices
  if ( ptr->NumGlobalNonzeros() < 42 && Eptr ) { 
      std::cout << " Matrix = " << std::endl ; 
      Eptr->Print( std::cout ) ; 
  } else {
      assert( Eptr );
  }
#endif

  Row.resize(ptr->NumMyNonzeros());
  Col.resize(ptr->NumMyNonzeros());
  Val.resize(ptr->NumMyNonzeros());

  int MaxNumEntries = ptr->MaxNumEntries();
  std::vector<int> Indices;
  std::vector<double> Values;
  Indices.resize(MaxNumEntries);
  Values.resize(MaxNumEntries);

  int count = 0;

  for (int i = 0; i < ptr->NumMyRows() ; ++i) {

    int GlobalRow = ptr->RowMatrixRowMap().GID(i);

    int NumEntries = 0;
    int ierr;
    ierr = ptr->ExtractMyRowCopy(i, MaxNumEntries,
				   NumEntries, &Values[0],
				   &Indices[0]);
    AMESOS_CHK_ERR(ierr);

    for (int j = 0 ; j < NumEntries ; ++j) {
      if (OnlyValues == false) {
	Row[count] = GlobalRow + 1;
	Col[count] = ptr->RowMatrixColMap().GID(Indices[j]) + 1;
      }
      
      // MS // Added on 15-Mar-05.
      if (AddToDiag_ && Indices[j] == i)
        Values[j] += AddToDiag_;

      Val[count] = Values[j];
      count++;
    }
  }

  MtxConvTime_ = AddTime("Total matrix conversion time", MtxConvTime_);
  
  assert (count <= ptr->NumMyNonzeros());

  return(0);
}