valarray1.C |
// mytest.C #include <iostream> #include <valarray> #include <slice> using namespace std; #include "Slice_iter.h" #include "Matrix.h" void print_valarray(valarray<double> &a) { int i; cout << "size=" << a.size() << endl; for (i = 0; i < a.size(); i++) cout << a[i] << " "; cout << endl; } double sup_norm(valarray<double> &a) { // valarray<double> t = a.apply(abs); valarray<double> t = abs(a); return t.sum(); } double naiseki(valarray<double> &a, valarray<double> &b) { valarray<double> ab = a * b; return ab.sum(); } int main() { int i, N; valarray<double> y; cin >> N; valarray<double> x(1.0, N), z(N); y.resize(N); for (i = 0; i < N; i++) y[i] = i; cout << "vector x" << endl; print_valarray(x); cout << "vector y" << endl; print_valarray(y); z = x + y; cout << "vector z:=x+y" << endl; print_valarray(z); x *= 0.2; cout << "vector x:=0.2*x" << endl; print_valarray(x); cout << "sup norm of x=" << sup_norm(x) << endl; x = 2; y = 3; cout << "vector x" << endl; print_valarray(x); cout << "vector y" << endl; print_valarray(y); cout << "naiseki of x, y=" << naiseki(x, y) << endl; return 0; } |
Slice_iter.h |
template<class T> class Slice_iter { valarray<T> *v; slice s; size_t curr; 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() { 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()(size_t i) { return ref(curr=i); } T& operator*() { return ref(curr); } // 関係演算子 ==, !=, < を定義 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; } } |
Matrix.h |
#include "Slice_iter.h" class Matrix { valarray<double> *v; size_t d1, d2; public: Matrix(size_t x, size_t y); 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); double operator()(size_t x, size_t y); Slice_iter<double> operator() (size_t i) { return row(i); } Cslice_iter<double> operator() (size_t i) const { return row(i); } Slice_iter<double> operator[] (size_t i) { return row(i); } Cslice_iter<double> operator[] (size_t i) const { return row(i); } Matrix& operator*=(double); valarray<double>& array() { return *v; } } Matrix::Matrix(size_t x, size_t y) { d1 = x; d2 = y; v = new valarray<double>(x*y); } double Matrix::operator()(size_t x, size_t y) { return row(x)[y]; } double mul(const valarray<double> &v1, const valarray<double> &v2) { double res = 0; for (int i = 0; i < v1.size(); i++) res += v1[i] * v2[i]; return res; } valarray<double> operator*(const Matrix &m, const valarray<double> &v) { 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; } 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)); } |