c++ valarray 实现矩阵与向量相乘

#include <iostream>
#include <valarray>

template<class T> class Slice_iter {
　　std::valarray<T>* v;
　　std::slice s;
　　size_t curr;
　　T& ref(size_t i) const {
　　　　return (*v)[s.start()+i*s.stride()];
　　}
public:
　　Slice_iter(std::valarray<T> *vv, std::slice ss):v(vv),s(ss),curr(0) {}
　　Slice_iter end() {
　　　　Slice_iter t=*this;
　　　　t.curr=s.size();
　　　　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(curr=i);}
　　T& operator* () {return ref(curr);}

　　size_t size() const {return s.size();}

　　template<class U>
　　friend bool operator==(const Slice_iter<U>&, const Slice_iter<U>&);
　　template<class U>
　　friend bool operator!=(const Slice_iter<U>&, const Slice_iter<U>&);
　　template<class U>
　　friend bool operator<(const Slice_iter<U>&, const Slice_iter<U>&);
};

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();
}

class Matrix {
　　std::valarray<double>* v;
　　size_t d1, d2;
public:
　　Matrix(size_t x, size_t y);
　　Matrix(const Matrix&);
　　Matrix& operator=(const Matrix&);
　　~Matrix() {
　　　　if(!v)
　　　　　　delete v;
　　}

　　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);
　　Slice_iter<double> column(size_t i);

　　Slice_iter<double> operator[] (size_t i) {return row(i);}

　　Matrix& operator*=(double);

　　std::valarray<double>& array() {return *v;}
};

inline Slice_iter<double> Matrix::row(size_t i) {
　　return Slice_iter<double>(v, std::slice(i*d2,d2,1));
}

inline Slice_iter<double> Matrix::column(size_t i) {
　　return Slice_iter<double>(v, std::slice(i, d1, d2));
}

Matrix::Matrix(size_t x, size_t y) {
　　d1=x;
　　d2=y;
　　v=new std::valarray<double> (x*y);
}

double mul(Slice_iter<double> v1, const std::valarray<double>& v2) {
　　double res=0;
　　for(int i=0; i<v1.size(); ++i)
　　　　res+=v1[i]*v2[i];

　　return res;
}

std::valarray<double> operator*(Matrix& m, const std::valarray<double>& v) {
　　std::valarray<double> res(m.dim1());
　　for(int i=0;i<m.dim1();++i)
　　　　res[i]=mul(m.row(i), v);

　　return res;
}

Matrix& Matrix::operator*=(double d) {
　　(*v) *= d;
　　return *this;
}

size_t slice_index(const std::slice& s, size_t i) {
　　return s.start()+i*s.stride();
}

void print_seq(const std::slice& s) {
　　for(int i=0;i<s.size(); ++i)
　　　　std::cout<<slice_index(s,i)<<" ";
}

int main() {
　　const double vd[] = {1,2,3,4,5,6,7,8};
　　std::valarray<double> v(vd, sizeof(vd)/sizeof(vd[0]));
　　Matrix MyMat(2, 4);
　　MyMat.array()=v;

　　std::cout<<"dim1 = "<<MyMat.dim1()<<" ";
　　std::cout<<"dim2 = "<<MyMat.dim2()<<" ";

　　for(size_t i=0;i<MyMat.dim1();++i) {
　　　　for(size_t j=0;j<MyMat.dim2();++j) {
　　　　　　std::cout<<MyMat.row(i)[j]<<" ";
　　　　}
　　　　std::cout<<" ";
　　}

　　std::slice seg(0, 4, 1);
　　std::valarray<double> bar=v[seg];
　　for(size_t i=0; i<bar.size(); ++i) {
　　　　std::cout<<bar[i]<<" ";
　　}

　　std::cout<<" ";

　　std::cout<<"The product of MyMat and bar is: ";
　　std::valarray<double> res=MyMat*bar;

　　for(size_t i=0;i<res.size();++i) {
　　　　std::cout<<res[i]<<" ";
　　}
　　std::cout<<" ";

　　return 0;
}