Anasazi Version of the Day
Tsqr_CombineBenchmarker.hpp
00001 #ifndef __Tsqr_CombineBenchmarker_hpp
00002 #define __Tsqr_CombineBenchmarker_hpp
00003 
00004 // @HEADER
00005 // ***********************************************************************
00006 //
00007 //                 Anasazi: Block Eigensolvers Package
00008 //                 Copyright (2010) Sandia Corporation
00009 //
00010 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00011 // license for use of this work by or on behalf of the U.S. Government.
00012 //
00013 // This library is free software; you can redistribute it and/or modify
00014 // it under the terms of the GNU Lesser General Public License as
00015 // published by the Free Software Foundation; either version 2.1 of the
00016 // License, or (at your option) any later version.
00017 //
00018 // This library is distributed in the hope that it will be useful, but
00019 // WITHOUT ANY WARRANTY; without even the implied warranty of
00020 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00021 // Lesser General Public License for more details.
00022 //
00023 // You should have received a copy of the GNU Lesser General Public
00024 // License along with this library; if not, write to the Free Software
00025 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00026 // USA
00027 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00028 //
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 #include "Tsqr_CombineBenchmark.hpp"
00033 
00034 #include <Tsqr_Random_NormalGenerator.hpp>
00035 #include <Tsqr_Random_MatrixGenerator.hpp>
00036 #include <Tsqr_verifyTimerConcept.hpp>
00037 
00038 #include <Tsqr_ApplyType.hpp>
00039 #include <Tsqr_Matrix.hpp>
00040 #include <Tsqr_ScalarTraits.hpp>
00041 #include <Tsqr_Util.hpp>
00042 
00043 #include <algorithm>
00044 #include <iostream>
00045 #include <limits>
00046 #include <sstream>
00047 #include <stdexcept>
00048 #include <vector>
00049 
00052 
00053 namespace TSQR {
00054   namespace Test {
00055 
00056     template< class Ordinal, class Scalar, class CombineType, class TimerType >
00057     class CombineBenchmarker {
00058     public:
00059       typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00060       typedef TSQR::Random::NormalGenerator< Ordinal, Scalar > normgen_type;
00061       typedef TSQR::Random::NormalGenerator< Ordinal, magnitude_type > normgen_mag_type;
00062       typedef TSQR::Random::MatrixGenerator< Ordinal, Scalar, normgen_type > matgen_type;
00063       typedef Matrix< Ordinal, Scalar > matrix_type;
00064       typedef std::pair< magnitude_type, magnitude_type > results_type;
00065 
00066     private:
00067       normgen_type gen_;
00068       normgen_mag_type magGen_;
00069       int numTrials_;
00070       bool debug_;
00071 
00072     public:
00073 
00074       CombineBenchmarker (normgen_type& gen, 
00075         normgen_mag_type& magGen, 
00076         const int numTrials,
00077         const bool debug = false) : 
00078   gen_ (gen), magGen_ (magGen), numTrials_ (numTrials), debug_ (debug)
00079       {
00080   TSQR::Test::verifyTimerConcept< TimerType >();
00081       }
00082 
00083       double 
00084       R1R2_benchmark (const Ordinal numCols)
00085       {
00086   using std::cerr;
00087   using std::endl;
00088   using std::invalid_argument;
00089   using std::ostringstream;
00090   using std::vector;
00091 
00092   if (numCols == 0)
00093     throw invalid_argument ("ncols == 0 is not allowed");
00094 
00095   // Generate two different sets of singular values.
00096   // Randomly perturb them, but make sure all are positive.
00097 
00098   vector< magnitude_type > sigma_R1 (numCols);
00099   vector< magnitude_type > sigma_R2 (numCols);
00100   generateSingularValues (sigma_R1, numCols);
00101   generateSingularValues (sigma_R2, numCols);
00102 
00103   matrix_type R1 (numCols, numCols, Scalar(0));
00104   matrix_type R2 (numCols, numCols, Scalar(0));
00105 
00106   matgen_type matgen (gen_);
00107   matgen.fill_random_R (numCols, R1.get(), R1.lda(), &sigma_R1[0]);
00108   matgen.fill_random_R (numCols, R2.get(), R2.lda(), &sigma_R2[0]);
00109 
00110   // Space to put the original test problem, expressed as one
00111   // dense matrix rather than in two blocks.  This will be a
00112   // deep copy of the test problem, since the test problem
00113   // matrices will be overwritten by the factorizations.
00114   matrix_type A_R1R2 (Ordinal(2) * numCols, numCols, Scalar(0));
00115 
00116   // Copy [R1; R2] into A_R1R2.
00117   copy_matrix (numCols, numCols, &A_R1R2(0, 0), A_R1R2.lda(), 
00118          R1.get(), R1.lda());
00119   copy_matrix (numCols, numCols, &A_R1R2(numCols, 0), A_R1R2.lda(), 
00120          R2.get(), R2.lda());
00121 
00122   // Space to put the explicit Q factor
00123   matrix_type Q_R1R2 (Ordinal(2) * numCols, numCols, Scalar(0));
00124 
00125   // Fill the explicit Q factor matrix with the first numCols
00126   // columns of the identity matrix.
00127   for (Ordinal k = 0; k < numCols; ++k)
00128     Q_R1R2(k, k) = Scalar(1);
00129 
00130   // tau factor array
00131   vector< Scalar > tau_R1R2 (numCols);
00132 
00133   // Workspace array for factorization and applying the Q factor.
00134   vector< Scalar > work (numCols);
00135 
00136   if (debug_)
00137     cerr << endl << "----------------------------------------" << endl
00138          << "TSQR::Combine test problem:" << endl
00139          << "qr( [R1; R2] ), with R1 and R2 " << numCols 
00140          << " by " << numCols << endl << endl;
00141 
00142   CombineType combiner; 
00143   TimerType timer ("Combine");
00144   timer.start ();
00145   for (int trialNum = 0; trialNum < numTrials_; ++trialNum)
00146     {
00147       combiner.factor_pair (numCols, R1.get(), R1.lda(), R2.get(), R2.lda(),
00148           &tau_R1R2[0], &work[0]);
00149       combiner.apply_pair (ApplyType("N"), numCols, numCols, 
00150          R2.get(), R2.lda(), &tau_R1R2[0], 
00151          &Q_R1R2(0, 0), Q_R1R2.lda(),
00152          &Q_R1R2(numCols, 0), Q_R1R2.lda(),
00153          &work[0]);
00154     }
00155   const double timingResult = timer.stop();
00156   if (debug_)
00157     cerr << "Results: " << numTrials_ << " consecutive runs took a total of " 
00158          << timingResult << " seconds" << endl;
00159   return timingResult;
00160       }
00161 
00162       double 
00163       RA_benchmark (const Ordinal numRows, const Ordinal numCols)
00164       {
00165   using std::cerr;
00166   using std::endl;
00167   using std::invalid_argument;
00168   using std::ostringstream;
00169   using std::vector;
00170 
00171   if (numRows < numCols)
00172     {
00173       ostringstream os;
00174       os << "# rows < # columns is not allowed.  You specified # rows = " 
00175          << numRows << " and # columns = " << numCols << ".";
00176       throw invalid_argument (os.str());
00177     }
00178   else if (numCols == 0)
00179     throw invalid_argument ("ncols == 0 is not allowed");
00180 
00181   // Generate two different sets of singular values.  Randomly
00182   // perturb them, but make sure all are positive.
00183   vector< magnitude_type > sigma_R3 (numCols);
00184   vector< magnitude_type > sigma_A (numCols);
00185   generateSingularValues (sigma_R3, numCols);
00186   generateSingularValues (sigma_A, numCols);
00187 
00188   // Generate the test problem [R3; A]
00189   matrix_type R3 (numCols, numCols, Scalar(0));
00190   matrix_type A (numRows, numCols, Scalar(0));
00191   matgen_type matgen (gen_);
00192   matgen.fill_random_R (numCols, R3.get(), R3.lda(), &sigma_R3[0]);
00193   matgen.fill_random_svd (numRows, numCols, A.get(), A.lda(), &sigma_A[0]);
00194 
00195   // Space to put the original test problem, expressed as one
00196   // dense matrix rather than in two blocks.  This will be a deep
00197   // copy of the test problem, since the test problem matrices
00198   // will be overwritten by the factorization.
00199   matrix_type A_R3A (numRows + numCols, numCols, Scalar(0));
00200 
00201   // Copy [R3; A] into A_R3A.
00202   copy_matrix (numCols, numCols, &A_R3A(0, 0), A_R3A.lda(), 
00203          R3.get(), R3.lda());
00204   copy_matrix (numRows, numCols, &A_R3A(numCols, 0), A_R3A.lda(), 
00205          A.get(), A.lda());
00206 
00207   // Space to put the explicit Q factor
00208   matrix_type Q_R3A (numRows + numCols, numCols, Scalar(0));
00209 
00210   // Fill the explicit Q factor matrix with the first numCols
00211   // columns of the identity matrix.
00212   for (Ordinal k = 0; k < numCols; ++k)
00213     Q_R3A(k, k) = Scalar(1);
00214 
00215   // tau factor array
00216   vector< Scalar > tau_R3A (numCols);
00217 
00218   // Workspace array for factorization and applying the Q factor.
00219   vector< Scalar > work (numCols);
00220 
00221   if (debug_)
00222     cerr << endl << "----------------------------------------" << endl
00223          << "TSQR::Combine test problem:" << endl
00224          << "qr( [R3; A] ), with R3 " << numCols << " by " << numCols 
00225          << " and A " << numRows << " by " << numCols << endl << endl;
00226 
00227   CombineType combiner; 
00228   TimerType timer ("Combine");
00229   timer.start ();
00230   for (int trialNum = 0; trialNum < numTrials_; ++trialNum)
00231     {
00232       combiner.factor_inner (numRows, numCols, R3.get(), R3.lda(),
00233            A.get(), A.lda(), &tau_R3A[0], &work[0]);
00234       combiner.apply_inner (ApplyType("N"), numRows, numCols, numCols,
00235           A.get(), A.lda(), &tau_R3A[0], 
00236           &Q_R3A(0, 0), Q_R3A.lda(),
00237           &Q_R3A(numCols, 0), Q_R3A.lda(), 
00238           &work[0]);
00239     }
00240   const double timingResult = timer.stop();
00241   if (debug_)
00242     cerr << "Results: " << numTrials_ << " consecutive runs took a total of " 
00243          << timingResult << " seconds" << endl;
00244   return timingResult;
00245       }
00246 
00247     private:
00248       void
00249       generateSingularValues (std::vector< magnitude_type >& sigma,
00250             const Ordinal numValues)
00251       {
00252   const magnitude_type machEps = 
00253     std::numeric_limits< magnitude_type >::epsilon();
00254   sigma.resize (numValues);
00255     
00256   // Relative amount by which to perturb each singular value.  The
00257   // perturbation will be multiplied by a normal(0,1) pseudorandom
00258   // number drawn from magGen.
00259   const magnitude_type perturbationFactor = magnitude_type(10) * machEps;
00260 
00261   sigma[0] = magnitude_type (1);
00262   for (Ordinal k = 1; k < numValues; ++k)
00263     {
00264       const magnitude_type perturbation = perturbationFactor * magGen_();
00265       const magnitude_type beforePerturb = sigma[k-1] / magnitude_type(2);
00266       const magnitude_type candidate = beforePerturb + perturbation;
00267 
00268       // If adding the perturbation to beforePerturb would result
00269       // in a nonpositive number, subtract instead.
00270       if (candidate <= magnitude_type(0))
00271         sigma[k] = beforePerturb - perturbation;
00272       else
00273         sigma[k] = candidate;
00274     }
00275       }
00276     };
00277 
00278   } // namespace Test
00279 } // namespace TSQR
00280 
00281 #endif // __Tsqr_CombineBenchmarker_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends