/* * * Template Numerical Toolkit (TNT) * * Mathematical and Computational Sciences Division * National Institute of Technology, * Gaithersburg, MD USA * * * This software was developed at the National Institute of Standards and * Technology (NIST) by employees of the Federal Government in the course * of their official duties. Pursuant to title 17 Section 105 of the * United States Code, this software is not subject to copyright protection * and is in the public domain. NIST assumes no responsibility whatsoever for * its use by other parties, and makes no guarantees, expressed or implied, * about its quality, reliability, or any other characteristic. * */ // C compatible matrix: row-oriented, 0-based [i][j] and 1-based (i,j) indexing // #ifndef TNT_CMAT_H #define TNT_CMAT_H #include "tnt_subscript.h" #include "tnt_vec.h" #include #include #include #include namespace TNT { template class Matrix { public: typedef Subscript size_type; typedef T value_type; typedef T element_type; typedef T* pointer; typedef T* iterator; typedef T& reference; typedef const T* const_iterator; typedef const T& const_reference; Subscript lbound() const { return 1;} protected: Subscript m_; Subscript n_; Subscript mn_; // total size T* v_; T** row_; T* vm1_ ; // these point to the same data, but are 1-based T** rowm1_; // internal helper function to create the array // of row pointers void initialize(Subscript M, Subscript N) { mn_ = M*N; m_ = M; n_ = N; v_ = new T[mn_]; row_ = new T*[M]; rowm1_ = new T*[M]; assert(v_ != NULL); assert(row_ != NULL); assert(rowm1_ != NULL); T* p = v_; vm1_ = v_ - 1; for (Subscript i=0; i &A) { initialize(A.m_, A.n_); copy(A.v_); } Matrix(Subscript M, Subscript N, const T& value = T()) { initialize(M,N); set(value); } Matrix(Subscript M, Subscript N, const T* v) { initialize(M,N); copy(v); } Matrix(Subscript M, Subscript N, const char *s) { initialize(M,N); //std::istrstream ins(s); std::istringstream ins(s); Subscript i, j; for (i=0; i> row_[i][j]; } // destructor // ~Matrix() { destroy(); } // reallocating // Matrix& newsize(Subscript M, Subscript N) { if (num_rows() == M && num_cols() == N) return *this; destroy(); initialize(M,N); return *this; } // assignments // Matrix& operator=(const Matrix &A) { if (v_ == A.v_) return *this; if (m_ == A.m_ && n_ == A.n_) // no need to re-alloc copy(A.v_); else { destroy(); initialize(A.m_, A.n_); copy(A.v_); } return *this; } Matrix& operator=(const T& scalar) { set(scalar); return *this; } Subscript dim(Subscript d) const { #ifdef TNT_BOUNDS_CHECK assert( d >= 1); assert( d <= 2); #endif return (d==1) ? m_ : ((d==2) ? n_ : 0); } Subscript num_rows() const { return m_; } Subscript num_cols() const { return n_; } inline T* operator[](Subscript i) { #ifdef TNT_BOUNDS_CHECK assert(0<=i); assert(i < m_) ; #endif return row_[i]; } inline const T* operator[](Subscript i) const { #ifdef TNT_BOUNDS_CHECK assert(0<=i); assert(i < m_) ; #endif return row_[i]; } inline reference operator()(Subscript i) { #ifdef TNT_BOUNDS_CHECK assert(1<=i); assert(i <= mn_) ; #endif return vm1_[i]; } inline const_reference operator()(Subscript i) const { #ifdef TNT_BOUNDS_CHECK assert(1<=i); assert(i <= mn_) ; #endif return vm1_[i]; } inline reference operator()(Subscript i, Subscript j) { #ifdef TNT_BOUNDS_CHECK assert(1<=i); assert(i <= m_) ; assert(1<=j); assert(j <= n_); #endif return rowm1_[i][j]; } inline const_reference operator() (Subscript i, Subscript j) const { #ifdef TNT_BOUNDS_CHECK assert(1<=i); assert(i <= m_) ; assert(1<=j); assert(j <= n_); #endif return rowm1_[i][j]; } }; /* *************************** I/O ********************************/ template std::ostream& operator<<(std::ostream &s, const Matrix &A) { Subscript M=A.num_rows(); Subscript N=A.num_cols(); s << M << " " << N << "\n"; for (Subscript i=0; i std::istream& operator>>(std::istream &s, Matrix &A) { Subscript M, N; s >> M >> N; if ( !(M == A.num_rows() && N == A.num_cols() )) { A.newsize(M,N); } for (Subscript i=0; i> A[i][j]; } return s; } // *******************[ basic matrix algorithms ]*************************** template Matrix operator+(const Matrix &A, const Matrix &B) { Subscript M = A.num_rows(); Subscript N = A.num_cols(); assert(M==B.num_rows()); assert(N==B.num_cols()); Matrix tmp(M,N); Subscript i,j; for (i=0; i Matrix operator-(const Matrix &A, const Matrix &B) { Subscript M = A.num_rows(); Subscript N = A.num_cols(); assert(M==B.num_rows()); assert(N==B.num_cols()); Matrix tmp(M,N); Subscript i,j; for (i=0; i Matrix mult_element(const Matrix &A, const Matrix &B) { Subscript M = A.num_rows(); Subscript N = A.num_cols(); assert(M==B.num_rows()); assert(N==B.num_cols()); Matrix tmp(M,N); Subscript i,j; for (i=0; i Matrix transpose(const Matrix &A) { Subscript M = A.num_rows(); Subscript N = A.num_cols(); Matrix S(N,M); Subscript i, j; for (i=0; i inline Matrix matmult(const Matrix &A, const Matrix &B) { #ifdef TNT_BOUNDS_CHECK assert(A.num_cols() == B.num_rows()); #endif Subscript M = A.num_rows(); Subscript N = A.num_cols(); Subscript K = B.num_cols(); Matrix tmp(M,K); T sum; for (Subscript i=0; i inline Matrix operator*(const Matrix &A, const Matrix &B) { return matmult(A,B); } template inline int matmult(Matrix& C, const Matrix &A, const Matrix &B) { assert(A.num_cols() == B.num_rows()); Subscript M = A.num_rows(); Subscript N = A.num_cols(); Subscript K = B.num_cols(); C.newsize(M,K); T sum; const T* row_i; const T* col_k; for (Subscript i=0; i Vector matmult(const Matrix &A, const Vector &x) { #ifdef TNT_BOUNDS_CHECK assert(A.num_cols() == x.dim()); #endif Subscript M = A.num_rows(); Subscript N = A.num_cols(); Vector tmp(M); T sum; for (Subscript i=0; i inline Vector operator*(const Matrix &A, const Vector &x) { return matmult(A,x); } } // namespace TNT #endif // CMAT_H