"The C++ Programming Language(3rd and SE)" Source Code

// Program to test slices and a simple N*M matrix class // pp 670-674 and 683-684 // No guarantees offered. Constructive comments to bs@research.att.com #include<iostream> #include<valarray> #include<algorithm> #include<numeric> // for inner_product using namespace std; // forward declarations to allow friend declarations: template<class T> class Slice_iter; template<class T> bool operator==(const Slice_iter<T>&, const Slice_iter<T>&); template<class T> bool operator!=(const Slice_iter<T>&, const Slice_iter<T>&); template<class T> bool operator< (const Slice_iter<T>&, const Slice_iter<T>&); template<class T> class Slice_iter { valarray<T>* v; slice s; size_t curr; // index of current element T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; } public: Slice_iter(valarray<T>* vv, slice ss) :v(vv), s(ss), curr(0) { } Slice_iter end() const { Slice_iter t = *this; t.curr = s.size(); // index of last-plus-one element return t; } Slice_iter& operator++() { curr++; return *this; } Slice_iter operator++(int) { Slice_iter t = *this; curr++; return t; } T& operator[](size_t i) { return ref(i); } // C style subscript T& operator()(size_t i) { return ref(i); } // Fortran-style subscript T& operator*() { return ref(curr); } // current element friend bool operator==<>(const Slice_iter& p, const Slice_iter& q); friend bool operator!=<>(const Slice_iter& p, const Slice_iter& q); friend bool operator< <>(const Slice_iter& p, const Slice_iter& q); }; template<class T> bool operator==(const Slice_iter<T>& p, const Slice_iter<T>& q) { return p.curr==q.curr && p.s.stride()==q.s.stride() && p.s.start()==q.s.start(); } template<class T> bool operator!=(const Slice_iter<T>& p, const Slice_iter<T>& q) { return !(p==q); } template<class T> bool operator<(const Slice_iter<T>& p, const Slice_iter<T>& q) { return p.curr<q.curr && p.s.stride()==q.s.stride() && p.s.start()==q.s.start(); } //------------------------------------------------------------- // forward declarations to allow friend declarations: template<class T> class Cslice_iter; template<class T> bool operator==(const Cslice_iter<T>&, const Cslice_iter<T>&); template<class T> bool operator!=(const Cslice_iter<T>&, const Cslice_iter<T>&); template<class T> bool operator< (const Cslice_iter<T>&, const Cslice_iter<T>&); template<class T> class Cslice_iter { valarray<T>* v; slice s; size_t curr; // index of current element const T& ref(size_t i) const { return (*v)[s.start()+i*s.stride()]; } public: Cslice_iter(valarray<T>* vv, slice ss): v(vv), s(ss), curr(0){} Cslice_iter end() const { Cslice_iter t = *this; t.curr = s.size(); // index of one plus last element return t; } Cslice_iter& operator++() { curr++; return *this; } Cslice_iter operator++(int) { Cslice_iter t = *this; curr++; return t; } const T& operator[](size_t i) const { return ref(i); } const T& operator()(size_t i) const { return ref(i); } const T& operator*() const { return ref(curr); } friend bool operator==<>(const Cslice_iter& p, const Cslice_iter& q); friend bool operator!=<>(const Cslice_iter& p, const Cslice_iter& q); friend bool operator< <>(const Cslice_iter& p, const Cslice_iter& q); }; template<class T> bool operator==(const Cslice_iter<T>& p, const Cslice_iter<T>& q) { return p.curr==q.curr && p.s.stride()==q.s.stride() && p.s.start()==q.s.start(); } template<class T> bool operator!=(const Cslice_iter<T>& p, const Cslice_iter<T>& q) { return !(p==q); } template<class T> bool operator<(const Cslice_iter<T>& p, const Cslice_iter<T>& q) { return p.curr<q.curr && p.s.stride()==q.s.stride() && p.s.start()==q.s.start(); } //------------------------------------------------------------- class Matrix { valarray<double>* v; // stores elements by column as described in 22.4.5 size_t d1, d2; // d1 == number of columns, d2 == number of rows public: Matrix(size_t x, size_t y); // note: no default constructor Matrix(const Matrix&); Matrix& operator=(const Matrix&); ~Matrix(); size_t size() const { return d1*d2; } size_t dim1() const { return d1; } size_t dim2() const { return d2; } Slice_iter<double> row(size_t i); Cslice_iter<double> row(size_t i) const; Slice_iter<double> column(size_t i); Cslice_iter<double> column(size_t i) const; double& operator()(size_t x, size_t y); // Fortran-style subscripts double operator()(size_t x, size_t y) const; Slice_iter<double> operator()(size_t i) { return column(i); } Cslice_iter<double> operator()(size_t i) const { return column(i); } Slice_iter<double> operator[](size_t i) { return column(i); } // C-style subscript Cslice_iter<double> operator[](size_t i) const { return column(i); } Matrix& operator*=(double); valarray<double>& array() { return *v; } }; inline Slice_iter<double> Matrix::row(size_t i) { return Slice_iter<double>(v,slice(i,d1,d2)); } inline Cslice_iter<double> Matrix::row(size_t i) const { return Cslice_iter<double>(v,slice(i,d1,d2)); } inline Slice_iter<double> Matrix::column(size_t i) { return Slice_iter<double>(v,slice(i*d2,d2,1)); } inline Cslice_iter<double> Matrix::column(size_t i) const { return Cslice_iter<double>(v,slice(i*d2,d2,1)); } Matrix::Matrix(size_t x, size_t y) { // check that x and y are sensible d1 = x; d2 = y; v = new valarray<double>(x*y); } Matrix::~Matrix() { delete v; } double& Matrix::operator()(size_t x, size_t y) { return column(x)[y]; } //------------------------------------------------------------- double mul(const Cslice_iter<double>& v1, const valarray<double>& v2) { double res = 0; for (size_t i = 0; i<v2.size(); i++) res+= v1[i]*v2[i]; return res; } valarray<double> operator*(const Matrix& m, const valarray<double>& v) { if (m.dim1()!=v.size()) cerr << "wrong number of elements in m*v\n"; valarray<double> res(m.dim2()); for (size_t i = 0; i<m.dim2(); i++) res[i] = mul(m.row(i),v); return res; } // alternative definition of m*v //valarray<double> operator*(const Matrix& m, valarray<double>& v) valarray<double> mul_mv(const Matrix& m, valarray<double>& v) { if (m.dim1()!=v.size()) cerr << "wrong number of elements in m*v\n"; valarray<double> res(m.dim2()); for (size_t i = 0; i<m.dim2(); i++) { const Cslice_iter<double>& ri = m.row(i); res[i] = inner_product(ri,ri.end(),&v[0],double(0)); } return res; } valarray<double> operator*(valarray<double>& v, const Matrix& m) { if (v.size()!=m.dim2()) cerr << "wrong number of elements in v*m\n"; valarray<double> res(m.dim1()); for (size_t i = 0; i<m.dim1(); i++) { const Cslice_iter<double>& ci = m.column(i); res[i] = inner_product(ci,ci.end(),&v[0],double(0)); } return res; } Matrix& Matrix::operator*=(double d) { (*v) *= d; return *this; } ostream& operator<<(ostream& os, Matrix& m) { for(int y=0; y<m.dim2(); y++) { for(int x=0; x<m.dim1(); x++) os<<m[x][y]<<"\t"; os << "\n"; } return os; } //------------------------------------------------------------- void f (int x_max, int y_max) // test layout and basic access { cout << "\nf(" << x_max <<"," << y_max << "):\n"; Matrix a(x_max, y_max); for(int x=0; x<x_max; x++) // initialize for(int y=0; y<y_max; y++) a[x][y]=x+y*10; cout<<"C-style access used to initialize:\n" << a; for(int x=0; x<x_max; x++) for(int y=0; y<y_max; y++) a(x,y)=x+y*10; cout <<"Fortran-style access used to initialize:\n" << a; cout << "addresses: \n"; for(int x=0; x<x_max; x++) for(int y=0; y<y_max; y++) cout<<"("<< x<<","<<y<<") at "<<&a[x][y]-&a[0][0]<<"\n"; cout <<"columns :\n"; for(int x=0; x<x_max; x++) { cout << "column " << x << ":\n"; for (Slice_iter<double> c = a.column(x); c!=c.end(); ++c) cout << "\t" << *c <<"\n"; } cout <<"rows :\n"; for(int y=0; y<y_max; y++) { cout << "row " << y << ":"; for(Slice_iter<double> r = a.row(y); r!=r.end(); ++r) cout << "\t" << *r ; cout <<"\n"; } } ostream& operator<<(ostream& os,