Kokkos Node API and Local Linear Algebra Kernels Version of the Day
Tsqr_Matrix.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_Tsqr_Matrix_hpp
00030 #define __TSQR_Tsqr_Matrix_hpp
00031 
00032 #include <Tsqr_Util.hpp>
00033 #include <Tsqr_MatView.hpp>
00034 
00035 #include <stdexcept>
00036 #include <sstream>
00037 #include <limits>
00038 #include <vector>
00039 
00042 
00043 namespace TSQR {
00044 
00058   template<class Ordinal, class Scalar>
00059   class Matrix {
00060   private:
00061     static bool 
00062     fits_in_size_t (const Ordinal& ord) 
00063     {
00064       const Ordinal result = static_cast< Ordinal > (static_cast< size_t > (ord));
00065       return (ord == result);
00066     }
00067 
00079     size_t
00080     verified_alloc_size (const Ordinal num_rows,
00081        const Ordinal num_cols) const 
00082     {
00083       if (! std::numeric_limits< Ordinal >::is_integer)
00084   throw std::logic_error("Ordinal must be an integer type");
00085 
00086       // Quick exit also checks for zero num_cols (which prevents
00087       // division by zero in the tests below).
00088       if (num_rows == 0 || num_cols == 0)
00089   return size_t(0);
00090 
00091       // If Ordinal is signed, make sure that num_rows and num_cols
00092       // are nonnegative.
00093       if (std::numeric_limits< Ordinal >::is_signed)
00094   {
00095     if (num_rows < 0)
00096       {
00097         std::ostringstream os;
00098         os << "# rows (= " << num_rows << ") < 0";
00099         throw std::logic_error (os.str());
00100       }
00101     else if (num_cols < 0)
00102       {
00103         std::ostringstream os;
00104         os << "# columns (= " << num_cols << ") < 0";
00105         throw std::logic_error (os.str());
00106       }
00107   }
00108 
00109       // If Ordinal is bigger than a size_t, do special range
00110       // checking.  The compiler warns (comparison of signed and
00111       // unsigned) if Ordinal is a signed type and we try to do
00112       // "numeric_limits<size_t>::max() <
00113       // std::numeric_limits<Ordinal>::max()", so instead we cast each
00114       // of num_rows and num_cols to size_t and back to Ordinal again,
00115       // and see if we get the same result.  If not, then we
00116       // definitely can't return a size_t product of num_rows and
00117       // num_cols.
00118       if (! fits_in_size_t (num_rows))
00119   {
00120     std::ostringstream os;
00121     os << "# rows (= " << num_rows << ") > max size_t value (= " 
00122        << std::numeric_limits<size_t>::max() << ")";
00123     throw std::range_error (os.str());
00124   }
00125       else if (! fits_in_size_t (num_cols))
00126   {
00127     std::ostringstream os;
00128     os << "# columns (= " << num_cols << ") > max size_t value (= "
00129        << std::numeric_limits<size_t>::max() << ")";
00130     throw std::range_error (os.str());
00131   }
00132 
00133       // Both num_rows and num_cols fit in a size_t, and are
00134       // nonnegative.  Now check whether their product also fits in a
00135       // size_t.  
00136       //
00137       // Note: This may throw a SIGFPE (floating-point exception) if
00138       // num_cols is zero.  Be sure to check first (above).
00139       if (static_cast<size_t>(num_rows) > 
00140     std::numeric_limits<size_t>::max() / static_cast<size_t>(num_cols))
00141   {
00142     std::ostringstream os;
00143     os << "num_rows (= " << num_rows << ") * num_cols (= "
00144        << num_cols << ") > max size_t value (= " 
00145        << std::numeric_limits<size_t>::max() << ")";
00146     throw std::range_error (os.str());
00147   }
00148       return static_cast<size_t>(num_rows) * static_cast<size_t>(num_cols);
00149     }
00150     
00151   public:
00152     typedef Scalar scalar_type;
00153     typedef Ordinal ordinal_type;
00154     typedef Scalar* pointer_type;
00155 
00157     Matrix (const Ordinal num_rows, 
00158       const Ordinal num_cols) :
00159       nrows_ (num_rows),
00160       ncols_ (num_cols),
00161       A_ (verified_alloc_size (num_rows, num_cols))
00162     {}
00163 
00165     Matrix (const Ordinal num_rows,
00166       const Ordinal num_cols,
00167       const Scalar& value) :
00168       nrows_ (num_rows),
00169       ncols_ (num_cols),
00170       A_ (verified_alloc_size (num_rows, num_cols), value)
00171     {}
00172 
00178     Matrix (const Matrix& in) :
00179       nrows_ (in.nrows()),
00180       ncols_ (in.ncols()),
00181       A_ (verified_alloc_size (in.nrows(), in.ncols()))
00182     {
00183       if (! in.empty())
00184   copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
00185     }
00186 
00188     Matrix () : nrows_(0), ncols_(0), A_(0) {}
00189 
00191     ~Matrix () {} 
00192 
00199     template<class MatrixViewType>
00200     Matrix (const MatrixViewType& in) :
00201       nrows_ (in.nrows()),
00202       ncols_ (in.ncols()),
00203       A_ (verified_alloc_size (in.nrows(), in.ncols()))
00204     {
00205       if (A_.size() != 0)
00206   copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
00207     }
00208 
00213     template<class MatrixViewType>
00214     void 
00215     copy (MatrixViewType& B)
00216     {
00217       const typename MatrixViewType::ordinal_type num_rows = B.nrows();
00218       const typename MatrixViewType::ordinal_type num_cols = B.ncols();
00219       if (num_rows != nrows() || num_cols != ncols())
00220   {
00221     std::ostringstream os;
00222     os << "Matrix::Copy: incompatible dimensions: attempt to assign " 
00223        << num_rows << " x " << num_cols << " matrix to an " 
00224        << (nrows()) << " x " << (ncols()) << "matrix";
00225     throw std::logic_error (os.str());
00226   }
00227       copy_matrix (nrows(), ncols(), get(), lda(), B.get(), B.lda());
00228     }
00229 
00231     void 
00232     fill (const Scalar value)
00233     {
00234       fill_matrix (nrows(), ncols(), get(), lda(), value);
00235     }
00236 
00241     Scalar& operator() (const Ordinal i, const Ordinal j) {
00242       return A_[i + j*lda()];
00243     }
00244 
00249     const Scalar& operator() (const Ordinal i, const Ordinal j) const {
00250       return A_[i + j*lda()];
00251     }
00252 
00254     Scalar& operator[] (const Ordinal i) {
00255       return A_[i];
00256     }
00257 
00262     template<class MatrixViewType>
00263     bool operator== (const MatrixViewType& B) const 
00264     {
00265       if (nrows() != B.nrows() || ncols() != B.ncols())
00266   return false;
00267     
00268       typedef typename MatrixViewType::ordinal_type second_ordinal_type;
00269       typedef typename MatrixViewType::scalar_type second_scalar_type;
00270       typedef typename MatrixViewType::pointer_type second_pointer_type;
00271 
00272       const ordinal_type A_nrows = nrows();
00273       const ordinal_type A_lda = lda();
00274       const ordinal_type A_ncols = ncols();
00275       const second_ordinal_type B_lda = B.lda();
00276       const scalar_type* A_j = get();
00277       const second_scalar_type* B_j = B.get();
00278 
00279       for (ordinal_type j = 0; j < A_ncols; ++j, A_j += A_lda, B_j += B_lda)
00280   for (ordinal_type i = 0; i < A_nrows; ++i)
00281     if (A_j[i] != B_j[i])
00282       return false;
00283       return true;
00284     }
00285 
00287     Ordinal nrows() const { return nrows_; }
00288 
00290     Ordinal ncols() const { return ncols_; }
00291 
00293     Ordinal lda() const { return nrows_; }
00294 
00296     bool empty() const { return nrows() == 0 || ncols() == 0; }
00297 
00299     Scalar* 
00300     get() 
00301     { 
00302       if (A_.size() > 0)
00303   return &A_[0]; 
00304       else
00305   return static_cast<Scalar*> (NULL);
00306     }
00307 
00309     const Scalar* 
00310     get() const
00311     { 
00312       if (A_.size() > 0)
00313   return &A_[0]; 
00314       else
00315   return static_cast<const Scalar*> (NULL);
00316     }
00317 
00319     MatView<Ordinal, Scalar> view () {
00320       return MatView<Ordinal, Scalar> (nrows(), ncols(), get(), lda());
00321     }
00322 
00324     ConstMatView<Ordinal, Scalar> const_view () const {
00325       typedef ConstMatView<Ordinal, Scalar> const_view_type;
00326       return const_view_type (nrows(), ncols(), 
00327             const_cast<const Scalar*> (get()), lda());
00328     }
00329 
00340     void
00341     reshape (const Ordinal num_rows, const Ordinal num_cols)
00342     {
00343       if (num_rows == nrows() && num_cols == ncols())
00344   return; // no need to reallocate or do anything else
00345 
00346       const size_t alloc_size = verified_alloc_size (num_rows, num_cols);
00347       nrows_ = num_rows;
00348       ncols_ = num_cols;
00349       A_.resize (alloc_size);
00350     }
00351 
00352   private:
00354     Ordinal nrows_;
00356     Ordinal ncols_;
00362     std::vector<Scalar> A_;
00363   };
00364 
00365 } // namespace TSQR
00366 
00367 #endif // __TSQR_Tsqr_Matrix_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends