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

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_PrintSparsity(const Epetra_RowMatrix& A, const char* InputFileName,
                         const int NumPDEEqns)
{

  int ltit;
  long long m,nc,nr,maxdim;
  double lrmrgn,botmrgn,xtit,ytit,ytitof,fnstit,siz = 0.0;
  double xl,xr, yb,yt, scfct,u2dot,frlw,delt,paperx;
  bool square = false;
  /*change square to .true. if you prefer a square frame around
    a rectangular matrix */
  double conv = 2.54;
  char munt = 'E'; /* put 'E' for centimeters, 'U' for inches */
  int ptitle = 0; /* position of the title, 0 under the drawing,
                     else above */
  FILE* fp = NULL;
  int NumMyRows;
  //int NumMyCols;
  long long NumGlobalRows;
  long long NumGlobalCols;
  int MyPID;
  int NumProc;
  char FileName[1024];
  char title[1024];

  const Epetra_Comm& Comm = A.Comm();

  /* --------------------- execution begins ---------------------- */

  if (strlen(A.Label()) != 0)
    strcpy(title, A.Label());
  else
    sprintf(title, "%s", "matrix");

  if (InputFileName == 0)
    sprintf(FileName, "%s.ps", title);
  else
    strcpy(FileName, InputFileName);

  MyPID = Comm.MyPID();
  NumProc = Comm.NumProc();

  NumMyRows = A.NumMyRows();
  //NumMyCols = A.NumMyCols();

  NumGlobalRows = A.NumGlobalRows64();
  NumGlobalCols = A.NumGlobalCols64();

  if (NumGlobalRows != NumGlobalCols)
    IFPACK_CHK_ERR(-1); // never tested

  /* to be changed for rect matrices */
  maxdim = (NumGlobalRows>NumGlobalCols)?NumGlobalRows:NumGlobalCols;
  maxdim /= NumPDEEqns;

  m = 1 + maxdim;
  nr = NumGlobalRows / NumPDEEqns + 1;
  nc = NumGlobalCols / NumPDEEqns + 1;

  if (munt == 'E') {
    u2dot = 72.0/conv;
    paperx = 21.0;
    siz = 10.0;
  }
  else {
    u2dot = 72.0;
    paperx = 8.5*conv;
    siz = siz*conv;
  }

  /* left and right margins (drawing is centered) */

  lrmrgn = (paperx-siz)/2.0;

  /* bottom margin : 2 cm */

  botmrgn = 2.0;
  /* c scaling factor */
  scfct = siz*u2dot/m;
  /* matrix frame line witdh */
  frlw = 0.25;
  /* font size for title (cm) */
  fnstit = 0.5;
  /* mfh 23 Jan 2013: title is always nonnull, since it's an array of
     fixed nonzero length.  The 'if' test thus results in a compiler
     warning. */
  /*if (title) ltit = strlen(title);*/
  /*else       ltit = 0;*/
  ltit = strlen(title);

  /* position of title : centered horizontally */
  /*                     at 1.0 cm vertically over the drawing */
  ytitof = 1.0;
  xtit = paperx/2.0;
  ytit = botmrgn+siz*nr/m + ytitof;
  /* almost exact bounding box */
  xl = lrmrgn*u2dot - scfct*frlw/2;
  xr = (lrmrgn+siz)*u2dot + scfct*frlw/2;
  yb = botmrgn*u2dot - scfct*frlw/2;
//.........这里部分代码省略.........

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);
//.........这里部分代码省略.........