Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Tsqr_CombineTest.cpp
00001 //@HEADER
00002 // ************************************************************************
00003 // 
00004 //          Kokkos: Node API and Parallel Node Kernels
00005 //              Copyright (2008) Sandia Corporation
00006 // 
00007 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 // 
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00038 // 
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #include "Tsqr_ConfigDefs.hpp"
00043 #include "Tsqr_CombineTest.hpp"
00044 
00045 #include <Teuchos_Assert.hpp>
00046 #include <Tsqr_Random_NormalGenerator.hpp>
00047 #include <Tsqr_Random_MatrixGenerator.hpp>
00048 
00049 #include <Tsqr_Combine.hpp>
00050 #include <Tsqr_LocalVerify.hpp>
00051 #include <Tsqr_Matrix.hpp>
00052 #include <Tsqr_ScalarTraits.hpp>
00053 #include <Tsqr_Util.hpp>
00054 
00055 #include <algorithm>
00056 #include <iostream>
00057 #include <limits>
00058 #include <sstream>
00059 #include <stdexcept>
00060 #include <vector>
00061 
00062 
00063 namespace TSQR {
00064   namespace Test {
00065 
00066     template<class Ordinal, class MagnitudeType, class NormalGenType>
00067     static void 
00068     generateSingularValues (NormalGenType& magGen,
00069           std::vector<MagnitudeType>& sigma,
00070           const Ordinal numValues)
00071     {
00072       typedef MagnitudeType magnitude_type;
00073       const magnitude_type machEps = 
00074   std::numeric_limits<magnitude_type>::epsilon();
00075       sigma.resize (numValues);
00076     
00077       // Relative amount by which to perturb each singular value.  The
00078       // perturbation will be multiplied by a normal(0,1) pseudorandom
00079       // number drawn from magGen.
00080       const magnitude_type perturbationFactor = magnitude_type(10) * machEps;
00081 
00082       sigma[0] = magnitude_type (1);
00083       for (Ordinal k = 1; k < numValues; ++k)
00084   {
00085     const magnitude_type perturbation = perturbationFactor * magGen();
00086     const magnitude_type beforePerturb = sigma[k-1] / magnitude_type(2);
00087     const magnitude_type candidate = beforePerturb + perturbation;
00088 
00089     // If adding the perturbation to beforePerturb would result
00090     // in a nonpositive number, subtract instead.
00091     if (candidate <= magnitude_type(0))
00092       sigma[k] = beforePerturb - perturbation;
00093     else
00094       sigma[k] = candidate;
00095   }
00096     }
00097 
00098     static void
00099     printCombineFieldNames ()
00100     {
00101       using std::cout;
00102       using std::endl;
00103 
00104       const char prefix[] = "%";
00105       cout << prefix
00106      << "method" 
00107      << ",kernel"
00108      << ",scalarType"
00109      << ",numRows"
00110      << ",numCols"
00111      << ",absFrobResid"
00112      << ",absFrobOrthog"
00113      << ",frobA"
00114      << endl;
00115     }
00116 
00117     template<class MagnitudeType>
00118     static void
00119     printR1R2results (const std::string& datatype,
00120           const int numCols,
00121           const std::vector<MagnitudeType>& results)
00122     {
00123       using std::cout;
00124       using std::endl;
00125 
00126       cout << "Combine"
00127      << "," << "R1R2"
00128      << "," << datatype
00129      << "," << (2*numCols)
00130      << "," << numCols
00131      << "," << results[0]
00132      << "," << results[1]
00133      << "," << results[2]
00134      << endl;
00135     }
00136 
00137     template<class MagnitudeType>
00138     static void
00139     printR3Aresults (const std::string& datatype,
00140          const int numRows,
00141          const int numCols,
00142          const std::vector<MagnitudeType>& results)
00143     {
00144       using std::cout;
00145       using std::endl;
00146 
00147       cout << "Combine"
00148      << "," << "R3A"
00149      << "," << datatype
00150      << "," << numRows
00151      << "," << numCols
00152      << "," << results[3]
00153      << "," << results[4]
00154      << "," << results[5]
00155      << endl;
00156     }
00157 
00158     template<class MagnitudeType>
00159     static void
00160     printResults (const std::string& datatype,
00161       const int numRows,
00162       const int numCols,
00163       const std::vector<MagnitudeType>& results,
00164       const bool printFieldNames)
00165     {
00166       if (printFieldNames)
00167   printCombineFieldNames();
00168       printR1R2results (datatype, numCols, results);
00169       printR3Aresults (datatype, numRows, numCols, results);
00170     }
00171 
00172     template<class MagnitudeType>
00173     static void
00174     printSimSeqTsqrResults (const std::string& datatype, 
00175           const int numRows, 
00176           const int numCols, 
00177           const std::vector<MagnitudeType>& results,
00178           const bool printFieldNames)
00179     {
00180       using std::cout;
00181       using std::endl;
00182 
00183       if (printFieldNames)
00184   {
00185     const char prefix[] = "%";
00186     cout << prefix
00187          << "method"
00188          << ",scalarType"
00189          << ",numRows"
00190          << ",numCols"
00191          << ",absFrobResid"
00192          << ",absFrobOrthog"
00193          << ",frobA"
00194          << endl;
00195   }
00196       cout << "CombineSimSeqTsqr"
00197      << "," << datatype
00198      << "," << numRows
00199      << "," << numCols
00200      << "," << results[0]
00201      << "," << results[1]
00202      << "," << results[2]
00203      << endl;
00204     }
00205 
00206     template< class MatrixViewType >
00207     static void
00208     printMatrix (std::ostream& out,
00209      const MatrixViewType& A)
00210     {
00211       print_local_matrix (out, A.nrows(), A.ncols(), A.get(), A.lda());
00212     }
00213 
00214     template< class MatrixViewType >
00215     static
00216     std::vector< typename ScalarTraits< typename MatrixViewType::scalar_type >::magnitude_type >
00217     localVerify (const MatrixViewType& A, 
00218      const MatrixViewType& Q, 
00219      const MatrixViewType& R)
00220     {
00221       return local_verify (A.nrows(), A.ncols(), A.get(), A.lda(), 
00222          Q.get(), Q.lda(), R.get(), R.lda());
00223     }
00224     
00236     template< class Ordinal, class Scalar >
00237     static std::vector< typename ScalarTraits< Scalar >::magnitude_type >
00238     verifyCombineTemplate (TSQR::Random::NormalGenerator< Ordinal, Scalar >& gen,
00239          TSQR::Random::NormalGenerator< Ordinal, typename ScalarTraits< Scalar >::magnitude_type >& magGen,
00240          const Ordinal numRows, 
00241          const Ordinal numCols,
00242          const bool debug)
00243     {
00244       using TSQR::Random::MatrixGenerator;
00245       using TSQR::Random::NormalGenerator;
00246       using std::cerr;
00247       using std::endl;
00248       using std::invalid_argument;
00249       using std::ostringstream;
00250       using std::pair;
00251       using std::vector;
00252 
00253       typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00254       typedef NormalGenerator< Ordinal, Scalar > normgen_type;
00255       typedef NormalGenerator< Ordinal, magnitude_type > normgen_mag_type;
00256       typedef MatrixGenerator< Ordinal, Scalar, normgen_type > matgen_type;
00257       typedef Matrix< Ordinal, Scalar > matrix_type;
00258       typedef vector< magnitude_type > results_type;
00259 
00260       if (numRows < numCols)
00261   {
00262     ostringstream os;
00263     os << "# rows < # columns is not allowed.  You specified # rows = " 
00264        << numRows << " and # columns = " << numCols << ".";
00265     throw invalid_argument (os.str());
00266   }
00267       else if (numCols == 0)
00268   throw invalid_argument ("ncols == 0 is not allowed");
00269 
00270       //
00271       // Generate four different sets of singular values.  Randomly
00272       // perturb them, but make sure all are positive.
00273       //
00274       vector< magnitude_type > sigma_R1 (numCols);
00275       vector< magnitude_type > sigma_R2 (numCols);
00276       vector< magnitude_type > sigma_R3 (numCols);
00277       vector< magnitude_type > sigma_A (numCols);
00278       generateSingularValues (magGen, sigma_R1, numCols);
00279       generateSingularValues (magGen, sigma_R2, numCols);
00280       generateSingularValues (magGen, sigma_R3, numCols);
00281       generateSingularValues (magGen, sigma_A, numCols);
00282 
00283       matrix_type R1 (numCols, numCols, Scalar(0));
00284       matrix_type R2 (numCols, numCols, Scalar(0));
00285       matrix_type R3 (numCols, numCols, Scalar(0));
00286       matrix_type A (numRows, numCols, Scalar(0));
00287       matgen_type matgen (gen);
00288       matgen.fill_random_R (numCols, R1.get(), R1.lda(), &sigma_R1[0]);
00289       matgen.fill_random_R (numCols, R2.get(), R2.lda(), &sigma_R2[0]);
00290       matgen.fill_random_R (numCols, R3.get(), R3.lda(), &sigma_R3[0]);
00291       matgen.fill_random_svd (numRows, numCols, A.get(), A.lda(), &sigma_A[0]);
00292 
00293       if (false && debug)
00294   {
00295     cerr << endl << "First test problem:" << endl;
00296     print_local_matrix (cerr, numCols, numCols, R1.get(), R1.lda());
00297     print_local_matrix (cerr, numCols, numCols, R2.get(), R2.lda());
00298     cerr << endl;
00299 
00300     cerr << endl << "Second test problem:" << endl;
00301     print_local_matrix (cerr, numCols, numCols, R3.get(), R3.lda());
00302     print_local_matrix (cerr, numRows, numCols, A.get(), A.lda());
00303     cerr << endl;
00304   }
00305 
00306       // Space to put the original test problem, expressed as one
00307       // dense matrix rather than in two blocks.  These will be deep
00308       // copies of the test problems, since the test problem matrices
00309       // will be overwritten by the factorizations.
00310       matrix_type A_R1R2 (Ordinal(2) * numCols, numCols, Scalar(0));
00311       matrix_type A_R3A (numRows + numCols, numCols, Scalar(0));
00312 
00313       // Copy [R1; R2] into A_R1R2.
00314       copy_matrix (numCols, numCols, &A_R1R2(0, 0), A_R1R2.lda(), 
00315        R1.get(), R1.lda());
00316       copy_matrix (numCols, numCols, &A_R1R2(numCols, 0), A_R1R2.lda(), 
00317        R2.get(), R2.lda());
00318 
00319       // Copy [R3; A] into A_R3A.
00320       copy_matrix (numCols, numCols, &A_R3A(0, 0), A_R3A.lda(), 
00321        R3.get(), R3.lda());
00322       copy_matrix (numRows, numCols, &A_R3A(numCols, 0), A_R3A.lda(), 
00323        A.get(), A.lda());
00324 
00325       // Space to put the explicit Q factors.
00326       matrix_type Q_R1R2 (Ordinal(2) * numCols, numCols, Scalar(0));
00327       matrix_type Q_R3A (numRows + numCols, numCols, Scalar(0));
00328 
00329       // Fill the explicit Q factor matrices with the first numCols
00330       // columns of the identity matrix.
00331       for (Ordinal k = 0; k < numCols; ++k)
00332   {
00333     Q_R1R2(k, k) = Scalar(1);
00334     Q_R3A(k, k) = Scalar(1);
00335   }
00336 
00337       // tau factor arrays, one for each factorization test.
00338       vector< Scalar > tau_R1R2 (numCols);
00339       vector< Scalar > tau_R3A (numCols);
00340 
00341       // Workspace array for factorization and applying the Q factor.
00342       // We recycle this workspace for all tests.
00343       vector< Scalar > work (numCols);
00344 
00345       if (debug)
00346   cerr << endl << "----------------------------------------" << endl
00347        << "TSQR::Combine first test problem:" << endl
00348        << "qr( [R1; R2] ), with R1 and R2 " << numCols 
00349        << " by " << numCols << endl << endl;
00350 
00351       Combine< Ordinal, Scalar > combiner;      
00352       combiner.factor_pair (numCols, R1.get(), R1.lda(), R2.get(), R2.lda(),
00353           &tau_R1R2[0], &work[0]);
00354       combiner.apply_pair (ApplyType("N"), numCols, numCols, 
00355          R2.get(), R2.lda(), &tau_R1R2[0], 
00356          &Q_R1R2(0, 0), Q_R1R2.lda(),
00357          &Q_R1R2(numCols, 0), Q_R1R2.lda(),
00358          &work[0]);
00359       if (debug)
00360   {
00361     cerr << "Results of first test problem:" << endl;
00362     cerr << "-- Copy of test problem:" << endl;
00363     print_local_matrix (cerr, A_R1R2.nrows(), A_R1R2.ncols(), 
00364             A_R1R2.get(), A_R1R2.lda());
00365     cerr << endl << "-- Q factor:" << endl;
00366     print_local_matrix (cerr, Q_R1R2.nrows(), Q_R1R2.ncols(), 
00367             Q_R1R2.get(), Q_R1R2.lda());
00368     cerr << endl << "-- R factor:" << endl;
00369     print_local_matrix (cerr, R1.nrows(), R1.ncols(), 
00370             R1.get(), R1.lda());
00371     cerr << endl;
00372   }
00373       const results_type firstResults = 
00374   local_verify (A_R1R2.nrows(), A_R1R2.ncols(), 
00375           A_R1R2.get(), A_R1R2.lda(),
00376           Q_R1R2.get(), Q_R1R2.lda(),
00377           R1.get(), R1.lda());
00378       if (debug)
00379   cerr << "\\| A - Q*R \\|_F = " << firstResults[0] << endl
00380        << "\\| I - Q'*Q \\|_F = " << firstResults[1] << endl
00381        << "\\| A \\|_A = " << firstResults[2] << endl;
00382 
00383       if (debug)
00384   cerr << endl << "----------------------------------------" << endl
00385        << "TSQR::Combine second test problem:" << endl
00386        << "qr( [R3; A] ), with R3 " << numCols << " by " << numCols 
00387        << " and A " << numRows << " by " << numCols << endl << endl;
00388 
00389       combiner.factor_inner (numRows, numCols, R3.get(), R3.lda(),
00390            A.get(), A.lda(), &tau_R3A[0], &work[0]);
00391       combiner.apply_inner (ApplyType("N"), numRows, numCols, numCols,
00392           A.get(), A.lda(), &tau_R3A[0], 
00393           &Q_R3A(0, 0), Q_R3A.lda(),
00394           &Q_R3A(numCols, 0), Q_R3A.lda(), 
00395           &work[0]);
00396       if (debug)
00397   {
00398     cerr << "Results of second test problem:" << endl;
00399     cerr << "-- Copy of test problem:" << endl;
00400     print_local_matrix (cerr, A_R3A.nrows(), A_R3A.ncols(), 
00401             A_R3A.get(), A_R3A.lda());
00402     cerr << endl << "-- Q factor:" << endl;
00403     print_local_matrix (cerr, Q_R3A.nrows(), Q_R3A.ncols(), 
00404             Q_R3A.get(), Q_R3A.lda());
00405     cerr << endl << "-- R factor:" << endl;
00406     print_local_matrix (cerr, R3.nrows(), R3.ncols(), 
00407             R3.get(), R3.lda());
00408     cerr << endl;
00409   }
00410       const results_type secondResults = 
00411   local_verify (A_R3A.nrows(), A_R3A.ncols(), 
00412           A_R3A.get(), A_R3A.lda(),
00413           Q_R3A.get(), Q_R3A.lda(),
00414           R3.get(), R3.lda());
00415       if (debug)
00416   cerr << "\\| A - Q*R \\|_F = " << secondResults[0] << endl
00417        << "\\| I - Q'*Q \\|_F = " << secondResults[1] << endl
00418        << "\\| A \\|_A = " << secondResults[2] << endl;
00419 
00420       vector< magnitude_type > finalResults;
00421       finalResults.push_back (firstResults[0]);
00422       finalResults.push_back (firstResults[1]);
00423       finalResults.push_back (firstResults[2]);
00424 
00425       finalResults.push_back (secondResults[0]);
00426       finalResults.push_back (secondResults[1]);
00427       finalResults.push_back (secondResults[2]);
00428       return finalResults;
00429     }
00430 
00431 
00432 
00433 
00436     template< class Ordinal, class Scalar >
00437     static std::vector< typename ScalarTraits< Scalar >::magnitude_type >
00438     verifyCombineSeqTemplate (TSQR::Random::NormalGenerator< Ordinal, Scalar >& gen,
00439             TSQR::Random::NormalGenerator< Ordinal, typename ScalarTraits< Scalar >::magnitude_type >& magGen,
00440             const Ordinal numRows, 
00441             const Ordinal numCols,
00442             const bool debug)
00443     {
00444       using TSQR::Random::MatrixGenerator;
00445       using TSQR::Random::NormalGenerator;
00446       using std::cerr;
00447       using std::endl;
00448       using std::invalid_argument;
00449       using std::ostringstream;
00450       using std::pair;
00451       using std::vector;
00452 
00453       typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00454       typedef NormalGenerator< Ordinal, Scalar > normgen_type;
00455       typedef NormalGenerator< Ordinal, magnitude_type > normgen_mag_type;
00456       typedef MatrixGenerator< Ordinal, Scalar, normgen_type > matgen_type;
00457       typedef Matrix< Ordinal, Scalar > matrix_type;
00458       typedef MatView< Ordinal, Scalar > mat_view_type;
00459       typedef vector< magnitude_type > results_type;
00460 
00461       if (numRows < numCols)
00462   {
00463     ostringstream os;
00464     os << "# rows < # columns is not allowed.  You specified # rows = " 
00465        << numRows << " and # columns = " << numCols << ".";
00466     throw invalid_argument (os.str());
00467   }
00468       else if (numCols == 0)
00469   throw invalid_argument ("ncols == 0 is not allowed");
00470 
00471       // Generate two different sets of singular values. 
00472       vector< magnitude_type > sigma_A1 (numCols);
00473       vector< magnitude_type > sigma_A2 (numCols);
00474       generateSingularValues (magGen, sigma_A1, numCols);
00475       generateSingularValues (magGen, sigma_A2, numCols);
00476 
00477       // Matrix consisting of two cache blocks.
00478       matrix_type A (Ordinal(2)*numRows, numCols, Scalar(0));
00479       // Views of the two cache blocks.
00480       mat_view_type A1 (numRows, numCols, &A(0,0), A.lda());
00481       mat_view_type A2 (numRows, numCols, &A(numRows,0), A.lda());
00482 
00483       // Fill the two cache blocks with random test problems.
00484       matgen_type matgen (gen);
00485       matgen.fill_random_svd (numRows, numCols, A1.get(), A1.lda(), &sigma_A1[0]);
00486       matgen.fill_random_svd (numRows, numCols, A2.get(), A2.lda(), &sigma_A2[0]);
00487 
00488       if (false && debug)
00489   {
00490     cerr << endl << "Test problem:" << endl;
00491     cerr << endl << "Original matrix:" << endl;
00492     printMatrix (cerr, A);
00493     cerr << endl << "First cache block:" << endl;
00494     printMatrix (cerr, A1);
00495     cerr << endl << "Second cache block:" << endl;
00496     printMatrix (cerr, A2);
00497     cerr << endl;
00498   }
00499 
00500       // Copy of the resulting test problem, stored as one dense
00501       // matrix rather than as two blocks.  We will use A_copy to
00502       // measure the residual error once we've completed the
00503       // factorization and computed the explicit Q factor.
00504       matrix_type A_copy (A);
00505 
00506       // Space to put the explicit Q factor.
00507       matrix_type Q (Ordinal(2) * numRows, numCols, Scalar(0));
00508 
00509       // Fill Q with the first numCols columns of the identity matrix.
00510       for (Ordinal k = 0; k < numCols; ++k)
00511   Q(k, k) = Scalar(1);
00512 
00513       // Two cache blocks (as views) of Q.
00514       mat_view_type Q1 (numRows, numCols, &Q(0,0), Q.lda());
00515       mat_view_type Q2 (numRows, numCols, &Q(numRows,0), Q.lda());
00516 
00517       // Two tau factor arrays, one for each cache block.
00518       vector< Scalar > tau1 (numCols);
00519       vector< Scalar > tau2 (numCols);
00520 
00521       // Workspace array for factorization and applying the Q factor.
00522       // We recycle this workspace for all tests.
00523       vector< Scalar > work (numCols);
00524 
00525       if (debug)
00526   cerr << endl << "----------------------------------------" << endl
00527        << "TSQR::Combine SequentialTsqr simulation with 2 cache blocks:" 
00528        << endl << "qr( [A1; A2] ), with A1 and A2 being each "
00529        << numRows << " by " << numCols << endl << endl;
00530 
00531       Combine< Ordinal, Scalar > combiner;      
00532       // qr( A1 )
00533       combiner.factor_first (numRows, numCols, A1.get(), A1.lda(), 
00534            &tau1[0], &work[0]);
00535       // View of numCols by numCols upper triangle of A1.
00536       mat_view_type R1 (numCols, numCols, A1.get(), A1.lda());
00537       // qr( [R1; A2] )
00538       combiner.factor_inner (numRows, numCols, R1.get(), R1.lda(),
00539            A2.get(), A2.lda(), &tau2[0], &work[0]);
00540       // Extract (a deep copy of) the R factor.
00541       matrix_type R (R1);
00542       // Zero out everything below the diagonal of R.
00543       for (Ordinal j = 0; j < numCols; ++j)
00544   for (Ordinal i = j+1; i < numCols; ++i)
00545     R(i,j) = Scalar(0);
00546 
00547       // Compute the explicit Q factor, by starting with A2 and
00548       // (working up the matrix A,) finishing with A1.
00549       combiner.apply_inner (ApplyType::NoTranspose, 
00550           numRows, numCols, numCols,
00551           A2.get(), A2.lda(), &tau2[0], 
00552           Q1.get(), Q1.lda(),
00553           Q2.get(), Q2.lda(), &work[0]);
00554       combiner.apply_first (ApplyType::NoTranspose,
00555           numRows, numCols, numCols,
00556           A1.get(), A.lda(), &tau1[0],
00557           Q1.get(), Q1.lda(), &work[0]);
00558       if (debug)
00559   {
00560     cerr << "Results of first test problem:" << endl;
00561     cerr << "-- Test matrix A:" << endl;
00562     printMatrix (cerr, A_copy);
00563     cerr << endl << "-- Q factor:" << endl;
00564     printMatrix (cerr, Q);
00565     cerr << endl << "-- R factor:" << endl;
00566     printMatrix (cerr, R);
00567     cerr << endl;
00568   }
00569       const results_type results = localVerify (A_copy, Q, R);
00570 
00571       if (debug)
00572   cerr << "\\| A - Q*R \\|_F = " << results[0] << endl
00573        << "\\| I - Q'*Q \\|_F = " << results[1] << endl
00574        << "\\| A \\|_F = " << results[2] << endl;
00575 
00576       return results;
00577     }
00578 
00579 
00580     void
00581     verifyCombine (const int numRows,
00582        const int numCols, 
00583        const bool testReal,
00584        const bool testComplex,
00585        const bool printFieldNames,
00586        const bool simulateSequentialTsqr,
00587        const bool debug)
00588     {
00589       using TSQR::Random::NormalGenerator;
00590       using std::cerr;
00591 #ifdef HAVE_TSQR_COMPLEX
00592       using std::complex;
00593 #endif // HAVE_TSQR_COMPLEX
00594       using std::cout;
00595       using std::endl;
00596       using std::pair;
00597       using std::string;
00598       using std::vector;
00599 
00600       //
00601       // We do tests one after another, using the seed from the
00602       // previous test in the current test, so that the pseudorandom
00603       // streams used by the tests are independent.
00604       //
00605 
00606       // Default seed for the next pseudorandom number generator.
00607       // This will be the same each time, so if you want
00608       // nondeterministic behavior, you should modify this routine to
00609       // let you supply the seed values yourself.
00610       vector<int> iseed(4);
00611       iseed[0] = 0;
00612       iseed[1] = 0;
00613       iseed[2] = 0;
00614       iseed[3] = 1;
00615 
00616       // Whether to print the field (i.e., column) names for the
00617       // output data.
00618       bool doPrintFieldNames = printFieldNames;
00619 
00620       if (! simulateSequentialTsqr)
00621   {
00622     if (testReal)
00623       {
00624         {
00625     NormalGenerator<int, float> normgenS (iseed);
00626     const vector<float> resultsS = 
00627       verifyCombineTemplate (normgenS, normgenS, numRows, 
00628            numCols, debug);
00629     // Only print field names (if at all) once per run, for
00630     // the first data type.
00631     printResults (string("float"), numRows, numCols, 
00632             resultsS, doPrintFieldNames);
00633     // Print field names at most once.
00634     doPrintFieldNames = false; 
00635     // Fetch the pseudorandom seed from the previous test.
00636     normgenS.getSeed (iseed);
00637         }
00638         {
00639     NormalGenerator<int, double> normgenD (iseed);
00640     const vector<double> resultsD = 
00641       verifyCombineTemplate (normgenD, normgenD, numRows, 
00642            numCols, debug);
00643     printResults (string("double"), numRows, numCols, 
00644             resultsD, doPrintFieldNames);
00645     doPrintFieldNames = false; 
00646     normgenD.getSeed (iseed);
00647         }
00648       }
00649 
00650     if (testComplex)
00651       {
00652 #ifdef HAVE_TSQR_COMPLEX
00653         {
00654     NormalGenerator<int, complex<float> > normgenC (iseed);
00655     NormalGenerator<int, float> normgenS (iseed);
00656     const vector<float> resultsC = 
00657       verifyCombineTemplate (normgenC, normgenS, numRows, 
00658            numCols, debug);
00659     printResults (string("complex<float>"), numRows, numCols, 
00660             resultsC, doPrintFieldNames);
00661     doPrintFieldNames = false; 
00662     // Even though normgenC and normgenS each updated the
00663     // random seed independently, for now we just fetch the
00664     // updated seed from normgenC.  This should still
00665     // produce reproducible results.
00666     normgenC.getSeed (iseed);
00667         }
00668         {
00669     NormalGenerator<int, complex<double> > normgenZ (iseed);
00670     NormalGenerator<int, double> normgenD (iseed);
00671     const vector<double> resultsZ = 
00672       verifyCombineTemplate (normgenZ, normgenD, numRows, 
00673            numCols, debug);
00674     printResults (string("complex<double>"), numRows, numCols,
00675             resultsZ, doPrintFieldNames);
00676     doPrintFieldNames = false;
00677     normgenZ.getSeed (iseed);
00678         }
00679 #else // NOT HAVE_TSQR_COMPLEX
00680         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00681          "Trilinos was not built with "
00682          "complex arithmetic support");
00683 #endif // HAVE_TSQR_COMPLEX
00684       }
00685   }
00686       else // simulateSequentialTsqr
00687   {
00688     if (testReal)
00689       {
00690         {
00691     NormalGenerator<int, float> normgenS (iseed);
00692     const vector<float> resultsS = 
00693       verifyCombineSeqTemplate (normgenS, normgenS, numRows, 
00694               numCols, debug);
00695     printSimSeqTsqrResults (string("float"), numRows, numCols, 
00696           resultsS, doPrintFieldNames);
00697     doPrintFieldNames = false;
00698     normgenS.getSeed (iseed);
00699         }
00700         {
00701     NormalGenerator<int, double> normgenD (iseed);
00702     const vector<double> resultsD = 
00703       verifyCombineSeqTemplate (normgenD, normgenD, numRows, 
00704               numCols, debug);
00705     printSimSeqTsqrResults (string("double"), numRows, numCols, 
00706           resultsD, doPrintFieldNames);
00707     doPrintFieldNames = false;
00708     normgenD.getSeed (iseed);
00709         }
00710       }
00711 
00712     if (testComplex)
00713       {
00714 #ifdef HAVE_TSQR_COMPLEX
00715         {
00716     NormalGenerator<int, complex<float> > normgenC (iseed);
00717     NormalGenerator<int, float> normgenS (iseed);
00718     const vector<float> resultsC = 
00719       verifyCombineSeqTemplate (normgenC, normgenS, numRows, 
00720               numCols, debug);
00721     printSimSeqTsqrResults (string("complex<float>"), numRows, numCols, 
00722           resultsC, doPrintFieldNames);
00723     doPrintFieldNames = false;
00724     normgenC.getSeed (iseed);
00725         }
00726         {
00727     NormalGenerator<int, complex<double> > normgenZ (iseed);
00728     NormalGenerator<int, double> normgenD (iseed);
00729     const vector<double> resultsZ = 
00730       verifyCombineSeqTemplate (normgenZ, normgenD, numRows, 
00731               numCols, debug);
00732     printSimSeqTsqrResults (string("complex<double>"), numRows, 
00733           numCols, resultsZ, doPrintFieldNames);
00734     doPrintFieldNames = false;
00735     normgenZ.getSeed (iseed);
00736         }
00737 #else // NOT HAVE_TSQR_COMPLEX
00738         TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00739          "Trilinos was not built with "
00740          "complex arithmetic support");
00741 #endif // HAVE_TSQR_COMPLEX
00742       }
00743   }
00744     }
00745   } // namespace Test
00746 } // namespace TSQR
00747 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends