Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Tsqr_Random_MatrixGenerator.hpp
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 #ifndef __TSQR_Random_MatrixGenerator_hpp
00030 #define __TSQR_Random_MatrixGenerator_hpp
00031 
00032 #include <Tsqr_Blas.hpp>
00033 #include <Tsqr_Lapack.hpp>
00034 #include <Tsqr_Matrix.hpp>
00035 #include <Tsqr_ScalarTraits.hpp>
00036 
00037 #include <algorithm>
00038 #include <limits>
00039 #include <sstream>
00040 #include <stdexcept>
00041 #include <vector>
00042 
00045 
00046 namespace TSQR {
00047   namespace Random {
00048 
00049     template< class Ordinal, class Scalar, class Generator >
00050     class MatrixGenerator {
00051     public:
00052       typedef Ordinal ordinal_type;
00053       typedef Scalar scalar_type;
00054       typedef typename ScalarTraits< Scalar >::magnitude_type magnitude_type;
00055       typedef Generator generator_type;
00056 
00057       MatrixGenerator (Generator& generator) : gen_ (generator) {}
00058 
00059       void
00060       fill_random (const Ordinal nrows, 
00061        const Ordinal ncols, 
00062        Scalar A[], 
00063        const Ordinal lda)
00064       {
00065   for (Ordinal j = 0; j < ncols; ++j)
00066     {
00067       Scalar* const A_j = &A[j*lda];
00068       for (Ordinal i = 0; i < nrows; ++i)
00069         A_j[i] = gen_();
00070     }
00071       }
00072 
00080       void
00081       explicit_Q (const Ordinal nrows, 
00082       const Ordinal ncols, 
00083       Scalar Q[], 
00084       const Ordinal ldq)
00085       {
00086   // Fill Q with random numbers
00087   this->fill_random (nrows, ncols, Q, ldq);
00088 
00089   // Get ready for QR factorization
00090   LAPACK< Ordinal, Scalar > lapack;
00091   std::vector< Scalar > tau (std::min(nrows, ncols));
00092 
00093   // Workspace query
00094   Scalar _lwork1, _lwork2;
00095   int info = 0;
00096   lapack.GEQRF (nrows, ncols, Q, ldq, &tau[0], &_lwork1, -1, &info);
00097   if (info != 0)
00098     throw std::logic_error("LAPACK GEQRF LWORK query failed");
00099 
00100   lapack.ORGQR (nrows, ncols, ncols, Q, ldq, &tau[0], &_lwork2, -1, &info);
00101   if (info != 0)
00102     throw std::logic_error("LAPACK ORGQR LWORK query failed");
00103 
00104   // Allocate workspace.  abs() returns a magnitude_type, and we
00105   // can compare those using std::max.  If Scalar is complex,
00106   // you can't compare it using max.
00107   const Ordinal lwork = checkedCast (std::max (ScalarTraits< Scalar >::abs (_lwork1), 
00108                  ScalarTraits< Scalar >::abs (_lwork2)));
00109   std::vector< Scalar > work (lwork);
00110 
00111   // Factor the input matrix
00112   lapack.GEQRF (nrows, ncols, Q, ldq, &tau[0], &work[0], lwork, &info);
00113   if (info != 0)
00114     throw std::runtime_error("LAPACK GEQRF failed");
00115 
00116   // Compute explicit Q factor in place
00117   lapack.ORGQR (nrows, ncols, ncols, Q, ldq, &tau[0], &work[0], lwork, &info);
00118   if (info != 0)
00119     throw std::runtime_error("LAPACK ORGQR failed");
00120       }
00121 
00122 
00131       void
00132       implicit_Q (const Ordinal nrows, 
00133       const Ordinal ncols, 
00134       Scalar Q[], 
00135       const Ordinal ldq,
00136       Scalar tau[])
00137       {
00138   // Fill Q with random numbers
00139   this->fill_random (nrows, ncols, Q, ldq);
00140 
00141   // Get ready for QR factorization
00142   LAPACK< Ordinal, Scalar > lapack;
00143 
00144   // Workspace query
00145   Scalar _lwork1;
00146   int info = 0;
00147   lapack.GEQRF (nrows, ncols, Q, ldq, tau, &_lwork1, -1, &info);
00148   if (info != 0)
00149     throw std::logic_error("LAPACK GEQRF LWORK query failed");
00150 
00151   // Allocate workspace.
00152   const Ordinal lwork = checkedCast (ScalarTraits< Scalar >::abs (_lwork1));
00153   std::vector< Scalar > work (lwork);
00154 
00155   // Factor the input matrix
00156   lapack.GEQRF (nrows, ncols, Q, ldq, tau, &work[0], lwork, &info);
00157   if (info != 0)
00158     throw std::runtime_error("LAPACK GEQRF failed");
00159       }
00160 
00161       template< class MatrixViewType >
00162       void
00163       implicit_Q (MatrixViewType& Q, 
00164       typename MatrixViewType::scalar_type tau[])
00165       {
00166   implicit_Q (Q.nrows(), Q.ncols(), Q.get(), Q.lda(), tau);
00167       }
00168 
00169       void
00170       fill_random_svd (const Ordinal nrows, 
00171            const Ordinal ncols, 
00172            Scalar A[], 
00173            const Ordinal lda,
00174            const magnitude_type singular_values[])
00175       {
00176   typedef Matrix< Ordinal, Scalar > matrix_type;
00177   typedef MatView< Ordinal, Scalar > matrix_view_type;
00178 
00179   matrix_type U (nrows, ncols, Scalar(0));
00180   matrix_type V (ncols, ncols, Scalar(0));
00181   std::vector<Scalar> tau_U (std::min (nrows, ncols));
00182   std::vector<Scalar> tau_V (ncols);
00183 
00184   // Fill A with zeros, and then make its diagonal the given set
00185   // of singular values.
00186   matrix_view_type A_view (nrows, ncols, A, lda);
00187   A_view.fill (Scalar (0));
00188   for (Ordinal j = 0; j < ncols; ++j)
00189     // Promote magnitude_type to Scalar here.
00190     A_view(j,j) = Scalar (singular_values[j]);
00191 
00192   // Generate random orthogonal U (nrows by ncols) and V (ncols by
00193   // ncols).  Keep them stored implicitly.
00194   implicit_Q (U, &tau_U[0]);
00195   implicit_Q (V, &tau_V[0]);
00196 
00197   // Workspace query for ORMQR.
00198   Scalar _lwork1, _lwork2;
00199   int info = 0;
00200   LAPACK< Ordinal, Scalar > lapack;
00201   lapack.ORMQR ("L", "N", nrows, ncols, ncols, U.get(), U.lda(), &tau_U[0], 
00202           A, lda, &_lwork1, -1, &info);
00203   if (info != 0)
00204     {
00205       std::ostringstream os;
00206       os << "LAPACK ORMQR LWORK query failed with INFO = " << info 
00207          << ": called ORMQR(\"L\", \"N\", " << nrows << ", " << ncols 
00208          << ", " << ncols << ", NULL, " << U.lda() << ", NULL, NULL, " 
00209          << lda << ", WORK, -1, &INFO)";
00210       throw std::logic_error(os.str());
00211     }
00212   if (ScalarTraits< Scalar >::is_complex)
00213     lapack.ORMQR ("R", "C", nrows, ncols, ncols, V.get(), V.lda(), &tau_V[0], 
00214       A, lda, &_lwork2, -1, &info);
00215   else
00216     lapack.ORMQR ("R", "T", nrows, ncols, ncols, V.get(), V.lda(), &tau_V[0], 
00217       A, lda, &_lwork2, -1, &info);
00218   if (info != 0)
00219     throw std::logic_error("LAPACK ORMQR LWORK query failed");
00220 
00221   // Allocate workspace.
00222   Ordinal lwork = checkedCast (std::max (ScalarTraits< Scalar >::abs (_lwork1), 
00223                  ScalarTraits< Scalar >::abs (_lwork2)));
00224   std::vector< Scalar > work (lwork);
00225 
00226   // Apply U to the left side of A, and V^H to the right side of A.
00227   lapack.ORMQR ("L", "N", nrows, ncols, ncols, U.get(), U.lda(), &tau_U[0], 
00228           A, lda, &work[0], lwork, &info);
00229   if (info != 0)
00230     throw std::runtime_error("LAPACK ORMQR failed (first time)");
00231   if (ScalarTraits< Scalar >::is_complex)
00232     lapack.ORMQR ("R", "C", nrows, ncols, ncols, V.get(), V.lda(), &tau_V[0], 
00233       A, lda, &work[0], lwork, &info);
00234   else
00235     lapack.ORMQR ("R", "T", nrows, ncols, ncols, V.get(), V.lda(), &tau_V[0], 
00236       A, lda, &work[0], lwork, &info);
00237   if (info != 0)
00238     throw std::runtime_error("LAPACK ORMQR failed (second time)");
00239       }
00240 
00241 
00257       void
00258       fill_random_R (const Ordinal n,
00259          Scalar R[],
00260          const Ordinal ldr,
00261          const magnitude_type singular_values[])
00262       {
00263   // Fill R with an n x n (not upper triangular) random matrix
00264   // having the given singular values.
00265   fill_random_svd (n, n, R, ldr, singular_values);
00266 
00267   // Compute the QR factorization in place of R (which isn't upper triangular yet).
00268   std::vector< Scalar > tau (n);
00269 
00270   // Workspace size query for QR factorization.
00271   Scalar _lwork1;
00272   int info = 0;
00273   LAPACK< Ordinal, Scalar > lapack;
00274   lapack.GEQRF (n, n, R, ldr, &tau[0], &_lwork1, -1, &info);
00275   if (info != 0)
00276     throw std::logic_error("LAPACK GEQRF LWORK query failed");
00277 
00278   // Allocate workspace
00279   Ordinal lwork = checkedCast (ScalarTraits< Scalar >::abs (_lwork1));
00280   std::vector< Scalar > work (lwork);
00281 
00282   // Compute QR factorization (implicit representation in place).
00283   lapack.GEQRF (n, n, R, ldr, &tau[0], &work[0], lwork, &info);
00284   if (info != 0)
00285     throw std::runtime_error("LAPACK GEQRF failed");
00286 
00287   // Zero out the stuff below the diagonal of R, leaving just the R factor.
00288   for (Ordinal j = 0; j < n; ++j)
00289     for (Ordinal i = j+1; i < n; ++i)
00290       R[i + j*ldr] = Scalar(0);
00291       }
00292 
00293     private:
00294       static Ordinal 
00295       checkedCast (const magnitude_type& x)
00296       {
00297   if (x < std::numeric_limits< Ordinal >::min() || x > std::numeric_limits< Ordinal >::max())
00298     throw std::range_error("Scalar input cannot be safely cast to an Ordinal");
00299   else if (std::numeric_limits< magnitude_type >::is_signed && 
00300      x < magnitude_type(0) &&
00301      ! std::numeric_limits< Ordinal >::is_signed)
00302     throw std::range_error("Scalar input is negative, but Ordinal is unsigned");
00303   else
00304     return static_cast< Ordinal > (x);
00305       }
00306 
00307       Generator& gen_;
00308     };
00309   } // namespace Random
00310 } // namespace TSQR
00311 
00312 #endif // __TSQR_Random_MatrixGenerator_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends