Anasazi Version of the Day
Tsqr_Matrix.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Anasazi: Block Eigensolvers Package
00005 //                 Copyright (2010) 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 
00055   template< class Ordinal, class Scalar >
00056   class Matrix {
00057   private:
00058     static bool 
00059     fits_in_size_t (const Ordinal& ord) 
00060     {
00061       const Ordinal result = static_cast< Ordinal > (static_cast< size_t > (ord));
00062       return (ord == result);
00063     }
00064 
00076     size_t
00077     verified_alloc_size (const Ordinal num_rows,
00078        const Ordinal num_cols) const 
00079     {
00080       if (! std::numeric_limits< Ordinal >::is_integer)
00081   throw std::logic_error("Ordinal must be an integer type");
00082 
00083       // Quick exit also checks for zero num_cols (which prevents
00084       // division by zero in the tests below).
00085       if (num_rows == 0 || num_cols == 0)
00086   return size_t(0);
00087 
00088       // If Ordinal is signed, make sure that num_rows and num_cols
00089       // are nonnegative.
00090       if (std::numeric_limits< Ordinal >::is_signed)
00091   {
00092     if (num_rows < 0)
00093       {
00094         std::ostringstream os;
00095         os << "# rows (= " << num_rows << ") < 0";
00096         throw std::logic_error (os.str());
00097       }
00098     else if (num_cols < 0)
00099       {
00100         std::ostringstream os;
00101         os << "# columns (= " << num_cols << ") < 0";
00102         throw std::logic_error (os.str());
00103       }
00104   }
00105 
00106       // If Ordinal is bigger than a size_t, do special range
00107       // checking.  The compiler warns (comparison of signed and
00108       // unsigned) if Ordinal is a signed type and we try to do
00109       // "numeric_limits<size_t>::max() <
00110       // std::numeric_limits<Ordinal>::max()", so instead we cast each
00111       // of num_rows and num_cols to size_t and back to Ordinal again,
00112       // and see if we get the same result.  If not, then we
00113       // definitely can't return a size_t product of num_rows and
00114       // num_cols.
00115       if (! fits_in_size_t (num_rows))
00116   {
00117     std::ostringstream os;
00118     os << "# rows (= " << num_rows << ") > max size_t value (= " 
00119        << std::numeric_limits<size_t>::max() << ")";
00120     throw std::range_error (os.str());
00121   }
00122       else if (! fits_in_size_t (num_cols))
00123   {
00124     std::ostringstream os;
00125     os << "# columns (= " << num_cols << ") > max size_t value (= "
00126        << std::numeric_limits<size_t>::max() << ")";
00127     throw std::range_error (os.str());
00128   }
00129 
00130       // Both num_rows and num_cols fit in a size_t, and are
00131       // nonnegative.  Now check whether their product also fits in a
00132       // size_t.  
00133       //
00134       // Note: This may throw a SIGFPE (floating-point exception) if
00135       // num_cols is zero.  Be sure to check first (above).
00136       if (static_cast<size_t>(num_rows) > 
00137     std::numeric_limits<size_t>::max() / static_cast<size_t>(num_cols))
00138   {
00139     std::ostringstream os;
00140     os << "num_rows (= " << num_rows << ") * num_cols (= "
00141        << num_cols << ") > max size_t value (= " 
00142        << std::numeric_limits<size_t>::max() << ")";
00143     throw std::range_error (os.str());
00144   }
00145       return static_cast<size_t>(num_rows) * static_cast<size_t>(num_cols);
00146     }
00147     
00148   public:
00149     typedef Scalar scalar_type;
00150     typedef Ordinal ordinal_type;
00151     typedef Scalar* pointer_type;
00152 
00153     Matrix (const Ordinal num_rows, 
00154       const Ordinal num_cols) :
00155       nrows_ (num_rows),
00156       ncols_ (num_cols),
00157       A_ (verified_alloc_size (num_rows, num_cols))
00158     {}
00159 
00160     Matrix (const Ordinal num_rows,
00161       const Ordinal num_cols,
00162       const Scalar& value) :
00163       nrows_ (num_rows),
00164       ncols_ (num_cols),
00165       A_ (verified_alloc_size (num_rows, num_cols), value)
00166     {}
00167 
00168     // We need an explicit copy constructor, because for some reason
00169     // the default copy constructor (with shallow copies of pointers,
00170     // eeek! double free()s!!!) overrides the generic "copy
00171     // constructors" below.
00172     Matrix (const Matrix& in) :
00173       nrows_ (in.nrows()),
00174       ncols_ (in.ncols()),
00175       A_ (verified_alloc_size (in.nrows(), in.ncols()))
00176     {
00177       if (A_.size() > 0)
00178   copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
00179     }
00180 
00181     Matrix () : nrows_(0), ncols_(0), A_(0) {}
00182 
00183     ~Matrix () {} 
00184 
00185     template< class MatrixViewType >
00186     Matrix (const MatrixViewType& in) :
00187       nrows_ (in.nrows()),
00188       ncols_ (in.ncols()),
00189       A_ (verified_alloc_size (in.nrows(), in.ncols()))
00190     {
00191       if (A_.size() != 0)
00192   copy_matrix (nrows(), ncols(), get(), lda(), in.get(), in.lda());
00193     }
00194 
00199     template< class MatrixViewType >
00200     void 
00201     copy (MatrixViewType& B)
00202     {
00203       const typename MatrixViewType::ordinal_type num_rows = B.nrows();
00204       const typename MatrixViewType::ordinal_type num_cols = B.ncols();
00205       if (num_rows != nrows() || num_cols != ncols())
00206   {
00207     std::ostringstream os;
00208     os << "Matrix::Copy: incompatible dimensions: attempt to assign " 
00209        << num_rows << " x " << num_cols << " matrix to an " 
00210        << (nrows()) << " x " << (ncols()) << "matrix";
00211     throw std::logic_error (os.str());
00212   }
00213       copy_matrix (nrows(), ncols(), get(), lda(), B.get(), B.lda());
00214     }
00215 
00216     void 
00217     fill (const Scalar value)
00218     {
00219       fill_matrix (nrows(), ncols(), get(), lda(), value);
00220     }
00221 
00225     Scalar& operator() (const Ordinal i, const Ordinal j) {
00226       return A_[i + j*lda()];
00227     }
00228 
00229     const Scalar& operator() (const Ordinal i, const Ordinal j) const {
00230       return A_[i + j*lda()];
00231     }
00232 
00234     Scalar& operator[] (const Ordinal i) {
00235       return A_[i];
00236     }
00237 
00238     template< class MatrixViewType >
00239     bool operator== (const MatrixViewType& B) const 
00240     {
00241       if (nrows() != B.nrows() || ncols() != B.ncols())
00242   return false;
00243     
00244       typedef typename MatrixViewType::ordinal_type second_ordinal_type;
00245       typedef typename MatrixViewType::scalar_type second_scalar_type;
00246       typedef typename MatrixViewType::pointer_type second_pointer_type;
00247 
00248       const ordinal_type A_nrows = nrows();
00249       const ordinal_type A_lda = lda();
00250       const ordinal_type A_ncols = ncols();
00251       const second_ordinal_type B_lda = B.lda();
00252       const scalar_type* A_j = get();
00253       const second_scalar_type* B_j = B.get();
00254 
00255       for (ordinal_type j = 0; j < A_ncols; ++j, A_j += A_lda, B_j += B_lda)
00256   for (ordinal_type i = 0; i < A_nrows; ++i)
00257     if (A_j[i] != B_j[i])
00258       return false;
00259       return true;
00260     }
00261 
00262     Ordinal nrows() const { return nrows_; }
00263     Ordinal ncols() const { return ncols_; }
00264     Ordinal lda() const { return nrows_; }
00265     bool empty() const { return nrows() == 0 || ncols() == 0; }
00266 
00267     Scalar* 
00268     get() 
00269     { 
00270       if (A_.size() > 0)
00271   return &A_[0]; 
00272       else
00273   return static_cast< Scalar* > (NULL);
00274     }
00275     const Scalar* 
00276     get() const
00277     { 
00278       if (A_.size() > 0)
00279   return &A_[0]; 
00280       else
00281   return static_cast< const Scalar* > (NULL);
00282     }
00283 
00284     MatView< Ordinal, Scalar > view () {
00285       return MatView< Ordinal, Scalar >(nrows(), ncols(), get(), lda());
00286     }
00287 
00288     ConstMatView< Ordinal, Scalar > const_view () const {
00289       return ConstMatView< Ordinal, Scalar >(nrows(), ncols(), 
00290                (const Scalar*) get(), lda());
00291     }
00292 
00298     void
00299     reshape (const Ordinal num_rows, const Ordinal num_cols)
00300     {
00301       if (num_rows == nrows() && num_cols == ncols())
00302   return; // no need to reallocate or do anything else
00303 
00304       const size_t alloc_size = verified_alloc_size (num_rows, num_cols);
00305       nrows_ = num_rows;
00306       ncols_ = num_cols;
00307       A_.resize (alloc_size);
00308     }
00309 
00310   private:
00311     Ordinal nrows_, ncols_;
00312     std::vector< Scalar > A_;
00313   };
00314 
00315 } // namespace TSQR
00316 
00317 #endif // __TSQR_Tsqr_Matrix_hpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends