C++ Epetra_Flops::ResetFlops方法代码示例

void runJadMatrixTests(Epetra_JadMatrix * A,  Epetra_MultiVector * b, Epetra_MultiVector * bt,
		    Epetra_MultiVector * xexact, bool StaticProfile, bool verbose, bool summary) {

  Epetra_MultiVector z(*b);
  Epetra_MultiVector r(*b);
  Epetra_SerialDenseVector resvec(b->NumVectors());

  //Timings
  Epetra_Flops flopcounter;
  A->SetFlopCounter(flopcounter);
  Epetra_Time timer(A->Comm());

  for (int j=0; j<2; j++) { // j = 0 is notrans, j = 1 is trans

    bool TransA = (j==1);
    A->SetUseTranspose(TransA);
    flopcounter.ResetFlops();
    timer.ResetStartTime();

    //10 matvecs
    for( int i = 0; i < 10; ++i )
      A->Apply(*xexact, z); // Compute z = A*xexact or z = A'*xexact

    double elapsed_time = timer.ElapsedTime();
    double total_flops = A->Flops();

    // Compute residual
    if (TransA)
      r.Update(-1.0, z, 1.0, *bt, 0.0); // r = bt - z
    else
      r.Update(-1.0, z, 1.0, *b, 0.0); // r = b - z

    r.Norm2(resvec.Values());

    if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
    double MFLOPs = total_flops/elapsed_time/1000000.0;
    if (verbose) cout << "Total MFLOPs for 10 " << " Jagged Diagonal MatVec's with (Trans = " << TransA
		      << ") " << MFLOPs << " (" << elapsed_time << " s)" <<endl;
    if (summary) {
      if (A->Comm().NumProc()==1) {
	if (TransA) cout << "TransMv" << '\t';
	else cout << "NoTransMv" << '\t';
      }
      cout << MFLOPs << endl;
    }
  }
  return;
}

//.........这里部分代码省略.........
  //
  // Add rows to the sparse matrix one at a time.
  //
  std::vector<double> Values(2);
  Values[0] = -1.0; Values[1] = -1.0;
  std::vector<int> Indices(2);
  const double two = 2.0;
  int NumEntries;

  for (int i = 0; i < NumMyElements; ++i)
    {
      if (MyGlobalElements[i] == 0)
        { // The first row of the matrix.
          Indices[0] = 1;
          NumEntries = 1;
        }
      else if (MyGlobalElements[i] == NumGlobalElements - 1)
        { // The last row of the matrix.
          Indices[0] = NumGlobalElements-2;
          NumEntries = 1;
        }
      else
        { // Any row of the matrix other than the first or last.
          Indices[0] = MyGlobalElements[i]-1;
          Indices[1] = MyGlobalElements[i]+1;
          NumEntries = 2;
        }
      ierr = A.InsertGlobalValues(MyGlobalElements[i], NumEntries, &Values[0], &Indices[0]);
      assert (ierr==0);
      // Insert the diagonal entry.
      ierr = A.InsertGlobalValues(MyGlobalElements[i], 1, &two, &MyGlobalElements[i]);
      assert(ierr==0);
    }

  // Finish up.  We can call FillComplete() with no arguments, because
  // the matrix is square.
  ierr = A.FillComplete ();
  assert (ierr==0);

  // Parameters for the power method.
  const int niters = NumGlobalElements*10;
  const double tolerance = 1.0e-2;

  //
  // Run the power method.  Keep track of the flop count and the total
  // elapsed time.
  //
  Epetra_Flops counter;
  A.SetFlopCounter(counter);
  Epetra_Time timer(Comm);
  double lambda = 0.0;
  ierr += powerMethod (lambda, A, niters, tolerance, verbose);
  double elapsedTime = timer.ElapsedTime ();
  double totalFlops =counter.Flops ();
  // Mflop/s: Million floating-point arithmetic operations per second.
  double Mflop_per_s = totalFlops / elapsedTime / 1000000.0;

  if (verbose) 
    cout << endl << endl << "Total Mflop/s for first solve = " 
         << Mflop_per_s << endl<< endl;

  // Increase the first (0,0) diagonal entry of the matrix.
  if (verbose) 
    cout << endl << "Increasing magnitude of first diagonal term, solving again"
         << endl << endl << endl;

  if (A.MyGlobalRow (0)) {
    int numvals = A.NumGlobalEntries (0);
    std::vector<double> Rowvals (numvals);
    std::vector<int> Rowinds (numvals);
    A.ExtractGlobalRowCopy (0, numvals, numvals, &Rowvals[0], &Rowinds[0]); // Get A(0,0)
    for (int i = 0; i < numvals; ++i) 
      if (Rowinds[i] == 0) 
        Rowvals[i] *= 10.0;

    A.ReplaceGlobalValues (0, numvals, &Rowvals[0], &Rowinds[0]);
  }

  //
  // Run the power method again.  Keep track of the flop count and the
  // total elapsed time.
  //
  lambda = 0.0;
  timer.ResetStartTime();
  counter.ResetFlops();
  ierr += powerMethod (lambda, A, niters, tolerance, verbose);
  elapsedTime = timer.ElapsedTime();
  totalFlops = counter.Flops();
  Mflop_per_s = totalFlops / elapsedTime / 1000000.0;

  if (verbose) 
    cout << endl << endl << "Total Mflop/s for second solve = " 
         << Mflop_per_s << endl << endl;

#ifdef EPETRA_MPI
  MPI_Finalize() ;
#endif

  return ierr;
}

//=========================================================================================
void runLUMatrixTests(Epetra_CrsMatrix * L,  Epetra_MultiVector * bL, Epetra_MultiVector * btL, Epetra_MultiVector * xexactL,
		      Epetra_CrsMatrix * U,  Epetra_MultiVector * bU, Epetra_MultiVector * btU, Epetra_MultiVector * xexactU,
		      bool StaticProfile, bool verbose, bool summary) {

  if (L->NoDiagonal()) {
    bL->Update(1.0, *xexactL, 1.0); // Add contribution of a unit diagonal to bL
    btL->Update(1.0, *xexactL, 1.0); // Add contribution of a unit diagonal to btL
  }
  if (U->NoDiagonal()) {
    bU->Update(1.0, *xexactU, 1.0); // Add contribution of a unit diagonal to bU
    btU->Update(1.0, *xexactU, 1.0); // Add contribution of a unit diagonal to btU
  }

  Epetra_MultiVector z(*bL);
  Epetra_MultiVector r(*bL);
  Epetra_SerialDenseVector resvec(bL->NumVectors());

  //Timings
  Epetra_Flops flopcounter;
  L->SetFlopCounter(flopcounter);
  U->SetFlopCounter(flopcounter);
  Epetra_Time timer(L->Comm());
  std::string statdyn =        "dynamic";
  if (StaticProfile) statdyn = "static ";

  for (int j=0; j<4; j++) { // j = 0/2 is notrans, j = 1/3 is trans

    bool TransA = (j==1 || j==3);
    std::string contig = "without";
    if (j>1) contig =    "with   ";

    if (j==2) {
      L->OptimizeStorage();
      U->OptimizeStorage();
    }

    flopcounter.ResetFlops();
    timer.ResetStartTime();

    //10 lower solves
    bool Upper = false;
    bool UnitDiagonal = L->NoDiagonal();  // If no diagonal, then unit must be used
    Epetra_MultiVector * b = TransA ? btL : bL;  // solve with the appropriate b vector
    for( int i = 0; i < 10; ++i )
      L->Solve(Upper, TransA, UnitDiagonal, *b, z); // Solve Lz = bL or L'z = bLt

    double elapsed_time = timer.ElapsedTime();
    double total_flops = L->Flops();

    // Compute residual
    r.Update(-1.0, z, 1.0, *xexactL, 0.0); // r = bt - z
    r.Norm2(resvec.Values());

    if (resvec.NormInf()>0.000001) {
      cout << "resvec = " << resvec << endl;
      cout << "z = " << z << endl;
      cout << "xexactL = " << *xexactL << endl;
      cout << "r = " << r << endl;
    }

    if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
    double MFLOPs = total_flops/elapsed_time/1000000.0;
    if (verbose) cout << "Total MFLOPs for 10 " << " Lower solves " << statdyn << " Profile (Trans = " << TransA
		      << ")  and " << contig << " opt storage = " << MFLOPs << " (" << elapsed_time << " s)" <<endl;
    if (summary) {
      if (L->Comm().NumProc()==1) {
	if (TransA) cout << "TransLSv" << statdyn<< "Prof" << contig << "OptStor" << '\t';
	else cout << "NoTransLSv" << statdyn << "Prof" << contig << "OptStor" << '\t';
      }
      cout << MFLOPs << endl;
    }
    flopcounter.ResetFlops();
    timer.ResetStartTime();

    //10 upper solves
    Upper = true;
    UnitDiagonal = U->NoDiagonal();  // If no diagonal, then unit must be used
    b = TransA ? btU : bU;  // solve with the appropriate b vector
    for( int i = 0; i < 10; ++i )
      U->Solve(Upper, TransA, UnitDiagonal, *b, z); // Solve Lz = bL or L'z = bLt

    elapsed_time = timer.ElapsedTime();
    total_flops = U->Flops();

    // Compute residual
    r.Update(-1.0, z, 1.0, *xexactU, 0.0); // r = bt - z
    r.Norm2(resvec.Values());

    if (resvec.NormInf()>0.001) {
      cout << "U = " << *U << endl;
      //cout << "resvec = " << resvec << endl;
      cout << "z = " << z << endl;
      cout << "xexactU = " << *xexactU << endl;
      //cout << "r = " << r << endl;
      cout << "b = " << *b << endl;
    }


    if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
//.........这里部分代码省略.........

void runMatrixTests(Epetra_CrsMatrix * A,  Epetra_MultiVector * b, Epetra_MultiVector * bt,
		    Epetra_MultiVector * xexact, bool StaticProfile, bool verbose, bool summary) {

  Epetra_MultiVector z(*b);
  Epetra_MultiVector r(*b);
  Epetra_SerialDenseVector resvec(b->NumVectors());

  //Timings
  Epetra_Flops flopcounter;
  A->SetFlopCounter(flopcounter);
  Epetra_Time timer(A->Comm());
  std::string statdyn =        "dynamic";
  if (StaticProfile) statdyn = "static ";

  for (int j=0; j<4; j++) { // j = 0/2 is notrans, j = 1/3 is trans

    bool TransA = (j==1 || j==3);
    std::string contig = "without";
    if (j>1) contig =    "with   ";

#ifdef EPETRA_SHORT_PERFTEST
    int kstart = 1;
#else
    int kstart = 0;
#endif
    for (int k=kstart; k<2; k++) { // Loop over old multiply vs. new multiply

      std::string oldnew = "old";
      if (k>0) oldnew =    "new";

      if (j==2) A->OptimizeStorage();

      flopcounter.ResetFlops();
      timer.ResetStartTime();

      if (k==0) {
	//10 matvecs
#ifndef EPETRA_SHORT_PERFTEST
	for( int i = 0; i < 10; ++i )
	  A->Multiply1(TransA, *xexact, z); // Compute z = A*xexact or z = A'*xexact using old Multiply method
#endif
      }
      else {
	//10 matvecs
	for( int i = 0; i < 10; ++i )
	  A->Multiply(TransA, *xexact, z); // Compute z = A*xexact or z = A'*xexact
      }

      double elapsed_time = timer.ElapsedTime();
      double total_flops = A->Flops();

      // Compute residual
      if (TransA)
	r.Update(-1.0, z, 1.0, *bt, 0.0); // r = bt - z
      else
	r.Update(-1.0, z, 1.0, *b, 0.0); // r = b - z

      r.Norm2(resvec.Values());

      if (verbose) cout << "ResNorm = " << resvec.NormInf() << ": ";
      double MFLOPs = total_flops/elapsed_time/1000000.0;
      if (verbose) cout << "Total MFLOPs for 10 " << oldnew << " MatVec's with " << statdyn << " Profile (Trans = " << TransA
			<< ")  and " << contig << " optimized storage = " << MFLOPs << " (" << elapsed_time << " s)" <<endl;
      if (summary) {
	if (A->Comm().NumProc()==1) {
	  if (TransA) cout << "TransMv" << statdyn<< "Prof" << contig << "OptStor" << '\t';
	  else cout << "NoTransMv" << statdyn << "Prof" << contig << "OptStor" << '\t';
	}
	cout << MFLOPs << endl;
      }
    }
  }
  return;
}

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


      GenerateCrsProblem(numNodesX, numNodesY, numProcsX, numProcsY, XUoff.Length(),
			 XUoff.Values(), YUoff.Values(), nrhs, comm, verbose, summary,
			 mapU, U, bU, btU, xexactU, StaticProfile, true);


      runLUMatrixTests(L, bL, btL, xexactL, U, bU, btU, xexactU, StaticProfile, verbose, summary);

      delete L;
      delete bL;
      delete btL;
      delete xexactL;
      delete mapL;

      delete U;
      delete bU;
      delete btU;
      delete xexactU;
      delete mapU;

      Epetra_MultiVector q(*map, nrhs);
      Epetra_MultiVector z(q);
      Epetra_MultiVector r(q);

      delete map;
      q.SetFlopCounter(flopcounter);
      z.SetFlopCounter(q);
      r.SetFlopCounter(q);

      resvec.Resize(nrhs);


      flopcounter.ResetFlops();
      timer.ResetStartTime();

      //10 norms
      for( int i = 0; i < 10; ++i )
	q.Norm2( resvec.Values() );

      elapsed_time = timer.ElapsedTime();
      total_flops = q.Flops();
      MFLOPs = total_flops/elapsed_time/1000000.0;
      if (verbose) cout << "\nTotal MFLOPs for 10 Norm2's= " << MFLOPs << endl;

      if (summary) {
	if (comm.NumProc()==1) cout << "Norm2" << '\t';
	cout << MFLOPs << endl;
      }

      flopcounter.ResetFlops();
      timer.ResetStartTime();

      //10 dot's
      for( int i = 0; i < 10; ++i )
	q.Dot(z, resvec.Values());

      elapsed_time = timer.ElapsedTime();
      total_flops = q.Flops();
      MFLOPs = total_flops/elapsed_time/1000000.0;
      if (verbose) cout << "Total MFLOPs for 10 Dot's  = " << MFLOPs << endl;

      if (summary) {
	if (comm.NumProc()==1) cout << "DotProd" << '\t';
	cout << MFLOPs << endl;
      }

//.........这里部分代码省略.........
  // variable needed for iteration
  double lambda = 0.0;
  // int niters = 10000;
  int niters = 200;
  double tolerance = 1.0e-1;

  /////////////////////////////////////////////////////////////////////////////////////////////////
	
  // Iterate

  Epetra_Flops flopcounter;
  A.SetFlopCounter(flopcounter);
  q.SetFlopCounter(A);
  z.SetFlopCounter(A);
  resid.SetFlopCounter(A);
	

  Epetra_Time timer(Comm);
  EPETRA_TEST_ERR(power_method(false, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
  double elapsed_time = timer.ElapsedTime();
  double total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
  double MFLOPs = total_flops/elapsed_time/1000000.0;

  if (verbose) cout << "\n\nTotal MFLOPs for first solve = " << MFLOPs << std::endl<< std::endl;

  /////////////////////////////////////////////////////////////////////////////////////////////////
	
  // Solve transpose problem

  if (verbose) cout << "\n\nUsing transpose of matrix and solving again (should give same result).\n\n"
		    << std::endl;
  // Iterate
  lambda = 0.0;
  flopcounter.ResetFlops();
  timer.ResetStartTime();
  EPETRA_TEST_ERR(power_method(true, A, q, z, resid, &lambda, niters, tolerance, verbose),ierr);
  elapsed_time = timer.ElapsedTime();
  total_flops = A.Flops() + q.Flops() + z.Flops() + resid.Flops();
  MFLOPs = total_flops/elapsed_time/1000000.0;

  if (verbose) cout << "\n\nTotal MFLOPs for transpose solve = " << MFLOPs << std::endl<< endl;

  /////////////////////////////////////////////////////////////////////////////////////////////////

  // Increase diagonal dominance

  if (verbose) cout << "\n\nIncreasing the magnitude of first diagonal term and solving again\n\n"
		    << endl;

  
  if (A.MyGlobalRow(0)) {
    int numvals = A.NumGlobalEntries(0);
    double * Rowvals = new double [numvals];
    int    * Rowinds = new int    [numvals];
    A.ExtractGlobalRowCopy(0, numvals, numvals, Rowvals, Rowinds); // Get A[0,0]

    for (int i=0; i<numvals; i++) if (Rowinds[i] == 0) Rowvals[i] *= 10.0;
    
    A.ReplaceGlobalValues(0, numvals, Rowvals, Rowinds);
    delete [] Rowvals;
    delete [] Rowinds;
  }
  // Iterate (again)
  lambda = 0.0;
  flopcounter.ResetFlops();
  timer.ResetStartTime();