Matrix.h

Go to the documentation of this file.

#ifndef   math_Matrix_H__
#define   math_Matrix_H__

#include "../Self.h"

#include <vector>
#include <algorithm>

#include <cstdlib>

namespace meow {
template<class Entry>
class Matrix {
public:
  typedef typename std::vector<Entry>::reference       EntryRef ;
  typedef typename std::vector<Entry>::const_reference EntryRefK;
private:
  struct Myself {
    size_t                rows_;
    size_t                cols_;
    std::vector<Entry> entries_;

    Myself():
    rows_(0), cols_(0), entries_(0) {
    }

    Myself(Myself const& b):
    rows_(b.rows_), cols_(b.cols_), entries_(b.entries_) {
    }

    Myself(size_t r, size_t c, Entry const& e):
    rows_(r), cols_(c), entries_(r * c, e) {
    }

    ~Myself() {
    }

    size_t index(size_t r, size_t c) const {
      return r * cols_ + c;
    }
  };

  Self<Myself> const self;
public:
  Matrix(): self() { }

  Matrix(Matrix const& m): self(m.self, Self<Myself>::COPY_FROM) {
  }

  Matrix(size_t r, size_t c, Entry const& e): self(Myself(r, c, e)) {
  }

  ~Matrix() { }

  Matrix& copyFrom(Matrix const& m) {
    self().copyFrom(m.self);
    return *this;
  }

  Matrix& referenceFrom(Matrix const& m) {
    self().referenceFrom(m.self);
    return *this;
  }

  void reset(size_t r, size_t c, Entry const& e) {
    self()->rows_ = r;
    self()->cols_ = c;
    self()->entries_.clear();
    self()->entries_.resize(r * c, e);
  }

  bool valid() const {
    return (rows() > 0 && cols() > 0);
  }

  size_t rows() const {
    return self->rows_;
  }

  size_t cols() const {
    return self->cols_;
  }

  size_t size() const {
    return rows() * cols();
  }

  size_t rows(size_t r, Entry const& e) {
    if (r != rows()) {
      self()->entries_.resize(r * cols(), e);
      self()->rows_ = r;
    }
    return rows();
  }

  size_t cols(size_t c, Entry const& e) {
    if (c != cols()) {
      Self<Myself> const old(self, Self<Myself>::COPY_FROM);
      self()->entries_.resize(rows() * c);
      self()->cols_ = c;
      for (size_t i = 0, I = rows(); i < I; i++) {
        size_t j, J1 = std::min(old->cols_, cols()), J2 = cols();
        for (j =  0; j < J1; j++)
          self()->entries_[self->index(i, j)] = old->entries_[old->index(i, j)];
        for (j = J1; j < J2; j++)
          self()->entries_[self->index(i, j)] = e;
      }
    }
    return cols();
  }

  size_t size(size_t r, size_t c, Entry const& e) {
    cols(c, e);
    rows(r, e);
    return rows() * cols();
  }

  Entry entry(size_t r, size_t c) const {
    return self->entries_[self->index(r, c)];
  }

  Entry entry(size_t r, size_t c, Entry const& e) {
    self()->entries_[self->index(r, c)] = e;
    return entry(r, c);
  }

  EntryRef entryGet(size_t r, size_t c) {
    return self()->entries_[self->index(r, c)];
  }

  void entries(ssize_t rFirst, ssize_t rLast,
               ssize_t cFirst, ssize_t cLast,
               Entry const& e) {
    for (ssize_t r = rFirst; r <= rLast; r++) {
      for (ssize_t c = cFirst; c <=cFirst; c++) {
        entry(r, c, e);
      }
    }
  }

  Matrix subMatrix(size_t rFirst, size_t rLast,
                   size_t cFirst, size_t cLast) const {
    if (rFirst > rLast || cFirst > cLast) return Matrix();
    if (rFirst == 0 && cFirst == 0) {
      Matrix ret(*this);
      ret.size(rLast + 1, cLast + 1, Entry(0));
      return ret;
    }
    Matrix ret(rLast - rFirst + 1, cLast - cFirst + 1, entry(rFirst, cFirst));
    for (size_t r = rFirst; r <= rLast; r++)
      for (size_t c = cFirst; c <= cLast; c++)
        ret.entry(r - rFirst, c - cFirst, entry(r, c));
    return ret;
  }

  Matrix row(size_t r) const {
    return subMatrix(r, r, 0, cols() - 1);
  }

  Matrix col(size_t c) const {
    return subMatrix(0, rows() - 1, c, c);
  }

  Matrix positive() const {
    return *this;
  }

  Matrix negative() const {
    Matrix ret(*this);
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = cols(); c < C; c++)
        ret.entry(r, c, -ret.entry(r, c));
    return ret;
  }

  Matrix add(Matrix const& m) const {
    if (rows() != m.rows() || cols() != m.cols()) return Matrix();
    Matrix ret(*this);
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = cols(); c < C; c++)
        ret.entry(r, c, ret.entry(r, c) + m.entry(r, c));
    return ret;
  }

  Matrix sub(Matrix const& m) const {
    if (rows() != m.rows() || cols() != m.cols()) return Matrix();
    Matrix ret(*this);
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = cols(); c < C; c++)
        ret.entry(r, c, ret.entry(r, c) - m.entry(r, c));
    return ret;
  }

  Matrix mul(Matrix const& m) const {
    if (cols() != m.rows()) return Matrix();
    Matrix ret(rows(), m.cols(), Entry(0));
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = m.cols(); c < C; c++)
        for (size_t k = 0, K = cols(); k < K; k++)
          ret.entry(r, c, ret.entry(r, c) + entry(r, k) * m.entry(k, c));
    return ret;
  }

  Matrix mul(Entry const& s) const {
    Matrix ret(*this);
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = cols(); c < C; c++)
        ret.entry(r, c, ret.entry(r, c) * s);
    return ret;
  }

  Matrix div(Entry const& s) const {
    Matrix ret(*this);
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = cols(); c < C; c++)
        ret.entry(r, c, ret.entry(r, c) / s);
    return ret;
  }

  Matrix identity() const {
    Matrix ret(*this);
    ret.identitied();
    return ret;
  }

  Matrix& identitied() {
    for (size_t r = 0, R = rows(); r < R; r++)
      for (size_t c = 0, C = cols(); c < C; c++)
        entry(r, c, (r == c ? Entry(1) : Entry(0)));
    return *this;
  }

  Matrix& diagonaled() {
    triangulared();
    for (size_t i = 0, I = rows(); i < I; ++i) {
      for (size_t j = i + 1, J = cols(); j < J; ++j) {
        entry(i, j, Entry(0));
      }
    }
    return *this;
  }

  Matrix diagonal() const {
    Matrix ret(*this);
    ret.diagonaled();
    return ret;
  }

  Matrix inverse() const {
    if (rows() != cols() || rows() == 0) return Matrix<Entry>();
    Matrix tmp(rows(), cols() * 2, Entry(0));
    for (size_t r = 0, R = rows(); r < R; r++) {
      for (size_t c = 0, C = cols(); c < C; c++) {
        tmp.entry(r, c, entry(r, c));
        tmp.entry(r, c + cols(), (r == c ? Entry(1) : Entry(0)));
      }
    }
    tmp.triangulared();
    for (ssize_t r = rows() - 1; r >= 0; r--) {
      if (tmp(r, r) == Entry(0)) return Matrix<Entry>();
      for (ssize_t r2 = r - 1; r2 >= 0; r2--) {
        Entry rat(-tmp.entry(r2, r) / tmp.entry(r, r));
        for (size_t c = r, C = tmp.cols(); c < C; c++) {
          tmp.entry(r2, c, tmp.entry(r2, c) + rat * tmp(r, c));
        }
      }
      Entry rat(tmp.entry(r, r));
      for (size_t c = cols(), C = tmp.cols(); c < C; c++) {
        tmp.entry(r, c - cols(), tmp.entry(r, c) / rat);
      }
    }
    tmp.size(cols(), rows(), Entry(0));
    return tmp;
  }

  Matrix& inversed() {
    copyFrom(inverse());
    return *this;
  }

  Matrix transpose() const {
    Matrix ret(cols(), rows(), Entry(0));
    for (size_t r = 0, R = cols(); r < R; r++)
      for (size_t c = 0, C = rows(); c < C; c++)
        ret.entry(r, c, entry(c, r));
    return ret;
  }

  Matrix& transposed() {
    copyFrom(transpose());
    return *this;
  }

  Matrix triangular() const {
    Matrix<Entry> ret(*this);
    ret.triangulared();
    return ret;
  }

  Matrix& triangulared() {
    for (size_t r = 0, c = 0, R = rows(), C = cols(); r < R && c < C; r++) {
      ssize_t maxR;
      for ( ; c < C; c++) {
        maxR = -1;
        for (size_t r2 = r; r2 < R; r2++)
          if (maxR == -1 || tAbs(entry(r2, c)) > tAbs(entry(maxR, c)))
            maxR = r2;
        if (entry(maxR, c) != Entry(0)) break;
      }
      if (c >= C) break;
      if (maxR != (ssize_t)r) {
        for (size_t c2 = c; c2 < C; c2++)
          std::swap(self()->entries_[self->index(   r, c2)],
                    self()->entries_[self->index(maxR, c2)]);
      }
      for (size_t r2 = r + 1; r2 < R; r2++) {
        Entry rati = -entry(r2, c) / entry(r, c);
        entry(r2, c, Entry(0));
        for (size_t c2 = c + 1; c2 < C; c2++)
          entry(r2, c2, entry(r2, c2) + entry(r, c2) * rati);
      }
    }
    return *this;
  }

  Matrix& operator=(Matrix const& m) {
    return copyFrom(m);
  }

  Entry operator()(size_t r, size_t c) const {
    return entry(r, c);
  }

  Entry operator()(size_t r, size_t c, Entry const& e) {
    return entry(r, c, e);
  }

  Matrix operator+() const {
    return positive();
  }

  Matrix operator-() const {
    return negative();
  }

  Matrix operator+(Matrix const& m) const {
    return add(m);
  }

  Matrix operator-(Matrix const& m) const {
    return sub(m);
  }

  Matrix operator*(Matrix const& m) const {
    return mul(m);
  }

  Matrix operator*(Entry const& s) const {
    return mul(s);
  }

  Matrix operator/(Entry const& s) const {
    return div(s);
  }
};

} // meow

#endif // math_Matrix_H__