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 (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_Tsqr_Matrix_hpp
00043 #define __TSQR_Tsqr_Matrix_hpp
00044 
00045 #include <Tsqr_Util.hpp>
00046 #include <Tsqr_MatView.hpp>
00047 
00048 #include <stdexcept>
00049 #include <sstream>
00050 #include <limits>
00051 #include <vector>
00052 
00055 
00056 namespace TSQR {
00057 
00071   template<class Ordinal, class Scalar>
00072   class Matrix {
00073   private:
00074     static bool 
00075     fits_in_size_t (const Ordinal& ord) 
00076     {
00077       const Ordinal result = static_cast< Ordinal > (static_cast< size_t > (ord));
00078       return (ord == result);
00079     }
00080 
00092     size_t
00093     verified_alloc_size (const Ordinal num_rows,
00094        const Ordinal num_cols) const 
00095     {
00096       if (! std::numeric_limits< Ordinal >::is_integer)
00097   throw std::logic_error("Ordinal must be an integer type");
00098 
00099       // Quick exit also checks for zero num_cols (which prevents
00100       // division by zero in the tests below).
00101       if (num_rows == 0 || num_cols == 0)
00102   return size_t(0);
00103 
00104       // If Ordinal is signed, make sure that num_rows and num_cols
00105       // are nonnegative.
00106       if (std::numeric_limits< Ordinal >::is_signed)
00107   {
00108     if (num_rows < 0)
00109       {
00110         std::ostringstream os;
00111         os << "# rows (= " << num_rows << ") < 0";
00112         throw std::logic_error (os.str());
00113       }
00114     else if (num_cols < 0)
00115       {
00116         std::ostringstream os;
00117         os << "# columns (= " << num_cols << ") < 0";
00118         throw std::logic_error (os.str());
00119       }
00120   }
00121 
00122       // If Ordinal is bigger than a size_t, do special range
00123       // checking.  The compiler warns (comparison of signed and
00124       // unsigned) if Ordinal is a signed type and we try to do
00125       // "numeric_limits<size_t>::max() <
00126       // std::numeric_limits<Ordinal>::max()", so instead we cast each
00127       // of num_rows and num_cols to size_t and back to Ordinal again,
00128       // and see if we get the same result.  If not, then we
00129       // definitely can't return a size_t product of num_rows and
00130       // num_cols.
00131       if (! fits_in_size_t (num_rows))
00132   {
00133     std::ostringstream os;
00134     os << "# rows (= " << num_rows << ") > max size_t value (= " 
00135        << std::numeric_limits<size_t>::max() << ")";
00136     throw std::range_error (os.str());
00137   }
00138       else if (! fits_in_size_t (num_cols))
00139   {
00140     std::ostringstream os;
00141     os << "# columns (= " << num_cols << ") > max size_t value (= "
00142        << std::numeric_limits<size_t>::max() << ")";
00143     throw std::range_error (os.str());
00144   }
00145 
00146       // Both num_rows and num_cols fit in a size_t, and are
00147       // nonnegative.  Now check whether their product also fits in a
00148       // size_t.  
00149       //
00150       // Note: This may throw a SIGFPE (floating-point exception) if
00151       // num_cols is zero.  Be sure to check first (above).
00152       if (static_cast<size_t>(num_rows) > 
00153     std::numeric_limits<size_t>::max() / static_cast<size_t>(num_cols))
00154   {
00155     std::ostringstream os;
00156     os << "num_rows (= " << num_rows << ") * num_cols (= "
00157        << num_cols << ") > max size_t value (= " 
00158        << std::numeric_limits<size_t>::max() << ")";
00159     throw std::range_error (os.str());
00160   }
00161       return static_cast<size_t>(num_rows) * static_cast<size_t>(num_cols);
00162     }
00163     
00164   public:
00165     typedef Scalar scalar_type;
00166     typedef Ordinal ordinal_type;
00167     typedef Scalar* pointer_type;
00168 
00170     Matrix (const Ordinal num_rows, 
00171       const Ordinal num_cols) :
00172       nrows_ (num_rows),
00173       ncols_ (num_cols),
00174       A_ (verified_alloc_size (num_rows, num_cols))
00175     {}
00176 
00178     Matrix (const Ordinal num_rows,
00179       const Ordinal num_cols,
00180       const Scalar& value) :
00181       nrows_ (num_rows),
00182       ncols_ (num_cols),
00183       A_ (verified_alloc_size (num_rows, num_cols), value)
00184     {}
00185 
00191     Matrix (const Matrix& in) :
00192       nrows_ (in.nrows()),
00193       ncols_ (in.ncols()),
00194       A_ (verified_alloc_size (in.nrows(), in.ncols()))
00195     {
00196       if (! in.empty())
00197   copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
00198     }
00199 
00201     Matrix () : nrows_(0), ncols_(0), A_(0) {}
00202 
00204     ~Matrix () {} 
00205 
00212     template<class MatrixViewType>
00213     Matrix (const MatrixViewType& in) :
00214       nrows_ (in.nrows()),
00215       ncols_ (in.ncols()),
00216       A_ (verified_alloc_size (in.nrows(), in.ncols()))
00217     {
00218       if (A_.size() != 0)
00219   copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
00220     }
00221 
00226     template<class MatrixViewType>
00227     void 
00228     copy (MatrixViewType& B)
00229     {
00230       const typename MatrixViewType::ordinal_type num_rows = B.nrows();
00231       const typename MatrixViewType::ordinal_type num_cols = B.ncols();
00232       if (num_rows != nrows() || num_cols != ncols())
00233   {
00234     std::ostringstream os;
00235     os << "Matrix::Copy: incompatible dimensions: attempt to assign " 
00236        << num_rows << " x " << num_cols << " matrix to an " 
00237        << (nrows()) << " x " << (ncols()) << "matrix";
00238     throw std::logic_error (os.str());
00239   }
00240       copy_matrix (nrows(), ncols(), get(), lda(), B.get(), B.lda());
00241     }
00242 
00244     void 
00245     fill (const Scalar value)
00246     {
00247       fill_matrix (nrows(), ncols(), get(), lda(), value);
00248     }
00249 
00254     Scalar& operator() (const Ordinal i, const Ordinal j) {
00255       return A_[i + j*lda()];
00256     }
00257 
00262     const Scalar& operator() (const Ordinal i, const Ordinal j) const {
00263       return A_[i + j*lda()];
00264     }
00265 
00267     Scalar& operator[] (const Ordinal i) {
00268       return A_[i];
00269     }
00270 
00275     template<class MatrixViewType>
00276     bool operator== (const MatrixViewType& B) const 
00277     {
00278       if (nrows() != B.nrows() || ncols() != B.ncols())
00279   return false;
00280     
00281       typedef typename MatrixViewType::ordinal_type second_ordinal_type;
00282       typedef typename MatrixViewType::scalar_type second_scalar_type;
00283       typedef typename MatrixViewType::pointer_type second_pointer_type;
00284 
00285       const ordinal_type A_nrows = nrows();
00286       const ordinal_type A_lda = lda();
00287       const ordinal_type A_ncols = ncols();
00288       const second_ordinal_type B_lda = B.lda();
00289       const scalar_type* A_j = get();
00290       const second_scalar_type* B_j = B.get();
00291 
00292       for (ordinal_type j = 0; j < A_ncols; ++j, A_j += A_lda, B_j += B_lda)
00293   for (ordinal_type i = 0; i < A_nrows; ++i)
00294     if (A_j[i] != B_j[i])
00295       return false;
00296       return true;
00297     }
00298 
00300     Ordinal nrows() const { return nrows_; }
00301 
00303     Ordinal ncols() const { return ncols_; }
00304 
00306     Ordinal lda() const { return nrows_; }
00307 
00309     bool empty() const { return nrows() == 0 || ncols() == 0; }
00310 
00312     Scalar* 
00313     get() 
00314     { 
00315       if (A_.size() > 0)
00316   return &A_[0]; 
00317       else
00318   return static_cast<Scalar*> (NULL);
00319     }
00320 
00322     const Scalar* 
00323     get() const
00324     { 
00325       if (A_.size() > 0)
00326   return &A_[0]; 
00327       else
00328   return static_cast<const Scalar*> (NULL);
00329     }
00330 
00332     MatView<Ordinal, Scalar> view () {
00333       return MatView<Ordinal, Scalar> (nrows(), ncols(), get(), lda());
00334     }
00335 
00337     ConstMatView<Ordinal, Scalar> const_view () const {
00338       typedef ConstMatView<Ordinal, Scalar> const_view_type;
00339       return const_view_type (nrows(), ncols(), 
00340             const_cast<const Scalar*> (get()), lda());
00341     }
00342 
00353     void
00354     reshape (const Ordinal num_rows, const Ordinal num_cols)
00355     {
00356       if (num_rows == nrows() && num_cols == ncols())
00357   return; // no need to reallocate or do anything else
00358 
00359       const size_t alloc_size = verified_alloc_size (num_rows, num_cols);
00360       nrows_ = num_rows;
00361       ncols_ = num_cols;
00362       A_.resize (alloc_size);
00363     }
00364 
00365   private:
00367     Ordinal nrows_;
00369     Ordinal ncols_;
00375     std::vector<Scalar> A_;
00376   };
00377 
00378 } // namespace TSQR
00379 
00380 #endif // __TSQR_Tsqr_Matrix_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends