Tsqr_ParTest.hpp

// @HEADER
// ***********************************************************************
//
//                 Anasazi: Block Eigensolvers Package
//                 Copyright (2010) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

#ifndef __TSQR_Test_DistTest_hpp
#define __TSQR_Test_DistTest_hpp

#include <Tsqr_Config.hpp>
#include <Tsqr_Random_NormalGenerator.hpp>
#include <Tsqr_verifyTimerConcept.hpp>

#include <Tsqr_generateStack.hpp>
#include <Tsqr_DistTsqr.hpp>
#include <Tsqr_GlobalVerify.hpp>
#include <Tsqr_printGlobalMatrix.hpp>

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <vector>


namespace TSQR {
  namespace Test {

    template< class Ordinal, class Scalar >
    class DistTsqrVerifier {
      TSQR::Random::NormalGenerator< Ordinal, Scalar > gen_;
      Teuchos::RCP< MessengerBase< Ordinal > > const ordinalComm_;
      Teuchos::RCP< MessengerBase< Scalar > > const scalarComm_;
      std::string scalarTypeName_;
      std::ostream& out_;
      std::ostream& err_;
      const bool testFactorExplicit_, testFactorImplicit_;
      const bool humanReadable_, printMatrices_, debug_;

    public:
      typedef Ordinal ordinal_type;
      typedef Scalar scalar_type;
      typedef typename ScalarTraits< scalar_type >::magnitude_type magnitude_type;
      typedef typename std::pair< magnitude_type, magnitude_type > result_type;
      typedef Matrix< ordinal_type, scalar_type > matrix_type;
      
      DistTsqrVerifier (const Teuchos::RCP< MessengerBase< Ordinal > >& ordinalComm,
      const Teuchos::RCP< MessengerBase< Scalar > >& scalarComm,
      const std::vector<int>& seed,
      const std::string& scalarTypeName,
      std::ostream& out,
      std::ostream& err,
      const bool testFactorExplicit,
      const bool testFactorImplicit,
      const bool humanReadable,
      const bool printMatrices,
      const bool debug) :
  gen_ (seed), 
  ordinalComm_ (ordinalComm),
  scalarComm_ (scalarComm),
  scalarTypeName_ (scalarTypeName), 
  out_ (out), 
  err_ (err),
  testFactorExplicit_ (testFactorExplicit),
  testFactorImplicit_ (testFactorImplicit),
  humanReadable_ (humanReadable), 
  printMatrices_ (printMatrices),
  debug_ (debug)
      {}

      DistTsqrVerifier (const Teuchos::RCP< MessengerBase< Ordinal > >& ordinalComm,
      const Teuchos::RCP< MessengerBase< Scalar > >& scalarComm,
      const std::string& scalarTypeName,
      std::ostream& out,
      std::ostream& err,
      const bool testFactorExplicit,
      const bool testFactorImplicit,
      const bool humanReadable,
      const bool printMatrices,
      const bool debug) :
  ordinalComm_ (ordinalComm),
  scalarComm_ (scalarComm),
  scalarTypeName_ (scalarTypeName), 
  out_ (out), 
  err_ (err),
  testFactorExplicit_ (testFactorExplicit),
  testFactorImplicit_ (testFactorImplicit),
  humanReadable_ (humanReadable), 
  printMatrices_ (printMatrices),
  debug_ (debug)
      {}

      void 
      getSeed (std::vector<int>& seed) const
      {
  gen_.getSeed (seed);
      }

      void 
      verify (const Ordinal numCols)
      {
  using std::endl;

  const int myRank = scalarComm_->rank();
  if (debug_)
    {
      scalarComm_->barrier();
      if (myRank == 0)
        err_ << "Verifying DistTsqr:" << endl;
      scalarComm_->barrier();
    }

  // Generate test problem.
  Matrix< Ordinal, Scalar > A_local, Q_local, R;
  testProblem (A_local, Q_local, R, numCols);
  if (debug_)
    {
      scalarComm_->barrier();
      if (myRank == 0)
        err_ << "-- Generated test problem." << endl;
      scalarComm_->barrier();
    }

  // Set up TSQR implementation.
  DistTsqr< Ordinal, Scalar > par (scalarComm_);
  if (debug_)
    {
      scalarComm_->barrier();
      if (myRank == 0)
        err_ << "-- DistTsqr object initialized" << endl << endl;
    }

  // Whether we've printed column headers yet (only matters for
  // non-humanReadable output).
  bool printedHeaders = false;

  // Test DistTsqr::factor() and DistTsqr::explicit_Q().
  if (testFactorImplicit_)
    {
      // Factor the matrix A (copied into R, which will be
      // overwritten on output)
      typedef typename DistTsqr< Ordinal, Scalar >::FactorOutput 
        factor_output_type;
      factor_output_type factorOutput = par.factor (R.view());
      if (debug_)
        {
    scalarComm_->barrier();
    if (myRank == 0)
      err_ << "-- Finished DistTsqr::factor" << endl;
        }
      // Compute the explicit Q factor
      par.explicit_Q (numCols, Q_local.get(), Q_local.lda(), factorOutput);
      if (debug_)
        {
    scalarComm_->barrier();
    if (myRank == 0)
      err_ << "-- Finished DistTsqr::explicit_Q" << endl;
        }
      // Verify the factorization
      result_type result = 
        global_verify (numCols, numCols, A_local.get(), A_local.lda(),
           Q_local.get(), Q_local.lda(), R.get(), R.lda(), 
           scalarComm_.get());
      if (debug_)
        {
    scalarComm_->barrier();
    if (myRank == 0)
      err_ << "-- Finished global_verify" << endl;
        }
      reportResults ("DistTsqr", numCols, result, (! printedHeaders));
      if (! printedHeaders)
        printedHeaders = true;
    }

  // Test DistTsqr::factorExplicit()
  if (testFactorExplicit_)
    {
      // Factor the matrix and compute the explicit Q factor, both
      // in a single operation.
      par.factorExplicit (R.view(), Q_local.view());
      if (debug_)
        {
    scalarComm_->barrier();
    if (myRank == 0)
      err_ << "-- Finished DistTsqr::factorExplicit" << endl;
        }

      if (printMatrices_)
        {
    if (myRank == 0)
      err_ << std::endl << "Computed Q factor:" << std::endl;
    printGlobalMatrix (err_, Q_local, scalarComm_.get(), ordinalComm_.get());
    if (myRank == 0)
      {
        err_ << std::endl << "Computed R factor:" << std::endl;
        print_local_matrix (err_, R.nrows(), R.ncols(), R.get(), R.lda());
        err_ << std::endl;
      }
        }

      // Verify the factorization
      result_type result = 
        global_verify (numCols, numCols, A_local.get(), A_local.lda(),
           Q_local.get(), Q_local.lda(), R.get(), R.lda(), 
           scalarComm_.get());
      if (debug_)
        {
    scalarComm_->barrier();
    if (myRank == 0)
      err_ << "-- Finished global_verify" << endl;
        }
      reportResults ("DistTsqrRB", numCols, result, (! printedHeaders));
      if (! printedHeaders)
        printedHeaders = true;
    }
      }

    private:
      void 
      reportResults (const std::string& method, 
         const Ordinal numCols, 
         const result_type& result,
         const bool printHeaders)
      {
  using std::endl;

  const int numProcs = scalarComm_->size();
  const int myRank = scalarComm_->rank();

  if (myRank == 0)
    {
      if (humanReadable_)
        {
    out_ << method << " accuracy results:" << endl
         << "Scalar type = " << scalarTypeName_ << endl
         << "Number of columns = " << numCols << endl
         << "Number of (MPI) processes = " << numProcs << endl
         << "Relative residual $\\|A - Q*R\\|_2 / \\|A\\|_2$ = " 
         << result.first << endl
         << "Relative orthogonality $\\|I - Q^T*Q\\|_2$ = " 
         << result.second << endl;
        }
      else
        {
    // Use scientific notation for floating-point numbers
    out_ << std::scientific;

    if (printHeaders)
      out_ << "%method,scalarType,numCols,numProcs,relResid,relOrthog" << endl;

    out_ << method
         << "," << scalarTypeName_
         << "," << numCols
         << "," << numProcs
         << "," << result.first
         << "," << result.second
         << endl;
        }
    }
      }

      void 
      testProblem (Matrix< Ordinal, Scalar >& A_local,
       Matrix< Ordinal, Scalar >& Q_local,
       Matrix< Ordinal, Scalar >& R,
       const Ordinal numCols)
      {
  const Ordinal numRowsLocal = numCols;

  // A_local: Space for the matrix A to factor -- local to each
  //   processor.
  //
  // A_global: Global matrix (only nonempty on Proc 0); only
  //   used temporarily.
  Matrix< Ordinal, Scalar > A_global;

  // This modifies A_local on all procs, and A_global on Proc 0.
  par_tsqr_test_problem (gen_, A_local, A_global, numCols, scalarComm_);
  
  if (printMatrices_)
    {
      const int myRank = scalarComm_->rank();

      if (myRank == 0)
        err_ << "Input matrix A:" << std::endl;
      printGlobalMatrix (err_, A_local, scalarComm_.get(), ordinalComm_.get());
      if (myRank == 0)
        err_ << std::endl;
    }

  // Copy the test problem input into R, since the factorization
  // will overwrite it in place with the final R factor.
  R.reshape (numCols, numCols);
  R.fill (Scalar(0));
  R.copy (A_local);

  // Prepare space in which to construct the explicit Q factor
  // (local component on this processor)
  Q_local.reshape (numRowsLocal, numCols);
  Q_local.fill (Scalar(0));
      }
    };


    template< class Ordinal, class Scalar, class TimerType >
    class DistTsqrBenchmarker {
      TSQR::Random::NormalGenerator< Ordinal, Scalar > gen_;
      Teuchos::RCP< MessengerBase< Scalar > > scalarComm_;
      Teuchos::RCP< MessengerBase< double > > doubleComm_; 
      std::string scalarTypeName_;

      std::ostream& out_;
      std::ostream& err_;
      const bool testFactorExplicit_, testFactorImplicit_;
      const bool humanReadable_, debug_;

    public:
      typedef Ordinal ordinal_type;
      typedef Scalar scalar_type;
      typedef typename ScalarTraits< scalar_type >::magnitude_type magnitude_type;
      typedef TimerType timer_type;

      DistTsqrBenchmarker (const Teuchos::RCP< MessengerBase< Scalar > >& scalarComm,
         const Teuchos::RCP< MessengerBase< double > >& doubleComm,
         const std::vector<int>& seed,
         const std::string& scalarTypeName,
         std::ostream& out,
         std::ostream& err,
         const bool testFactorExplicit,
         const bool testFactorImplicit,
         const bool humanReadable,
         const bool debug) :
  gen_ (seed), 
  scalarComm_ (scalarComm),
  doubleComm_ (doubleComm),
  scalarTypeName_ (scalarTypeName), 
  out_ (out), 
  err_ (err),
  testFactorExplicit_ (testFactorExplicit),
  testFactorImplicit_ (testFactorImplicit),
  humanReadable_ (humanReadable), 
  debug_ (debug)
      {}

      DistTsqrBenchmarker (const Teuchos::RCP< MessengerBase< Scalar > >& scalarComm,
         const Teuchos::RCP< MessengerBase< double > >& doubleComm,
         const std::string& scalarTypeName,
         std::ostream& out,
         std::ostream& err,
         const bool testFactorExplicit,
         const bool testFactorImplicit,
         const bool humanReadable,
         const bool debug) :
  scalarComm_ (scalarComm),
  doubleComm_ (doubleComm),
  scalarTypeName_ (scalarTypeName), 
  out_ (out), 
  err_ (err),
  testFactorExplicit_ (testFactorExplicit),
  testFactorImplicit_ (testFactorImplicit),
  humanReadable_ (humanReadable), 
  debug_ (debug)
      {}

      void 
      getSeed (std::vector<int>& seed) const
      {
  gen_.getSeed (seed);
      }

      void 
      benchmark (const int numTrials, const Ordinal numCols)
      {
  using std::endl;

  // Set up test problem.
  Matrix< Ordinal, Scalar > A_local, Q_local, R;
  testProblem (A_local, Q_local, R, numCols);

  // Set up TSQR implementation.
  DistTsqr< Ordinal, Scalar > par (scalarComm_);

  // Whether we've printed column headers yet (only matters for
  // non-humanReadable output).
  bool printedHeaders = false;

  if (testFactorImplicit_)
    {
      std::string timerName ("DistTsqr");

      // Benchmark DistTsqr (factor() and explicit_Q()) for
      // numTrials trials.
      timer_type timer (timerName);
      timer.start();
      for (int trialNum = 0; trialNum < numTrials; ++trialNum)
        {
    // Factor the matrix A (copied into R, which will be
    // overwritten on output)
    typedef typename DistTsqr< Ordinal, Scalar >::FactorOutput 
      factor_output_type;
    factor_output_type factorOutput = par.factor (R.view());

    // Compute the explicit Q factor
    par.explicit_Q (numCols, Q_local.get(), Q_local.lda(), factorOutput);
        }
      // Cumulative timing on this MPI process.
      // "Cumulative" means the elapsed time of numTrials executions.
      const double localCumulativeTiming = timer.stop();

      // reportResults() must be called on all processes, since this
      // figures out the min and max timings over all processes.
      reportResults (timerName, numTrials, numCols, localCumulativeTiming, 
         (! printedHeaders));
      if (! printedHeaders)
        printedHeaders = true;
    }

  if (testFactorExplicit_)
    {
      std::string timerName ("DistTsqrRB");

      // Benchmark DistTsqr::factorExplicit() for numTrials trials.
      timer_type timer (timerName);
      timer.start();
      for (int trialNum = 0; trialNum < numTrials; ++trialNum)
        {
    par.factorExplicit (R.view(), Q_local.view());
        }
      // Cumulative timing on this MPI process.
      // "Cumulative" means the elapsed time of numTrials executions.
      const double localCumulativeTiming = timer.stop();

      // Report cumulative (not per-invocation) timing results
      reportResults (timerName, numTrials, numCols, 
         localCumulativeTiming, (! printedHeaders));
      if (! printedHeaders)
        printedHeaders = true;

      // Per-invocation timings (for factorExplicit() benchmark
      // only).  localTimings were computed on this MPI process;
      // globalTimings are statistical summaries of those over
      // all MPI processes.  We only collect that data for
      // factorExplicit().
      std::vector< TimeStats > localTimings;
      std::vector< TimeStats > globalTimings;
      par.getFactorExplicitTimings (localTimings);
      for (std::vector< TimeStats >::size_type k = 0; k < localTimings.size(); ++k)
        globalTimings.push_back (TimeStats::globalTimeStats (doubleComm_.get(), localTimings[k]));
      std::vector< std::string > timingLabels;
      par.getFactorExplicitTimingLabels (timingLabels);

      if (humanReadable_)
        out_ << timerName << " per-invocation benchmark results:" << endl;

      const std::string labelLabel ("label,scalarType");
      for (std::vector< std::string >::size_type k = 0; k < timingLabels.size(); ++k)
        {
    // Only print column headers once
    const bool printHeaders = (k == 0);
    globalTimings[k].print (out_, humanReadable_, 
          timingLabels[k] + "," + scalarTypeName_, 
          labelLabel, printHeaders);
        }
    }
      }

    private:
      void 
      reportResults (const std::string& method,
         const int numTrials,
         const ordinal_type numCols,
         const double localTiming,
         const bool printHeaders)
      {
  using std::endl;

  // Find min and max timing over all MPI processes
  TimeStats localStats;
  localStats.update (localTiming);
  TimeStats globalStats = TimeStats::globalTimeStats (doubleComm_.get(), localStats);

  // Only Rank 0 prints the final results.
  const bool printResults = (doubleComm_->rank() == 0);
  if (printResults)
    {
      const int numProcs = doubleComm_->size();
      if (humanReadable_)
        {
    out_ << method << " cumulative benchmark results (total time over all trials):" << endl
         << "Scalar type = " << scalarTypeName_ << endl
         << "Number of columns = " << numCols << endl
         << "Number of (MPI) processes = " << numProcs << endl
         << "Number of trials = " << numTrials << endl
         << "Min timing (in seconds) = " << globalStats.min() << endl
         << "Mean timing (in seconds) = " << globalStats.mean() << endl
         << "Max timing (in seconds) = " << globalStats.max() << endl
         << endl;
        }
      else
        {
    // Use scientific notation for floating-point numbers
    out_ << std::scientific;

    if (printHeaders)
      out_ << "%method,scalarType,numCols,numProcs,numTrials,min,mean,max" << endl;

    out_ << method
         << "," << scalarTypeName_
         << "," << numCols
         << "," << numProcs
         << "," << numTrials
         << "," << globalStats.min()
         << "," << globalStats.mean()
         << "," << globalStats.max()
         << endl;
        }
    }
      }

      void 
      testProblem (Matrix< Ordinal, Scalar >& A_local,
       Matrix< Ordinal, Scalar >& Q_local,
       Matrix< Ordinal, Scalar >& R,
       const Ordinal numCols)
      {
  const Ordinal numRowsLocal = numCols;

  // A_local: Space for the matrix A to factor -- local to each
  //   processor.
  //
  // A_global: Global matrix (only nonempty on Proc 0); only
  //   used temporarily.
  Matrix< Ordinal, Scalar > A_global;

  // This modifies A_local on all procs, and A_global on Proc 0.
  par_tsqr_test_problem (gen_, A_local, A_global, numCols, scalarComm_);

  // Copy the test problem input into R, since the factorization
  // will overwrite it in place with the final R factor.
  R.reshape (numCols, numCols);
  R.copy (A_local);

  // Prepare space in which to construct the explicit Q factor
  // (local component on this processor)
  Q_local.reshape (numRowsLocal, numCols);
  Q_local.fill (Scalar(0));
      }

      static void
      conceptChecks () 
      {
  verifyTimerConcept< timer_type >();
      }
    };


  } // namespace Test
} // namespace TSQR

#endif // __TSQR_Test_DistTest_hpp