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