D.11 valarray

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