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