/* warning: this small multidimensional matrix library uses a few features not taught in ENGR112 and not explained in elementary textbooks (c) Bjarne Stroustrup, Texas A&M University. Use as you like as long as you acknowledge the source. */ #ifndef MATRIX_LIB #define MATRIX_LIB #include #include //#include namespace Numeric_lib { //----------------------------------------------------------------------------- struct Matrix_error { std::string name; Matrix_error(const char* q) :name(q) { } Matrix_error(std::string n) :name(n) { } }; //----------------------------------------------------------------------------- inline void error(const char* p) { throw Matrix_error(p); } //----------------------------------------------------------------------------- typedef long Index; // I still dislike unsigned //----------------------------------------------------------------------------- // The general Matrix template is simply a prop for its specializations: template class Matrix { // multidimensional matrix class // ( ) does multidimensional subscripting // [ ] does C style "slicing": gives an N-1 dimensional matrix from an N dimensional one // row() is equivalent to [ ] // column() is not (yet) implemented because it requires strides. // = has copy semantics // ( ) and [ ] are range checked // slice() to give sub-ranges private: Matrix(); // this should never be compiled }; //----------------------------------------------------------------------------- template class Row ; // forward declaration //----------------------------------------------------------------------------- // function objects for various apply() operations: template struct Assign { void operator()(T& a, const T& c) { a = c; } }; template struct Add_assign { void operator()(T& a, const T& c) { a += c; } }; template struct Mul_assign { void operator()(T& a, const T& c) { a *= c; } }; template struct Minus_assign { void operator()(T& a, const T& c) { a -= c; } }; template struct Div_assign { void operator()(T& a, const T& c) { a /= c; } }; template struct Mod_assign { void operator()(T& a, const T& c) { a %= c; } }; template struct Or_assign { void operator()(T& a, const T& c) { a |= c; } }; template struct Xor_assign { void operator()(T& a, const T& c) { a ^= c; } }; template struct And_assign { void operator()(T& a, const T& c) { a &= c; } }; template struct Not_assign { void operator()(T& a) { a = !a; } }; template struct Not { T operator()(T& a) { return !a; } }; template struct Unary_minus { T operator()(T& a) { return -a; } }; template struct Complement { T operator()(T& a) { return ~a; } }; //----------------------------------------------------------------------------- // Matrix_base represents the common part of the Matrix classes: template class Matrix_base { // matrixs store their memory (elements) in Matrix_base and have copy semantics // Matrix_base does element-wise operations protected: T* elem; // vector? no: we couldn't easily provide a vector for a slice const Index sz; mutable bool owns; mutable bool xfer; public: Matrix_base(Index n) :elem(new T[n]()), sz(n), owns(true), xfer(false) // matrix of n elements (default initialized) { // std::cerr << "new[" << n << "]->" << elem << "\n"; } Matrix_base(Index n, T* p) :elem(p), sz(n), owns(false), xfer(false) // descriptor for matrix of n elements owned by someone else { } ~Matrix_base() { if (owns) { // std::cerr << "delete[" << sz << "] " << elem << "\n"; delete[]elem; } } // if necessay, we can get to the raw matrix: T* data() { return elem; } const T* data() const { return elem; } Index size() const { return sz; } void copy_elements(const Matrix_base& a) { if (sz!=a.sz) error("copy_elements()"); for (Index i=0; i void base_apply(F f) { for (Index i = 0; i void base_apply(F f, const T& c) { for (Index i = 0; i class Matrix : public Matrix_base { const Index d1; protected: // for use by Row: Matrix(Index n1, T* p) : Matrix_base(n1,p), d1(n1) { // std::cerr << "construct 1D Matrix from data\n"; } public: Matrix(Index n1) : Matrix_base(n1), d1(n1) { } Matrix(Row& a) : Matrix_base(a.dim1(),a.p), d1(a.dim1()) { // std::cerr << "construct 1D Matrix from Row\n"; } // copy constructor: let the base do the copy: Matrix(const Matrix& a) : Matrix_base(a.size(),0), d1(a.d1) { // std::cerr << "copy ctor\n"; this->base_copy(a); } template Matrix(const T (&a)[n]) : Matrix_base(n), d1(n) // deduce "n" (and "T"), Matrix_base allocates T[n] { // std::cerr << "matrix ctor\n"; for (Index i = 0; ielem[i]=a[i]; } Matrix(const T* p, Index n) : Matrix_base(n), d1(n) // Matrix_base allocates T[n] { // std::cerr << "matrix ctor\n"; for (Index i = 0; ielem[i]=p[i]; } template Matrix(const Matrix& a, F f) : Matrix_base(a.size()), d1(a.d1) // construct a new Matrix with element's that are functions of a's elements: // does not modify a unless f has been specifically programmed to modify its argument // T f(const T&) would be a typical type for f { for (Index i = 0; isz; ++i) this->elem[i] = f(a.elem[i]); } template Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base(a.size()), d1(a.d1) // construct a new Matrix with element's that are functions of a's elements: // does not modify a unless f has been specifically programmed to modify its argument // T f(const T&, const Arg&) would be a typical type for f { for (Index i = 0; isz; ++i) this->elem[i] = f(a.elem[i],t1); } Matrix& operator=(const Matrix& a) // copy assignment: let the base do the copy { // std::cerr << "copy assignment (" << this->size() << ',' << a.size()<< ")\n"; if (d1!=a.d1) error("length error in 1D="); this->base_assign(a); return *this; } ~Matrix() { } Index dim1() const { return d1; } // number of elements in a row Matrix xfer() // make an Matrix to move elements out of a scope { Matrix x(dim1(),this->data()); // make a descriptor this->base_xfer(x); // transfer (temporary) ownership to x return x; } void range_check(Index n1) const { // std::cerr << "range check: (" << d1 << "): " << n1 << "\n"; if (n1<0 || d1<=n1) error("1D range error: dimension 1"); } // subscripting: T& operator()(Index n1) { range_check(n1); return this->elem[n1]; } const T& operator()(Index n1) const { range_check(n1); return this->elem[n1]; } // slicing (the same as subscripting for 1D matrixs): T& operator[](Index n) { return row(n); } const T& operator[](Index n) const { return row(n); } T& row(Index n) { range_check(n); return this->elem[n]; } const T& row(Index n) const { range_check(n); return this->elem[n]; } Row slice(Index n) // the last elements from a[n] onwards { if (n<0) n=0; else if(d1(d1-n,this->elem+n); } const Row slice(Index n) const // the last elements from a[n] onwards { if (n<0) n=0; else if(d1(d1-n,this->elem+n); } Row slice(Index n, Index m) // m elements starting with a[n] { if (n<0) n=0; else if(d1(m,this->elem+n); } const Row slice(Index n, Index m) const // m elements starting with a[n] { if (n<0) n=0; else if(d1(m,this->elem+n); } // element-wise operations: template Matrix& apply(F f) { this->base_apply(f); return *this; } template Matrix& apply(F f,const T& c) { this->base_apply(f,c); return *this; } Matrix& operator=(const T& c) { this->base_apply(Assign(),c); return *this; } Matrix& operator*=(const T& c) { this->base_apply(Mul_assign(),c); return *this; } Matrix& operator/=(const T& c) { this->base_apply(Div_assign(),c); return *this; } Matrix& operator%=(const T& c) { this->base_apply(Mod_assign(),c); return *this; } Matrix& operator+=(const T& c) { this->base_apply(Add_assign(),c); return *this; } Matrix& operator-=(const T& c) { this->base_apply(Minus_assign(),c); return *this; } Matrix& operator&=(const T& c) { this->base_apply(And_assign(),c); return *this; } Matrix& operator|=(const T& c) { this->base_apply(Or_assign(),c); return *this; } Matrix& operator^=(const T& c) { this->base_apply(Xor_assign(),c); return *this; } Matrix operator!() { return xfer(Matrix(*this,Not())); } Matrix operator-() { return xfer(Matrix(*this,Unary_minus())); } Matrix operator~() { return xfer(Matrix(*this,Complement())); } template Matrix apply_new(F f) { return xfer(Matrix(*this,f)); } void swap_rows(Index i, Index j) // swap_rows() uses a row's worth of memory for better run-time performance // if you want pairwise swap, just write it yourself { if (i == j) return; /* Matrix temp = (*this)[i]; (*this)[i] = (*this)[j]; (*this)[j] = temp; */ Index max = (*this)[i].size(); for (Index ii=0; ii class Matrix : public Matrix_base { const Index d1; const Index d2; protected: // for use by Row: Matrix(Index n1, Index n2, T* p) : Matrix_base(n1*n2,p), d1(n1), d2(n2) { // std::cerr << "construct 3D Matrix from data\n"; } public: Matrix(Index n1, Index n2) : Matrix_base(n1*n2), d1(n1), d2(n2) { } Matrix(Row& a) : Matrix_base(a.dim1()*a.dim2(),a.p), d1(a.dim1()), d2(a.dim2()) { // std::cerr << "construct 2D Matrix from Row\n"; } // copy constructor: let the base do the copy: Matrix(const Matrix& a) : Matrix_base(a.size(),0), d1(a.d1), d2(a.d2) { // std::cerr << "copy ctor\n"; this->base_copy(a); } template Matrix(const T (&a)[n1][n2]) : Matrix_base(n1*n2), d1(n1), d2(n2) // deduce "n1", "n2" (and "T"), Matrix_base allocates T[n1*n2] { // std::cerr << "matrix ctor\n"; for (Index i = 0; ielem[i*n2+j]=a[i][j]; } template Matrix(const Matrix& a, F f) : Matrix_base(a.size()), d1(a.d1), d2(a.d2) // construct a new Matrix with element's that are functions of a's elements: // does not modify a unless f has been specifically programmed to modify its argument // T f(const T&) would be a typical type for f { for (Index i = 0; isz; ++i) this->elem[i] = f(a.elem[i]); } template Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base(a.size()), d1(a.d1), d2(a.d2) // construct a new Matrix with element's that are functions of a's elements: // does not modify a unless f has been specifically programmed to modify its argument // T f(const T&, const Arg&) would be a typical type for f { for (Index i = 0; isz; ++i) this->elem[i] = f(a.elem[i],t1); } Matrix& operator=(const Matrix& a) // copy assignment: let the base do the copy { // std::cerr << "copy assignment (" << this->size() << ',' << a.size()<< ")\n"; if (d1!=a.d1 || d2!=a.d2) error("length error in 2D ="); this->base_assign(a); return *this; } ~Matrix() { } Index dim1() const { return d1; } // number of elements in a row Index dim2() const { return d2; } // number of elements in a column Matrix xfer() // make an Matrix to move elements out of a scope { Matrix x(dim1(),dim2(),this->data()); // make a descriptor this->base_xfer(x); // transfer (temporary) ownership to x return x; } void range_check(Index n1, Index n2) const { // std::cerr << "range check: (" << d1 << "," << d2 << "): " << n1 << " " << n2 << "\n"; if (n1<0 || d1<=n1) error("2D range error: dimension 1"); if (n2<0 || d2<=n2) error("2D range error: dimension 2"); } // subscripting: T& operator()(Index n1, Index n2) { range_check(n1,n2); return this->elem[n1*d2+n2]; } const T& operator()(Index n1, Index n2) const { range_check(n1,n2); return this->elem[n1*d2+n2]; } // slicing (return a row): Row operator[](Index n) { return row(n); } const Row operator[](Index n) const { return row(n); } Row row(Index n) { range_check(n,0); return Row(d2,&this->elem[n*d2]); } const Row row(Index n) const { range_check(n,0); return Row(d2,&this->elem[n*d2]); } Row slice(Index n) // rows [n:d1) { if (n<0) n=0; else if(d1(d1-n,d2,this->elem+n*d2); } const Row slice(Index n) const // rows [n:d1) { if (n<0) n=0; else if(d1(d1-n,d2,this->elem+n*d2); } Row slice(Index n, Index m) // the rows [n:m) { if (n<0) n=0; if(d1(m-n,d2,this->elem+n*d2); } const Row slice(Index n, Index m) const // the rows [n:sz) { if (n<0) n=0; if(d1(m-n,d2,this->elem+n*d2); } // Column column(Index n); // not (yet) implemented: requies strides and operations on columns // element-wise operations: template Matrix& apply(F f) { this->base_apply(f); return *this; } template Matrix& apply(F f,const T& c) { this->base_apply(f,c); return *this; } Matrix& operator=(const T& c) { this->base_apply(Assign(),c); return *this; } Matrix& operator*=(const T& c) { this->base_apply(Mul_assign(),c); return *this; } Matrix& operator/=(const T& c) { this->base_apply(Div_assign(),c); return *this; } Matrix& operator%=(const T& c) { this->base_apply(Mod_assign(),c); return *this; } Matrix& operator+=(const T& c) { this->base_apply(Add_assign(),c); return *this; } Matrix& operator-=(const T& c) { this->base_apply(Minus_assign(),c); return *this; } Matrix& operator&=(const T& c) { this->base_apply(And_assign(),c); return *this; } Matrix& operator|=(const T& c) { this->base_apply(Or_assign(),c); return *this; } Matrix& operator^=(const T& c) { this->base_apply(Xor_assign(),c); return *this; } Matrix operator!() { return xfer(Matrix(*this,Not())); } Matrix operator-() { return xfer(Matrix(*this,Unary_minus())); } Matrix operator~() { return xfer(Matrix(*this,Complement())); } template Matrix apply_new(F f) { return xfer(Matrix(*this,f)); } void swap_rows(Index i, Index j) // swap_rows() uses a row's worth of memory for better run-time performance // if you want pairwise swap, just write it yourself { if (i == j) return; /* Matrix temp = (*this)[i]; (*this)[i] = (*this)[j]; (*this)[j] = temp; */ Index max = (*this)[i].size(); for (Index ii=0; ii class Matrix : public Matrix_base { const Index d1; const Index d2; const Index d3; protected: // for use by Row: Matrix(Index n1, Index n2, Index n3, T* p) : Matrix_base(n1*n2*n3,p), d1(n1), d2(n2), d3(n3) { // std::cerr << "construct 3D Matrix from data\n"; } public: Matrix(Index n1, Index n2, Index n3) : Matrix_base(n1*n2*n3), d1(n1), d2(n2), d3(n3) { } Matrix(Row& a) : Matrix_base(a.dim1()*a.dim2()*a.dim3(),a.p), d1(a.dim1()), d2(a.dim2()), d3(a.dim3()) { // std::cerr << "construct 3D Matrix from Row\n"; } // copy constructor: let the base do the copy: Matrix(const Matrix& a) : Matrix_base(a.size(),0), d1(a.d1), d2(a.d2), d3(a.d3) { // std::cerr << "copy ctor\n"; this->base_copy(a); } template Matrix(const T (&a)[n1][n2][n3]) : Matrix_base(n1*n2), d1(n1), d2(n2), d3(n3) // deduce "n1", "n2", "n3" (and "T"), Matrix_base allocates T[n1*n2*n3] { // std::cerr << "matrix ctor\n"; for (Index i = 0; ielem[i*n2*n3+j*n3+k]=a[i][j][k]; } template Matrix(const Matrix& a, F f) : Matrix_base(a.size()), d1(a.d1), d2(a.d2), d3(a.d3) // construct a new Matrix with element's that are functions of a's elements: // does not modify a unless f has been specifically programmed to modify its argument // T f(const T&) would be a typical type for f { for (Index i = 0; isz; ++i) this->elem[i] = f(a.elem[i]); } template Matrix(const Matrix& a, F f, const Arg& t1) : Matrix_base(a.size()), d1(a.d1), d2(a.d2), d3(a.d3) // construct a new Matrix with element's that are functions of a's elements: // does not modify a unless f has been specifically programmed to modify its argument // T f(const T&, const Arg&) would be a typical type for f { for (Index i = 0; isz; ++i) this->elem[i] = f(a.elem[i],t1); } Matrix& operator=(const Matrix& a) // copy assignment: let the base do the copy { // std::cerr << "copy assignment (" << this->size() << ',' << a.size()<< ")\n"; if (d1!=a.d1 || d2!=a.d2 || d3!=a.d3) error("length error in 2D ="); this->base_assign(a); return *this; } ~Matrix() { } Index dim1() const { return d1; } // number of elements in a row Index dim2() const { return d2; } // number of elements in a column Index dim3() const { return d3; } // number of elements in a depth Matrix xfer() // make an Matrix to move elements out of a scope { Matrix x(dim1(),dim2(),dim3(),this->data()); // make a descriptor this->base_xfer(x); // transfer (temporary) ownership to x return x; } void range_check(Index n1, Index n2, Index n3) const { // std::cerr << "range check: (" << d1 << "," << d2 << "): " << n1 << " " << n2 << "\n"; if (n1<0 || d1<=n1) error("3D range error: dimension 1"); if (n2<0 || d2<=n2) error("3D range error: dimension 2"); if (n3<0 || d3<=n3) error("3D range error: dimension 3"); } // subscripting: T& operator()(Index n1, Index n2, Index n3) { range_check(n1,n2,n3); return this->elem[d2*d3*n1+d3*n2+n3]; }; const T& operator()(Index n1, Index n2, Index n3) const { range_check(n1,n2,n3); return this->elem[d2*d3*n1+d3*n2+n3]; }; // slicing (return a row): Row operator[](Index n) { return row(n); } const Row operator[](Index n) const { return row(n); } Row row(Index n) { range_check(n,0,0); return Row(d2,d3,&this->elem[n*d2*d3]); } const Row row(Index n) const { range_check(n,0,0); return Row(d2,d3,&this->elem[n*d2*d3]); } Row slice(Index n) // rows [n:d1) { if (n<0) n=0; else if(d1(d1-n,d2,d3,this->elem+n*d2*d3); } const Row slice(Index n) const // rows [n:d1) { if (n<0) n=0; else if(d1(d1-n,d2,d3,this->elem+n*d2*d3); } Row slice(Index n, Index m) // the rows [n:m) { if (n<0) n=0; if(d1(m-n,d2,d3,this->elem+n*d2*d3); } const Row slice(Index n, Index m) const // the rows [n:sz) { if (n<0) n=0; if(d1(m-n,d2,d3,this->elem+n*d2*d3); } // Column column(Index n); // not (yet) implemented: requies strides and operations on columns // element-wise operations: template Matrix& apply(F f) { this->base_apply(f); return *this; } template Matrix& apply(F f,const T& c) { this->base_apply(f,c); return *this; } Matrix& operator=(const T& c) { this->base_apply(Assign(),c); return *this; } Matrix& operator*=(const T& c) { this->base_apply(Mul_assign(),c); return *this; } Matrix& operator/=(const T& c) { this->base_apply(Div_assign(),c); return *this; } Matrix& operator%=(const T& c) { this->base_apply(Mod_assign(),c); return *this; } Matrix& operator+=(const T& c) { this->base_apply(Add_assign(),c); return *this; } Matrix& operator-=(const T& c) { this->base_apply(Minus_assign(),c); return *this; } Matrix& operator&=(const T& c) { this->base_apply(And_assign(),c); return *this; } Matrix& operator|=(const T& c) { this->base_apply(Or_assign(),c); return *this; } Matrix& operator^=(const T& c) { this->base_apply(Xor_assign(),c); return *this; } Matrix operator!() { return xfer(Matrix(*this,Not())); } Matrix operator-() { return xfer(Matrix(*this,Unary_minus())); } Matrix operator~() { return xfer(Matrix(*this,Complement())); } template Matrix apply_new(F f) { return xfer(Matrix(*this,f)); } void swap_rows(Index i, Index j) // swap_rows() uses a row's worth of memory for better run-time performance // if you want pairwise swap, just write it yourself { if (i == j) return; Matrix temp = (*this)[i]; (*this)[i] = (*this)[j]; (*this)[j] = temp; } }; //----------------------------------------------------------------------------- template Matrix scale_and_add(const Matrix& a, T c, const Matrix& b) // Fortran "saxpy()" ("fma" for "fused multiply-add"). // will the copy constructor be called twice and defeat the xfer optimization? { if (a.size() != b.size()) error("sizes wrong for scale_and_add()"); Matrix res(a.size()); for (Index i = 0; i T dot_product(const Matrix&a , const Matrix& b) { if (a.size() != b.size()) error("sizes wrong for dot product"); T sum = 0; for (Index i = 0; i Matrix xfer(Matrix& a) { return a.xfer(); } //----------------------------------------------------------------------------- template A apply(F f, A x) { A res(x,f); return xfer(res); } template A apply(F f, A x, Arg a) { A res(x,f,a); return xfer(res); } //----------------------------------------------------------------------------- // The default values for T and D have been declared before. template class Row { // general version exists only to allow specializations private: Row(); }; //----------------------------------------------------------------------------- template class Row : public Matrix { public: Row(Index n, T* p) : Matrix(n,p) { } Matrix& operator=(const T& c) { this->base_apply(Assign(),c); return *this; } Matrix& operator=(const Matrix& a) { return *static_cast*>(this)=a; } }; //----------------------------------------------------------------------------- template class Row : public Matrix { public: Row(Index n1, Index n2, T* p) : Matrix(n1,n2,p) { } Matrix& operator=(const T& c) { this->base_apply(Assign(),c); return *this; } Matrix& operator=(const Matrix& a) { return *static_cast*>(this)=a; } }; //----------------------------------------------------------------------------- template class Row : public Matrix { public: Row(Index n1, Index n2, Index n3, T* p) : Matrix(n1,n2,n3,p) { } Matrix& operator=(const T& c) { this->base_apply(Assign(),c); return *this; } Matrix& operator=(const Matrix& a) { return *static_cast*>(this)=a; } }; //----------------------------------------------------------------------------- template Matrix scale_and_add(const Matrix& a, const Matrix c, const Matrix& b) { Matrix res(a.size()); if (a.size() != b.size()) error("sizes wrong for scale_and_add"); for (Index i = 0; i Matrix operator*(const Matrix& m, const T& c) { Matrix r(m); return r*=c; } template Matrix operator/(const Matrix& m, const T& c) { Matrix r(m); return r/=c; } template Matrix operator%(const Matrix& m, const T& c) { Matrix r(m); return r%=c; } template Matrix operator+(const Matrix& m, const T& c) { Matrix r(m); return r+=c; } template Matrix operator-(const Matrix& m, const T& c) { Matrix r(m); return r-=c; } template Matrix operator&(const Matrix& m, const T& c) { Matrix r(m); return r&=c; } template Matrix operator|(const Matrix& m, const T& c) { Matrix r(m); return r|=c; } template Matrix operator^(const Matrix& m, const T& c) { Matrix r(m); return r^=c; } //----------------------------------------------------------------------------- } #endif