C++ Epetra_Vector::SetFlopCounter方法代码示例

int
powerMethod (double & lambda, 
             Epetra_CrsMatrix& A, 
             const int niters, 
             const double tolerance,
             const bool verbose)
{
  // In the power iteration, z = A*q.  Thus, q must be in the domain
  // of A, and z must be in the range of A.  The residual vector is of
  // course in the range of A.
  Epetra_Vector q (A.OperatorDomainMap ());
  Epetra_Vector z (A.OperatorRangeMap ());
  Epetra_Vector resid (A.OperatorRangeMap ());

  Epetra_Flops* counter = A.GetFlopCounter();
  if (counter != 0) {
    q.SetFlopCounter(A);
    z.SetFlopCounter(A);
    resid.SetFlopCounter(A);
  }

  // Initialize the starting vector z with random data.
  z.Random();

  double normz, residual;
  int ierr = 1;
  for (int iter = 0; iter < niters; ++iter)
    {
      z.Norm2 (&normz);        // normz  := ||z||_2
      q.Scale (1.0/normz, z);  // q      := z / normz
      A.Multiply(false, q, z); // z      := A * q
      q.Dot(z, &lambda);       // lambda := dot (q, z)

      // Compute the residual vector and display status output every
      // 100 iterations, or if we have reached the maximum number of
      // iterations.
      if (iter % 100 == 0 || iter + 1 == niters)
        {
          resid.Update (1.0, z, -lambda, q, 0.0); // resid := A*q - lambda*q
          resid.Norm2 (&residual);                // residual := ||resid||_2
          if (verbose) 
            cout << "Iter = " << iter << "  Lambda = " << lambda 
                 << "  Residual of A*q - lambda*q = " << residual << endl;
        } 
      if (residual < tolerance) { // We've converged!
        ierr = 0;
        break;
      }
    }
  return ierr;
}

//.........这里部分代码省略.........

  double elapsed_time, total_flops, MFLOPs;
  Epetra_Time timer(Comm);

  int LevelFill = 0;
  if (argc > 2)  LevelFill = atoi(argv[2]);
  if (verbose) cout << "Using Level Fill = " << LevelFill << endl;
  int Overlap = 0;
  if (argc > 3) Overlap = atoi(argv[3]);
  if (verbose) cout << "Using Level Overlap = " << Overlap << endl;
  double Athresh = 0.0;
  if (argc > 4) Athresh = atof(argv[4]);
  if (verbose) cout << "Using Absolute Threshold Value of = " << Athresh << endl;

  double Rthresh = 1.0;
  if (argc > 5) Rthresh = atof(argv[5]);
  if (verbose) cout << "Using Relative Threshold Value of = " << Rthresh << endl;

  Ifpack_IlukGraph * IlukGraph = 0;
  Ifpack_CrsRiluk * ILUK = 0;

  if (LevelFill>-1) {
    elapsed_time = timer.ElapsedTime();
    IlukGraph = new Ifpack_IlukGraph(Ap->Graph(), LevelFill, Overlap);
    assert(IlukGraph->ConstructFilledGraph()==0);
    elapsed_time = timer.ElapsedTime() - elapsed_time;
    if (verbose) cout << "Time to construct ILUK graph = " << elapsed_time << endl;


    Epetra_Flops fact_counter;
  
    elapsed_time = timer.ElapsedTime();
    ILUK = new Ifpack_CrsRiluk(*IlukGraph);
    ILUK->SetFlopCounter(fact_counter);
    ILUK->SetAbsoluteThreshold(Athresh);
    ILUK->SetRelativeThreshold(Rthresh);
    //assert(ILUK->InitValues()==0);
    int initerr = ILUK->InitValues(*Ap);
    if (initerr!=0) {
      cout << endl << Comm << endl << "  InitValues error = " << initerr;
      if (initerr==1) cout << "  Zero diagonal found, warning error only";
      cout << endl << endl;
    }
    assert(ILUK->Factor()==0);
    elapsed_time = timer.ElapsedTime() - elapsed_time;
    total_flops = ILUK->Flops();
    MFLOPs = total_flops/elapsed_time/1000000.0;
    if (verbose) cout << "Time to compute preconditioner values = " 
		    << elapsed_time << endl
		    << "MFLOPS for Factorization = " << MFLOPs << endl;
    //cout << *ILUK << endl;
  double Condest;
  ILUK->Condest(false, Condest);

  if (verbose) cout << "Condition number estimate for this preconditioner = " << Condest << endl;
  }
  int Maxiter = 100;
  double Tolerance = 1.0E-8;

  Epetra_Flops counter;
  Ap->SetFlopCounter(counter);
  xp->SetFlopCounter(*Ap);
  bp->SetFlopCounter(*Ap);
  if (ILUK!=0) ILUK->SetFlopCounter(*Ap);

  elapsed_time = timer.ElapsedTime();