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

// ============================================================================
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");

  long long Rows = Matrix.NumGlobalRows64();
  long long Cols = Matrix.NumGlobalRows64();
  long long Nonzeros = Matrix.NumGlobalNonzeros64();

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

        long long GRID = Matrix.RowMatrixRowMap().GID64(i);

        for (int j = 0; j < NumMyEntries; ++j)
          of << GRID << " " << Matrix.RowMatrixColMap().GID64(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 RowMatrixToMatrixMarketFile( const char *filename, const Epetra_RowMatrix & A, 
				 const char * matrixName,
				 const char *matrixDescription, 
				 bool writeHeader) {
  long long M = A.NumGlobalRows64();
  long long N = A.NumGlobalCols64();
  long long nz = A.NumGlobalNonzeros64();

  FILE * handle = 0;

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

    handle = fopen(filename,"w");
    if (!handle) {EPETRA_CHK_ERR(-1);}
    MM_typecode matcode;
    mm_initialize_typecode(&matcode);
    mm_set_matrix(&matcode);
    mm_set_coordinate(&matcode);
    mm_set_real(&matcode);

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

  if (A.RowMatrixRowMap().Comm().MyPID()==0) // Only PE 0 opened a file
    if (fclose(handle)!=0) {EPETRA_CHK_ERR(-1);}
  return(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.NumGlobalCols64()==B.NumGlobalCols64(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalDiagonals64()==B.NumGlobalDiagonals64(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalNonzeros64()==B.NumGlobalNonzeros64(),ierr);
  EPETRA_TEST_ERR(!A.NumGlobalRows64()==B.NumGlobalRows64(),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 Ifpack_AnalyzeMatrixElements(const Epetra_RowMatrix& A,
                                 const bool abs, const int steps)
{

  bool verbose = (A.Comm().MyPID() == 0);
  double min_val =  DBL_MAX;
  double max_val = -DBL_MAX;

  std::vector<int>    colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

  for (int i = 0 ; i < A.NumMyRows() ; ++i) {

    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    for (int j = 0 ; j < Nnz ; ++j) {
      double v = colVal[j];
      if (abs)
        if (v < 0) v = -v;
      if (v < min_val)
        min_val = v;
      if (v > max_val)
        max_val = v;
    }
  }

  if (verbose) {
    cout << endl;
    Ifpack_PrintLine();
    cout << "Label of matrix = " << A.Label() << endl;
    cout << endl;
  }

  double delta = (max_val - min_val) / steps;
  for (int k = 0 ; k < steps ; ++k) {

    double below = delta * k + min_val;
    double above = below + delta;
    int MyBelow = 0, GlobalBelow;

    for (int i = 0 ; i < A.NumMyRows() ; ++i) {

      int Nnz;
      IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                        &colVal[0],&colInd[0]));

      for (int j = 0 ; j < Nnz ; ++j) {
        double v = colVal[j];
        if (abs)
          if (v < 0) v = -v;
        if (v >= below && v < above) MyBelow++;
      }

    }
    A.Comm().SumAll(&MyBelow, &GlobalBelow, 1);
    if (verbose) {
      printf("Elements in [%+7e, %+7e) = %10d ( = %5.2f %%)\n",
             below, above, GlobalBelow,
             100.0 * GlobalBelow / A.NumGlobalNonzeros64());
    }
  }

  if (verbose) {
    Ifpack_PrintLine();
    cout << endl;
  }

  return(0);
}

int Ifpack_Analyze(const Epetra_RowMatrix& A, const bool Cheap,
                   const int NumPDEEqns)
{

  int NumMyRows = A.NumMyRows();
  long long NumGlobalRows = A.NumGlobalRows64();
  long long NumGlobalCols = A.NumGlobalCols64();
  long long MyBandwidth = 0, GlobalBandwidth;
  long long MyLowerNonzeros = 0, MyUpperNonzeros = 0;
  long long GlobalLowerNonzeros, GlobalUpperNonzeros;
  long long MyDiagonallyDominant = 0, GlobalDiagonallyDominant;
  long long MyWeaklyDiagonallyDominant = 0, GlobalWeaklyDiagonallyDominant;
  double MyMin, MyAvg, MyMax;
  double GlobalMin, GlobalAvg, GlobalMax;
  long long GlobalStorage;

  bool verbose = (A.Comm().MyPID() == 0);

  GlobalStorage = sizeof(int*) * NumGlobalRows +
    sizeof(int) * A.NumGlobalNonzeros64() +
    sizeof(double) * A.NumGlobalNonzeros64();

  if (verbose) {
    print();
    Ifpack_PrintLine();
    print<const char*>("Label", A.Label());
    print<long long>("Global rows", NumGlobalRows);
    print<long long>("Global columns", NumGlobalCols);
    print<long long>("Stored nonzeros", A.NumGlobalNonzeros64());
    print<long long>("Nonzeros / row", A.NumGlobalNonzeros64() / NumGlobalRows);
    print<double>("Estimated storage (Mbytes)", 1.0e-6 * GlobalStorage);
  }

  long long NumMyActualNonzeros = 0, NumGlobalActualNonzeros;
  long long NumMyEmptyRows = 0, NumGlobalEmptyRows;
  long long NumMyDirichletRows = 0, NumGlobalDirichletRows;

  std::vector<int> colInd(A.MaxNumEntries());
  std::vector<double> colVal(A.MaxNumEntries());

  Epetra_Vector Diag(A.RowMatrixRowMap());
  Epetra_Vector RowSum(A.RowMatrixRowMap());
  Diag.PutScalar(0.0);
  RowSum.PutScalar(0.0);

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

    long long GRID = A.RowMatrixRowMap().GID64(i);
    int Nnz;
    IFPACK_CHK_ERR(A.ExtractMyRowCopy(i,A.MaxNumEntries(),Nnz,
                                      &colVal[0],&colInd[0]));

    if (Nnz == 0)
      NumMyEmptyRows++;

    if (Nnz == 1)
      NumMyDirichletRows++;

    for (int j = 0 ; j < Nnz ; ++j) {

      double v = colVal[j];
      if (v < 0) v = -v;
      if (colVal[j] != 0.0)
        NumMyActualNonzeros++;

      long long GCID = A.RowMatrixColMap().GID64(colInd[j]);

      if (GCID != GRID)
        RowSum[i] += v;
      else
        Diag[i] = v;

      if (GCID < GRID)
        MyLowerNonzeros++;
      else if (GCID > GRID)
        MyUpperNonzeros++;
      long long b = GCID - GRID;
      if (b < 0) b = -b;
      if (b > MyBandwidth)
        MyBandwidth = b;
    }

    if (Diag[i] > RowSum[i])
      MyDiagonallyDominant++;

    if (Diag[i] >= RowSum[i])
      MyWeaklyDiagonallyDominant++;

    RowSum[i] += Diag[i];
  }

  // ======================== //
  // summing up global values //
  // ======================== //

  A.Comm().SumAll(&MyDiagonallyDominant,&GlobalDiagonallyDominant,1);
  A.Comm().SumAll(&MyWeaklyDiagonallyDominant,&GlobalWeaklyDiagonallyDominant,1);
  A.Comm().SumAll(&NumMyActualNonzeros, &NumGlobalActualNonzeros, 1);
  A.Comm().SumAll(&NumMyEmptyRows, &NumGlobalEmptyRows, 1);
  A.Comm().SumAll(&NumMyDirichletRows, &NumGlobalDirichletRows, 1);
//.........这里部分代码省略.........