Program Listing for File matrix.hpp

↰ Return to documentation for file (source/pele/matrix.hpp)

#ifndef _PELE_MATRIX_H_
#define _PELE_MATRIX_H_
#include "pele/array.hpp"
namespace pele {

template <typename dtype>
class MatrixAdapter : public pele::Array<dtype> {
 public:
  size_t _dim2;

  MatrixAdapter(size_t dim1, size_t dim2, dtype val = 0)
      : pele::Array<dtype>(dim1 * dim2, val), _dim2(dim2) {}

  MatrixAdapter(pele::Array<double> v, size_t dim2)
      : pele::Array<dtype>(v), _dim2(dim2) {
    if (v.size() % dim2 != 0) {
      throw std::invalid_argument("v.size() is not divisible by dim2");
    }
  }

  MatrixAdapter(double *data, size_t dim1, size_t dim2)
      : pele::Array<dtype>(data, dim1 * dim2), _dim2(dim2) {}

  inline dtype const &operator()(size_t i, size_t j) const {
    return this->operator[](i * _dim2 + j);
  }
  inline dtype &operator()(size_t i, size_t j) {
    return this->operator[](i * _dim2 + j);
  }

  inline std::pair<size_t, size_t> shape() const {
    return std::pair<size_t, size_t>(this->size() / _dim2, _dim2);
  }
};

template <class dtype>
MatrixAdapter<dtype> hacky_mat_mul(MatrixAdapter<dtype> const &A,
                                   MatrixAdapter<dtype> const &B) {
  assert(A.shape().second == B.shape().first);
  size_t const L = A.shape().second;
  size_t const N = A.shape().first;
  size_t const M = B.shape().second;

  MatrixAdapter<dtype> C(N, M, 0);
  for (size_t i = 0; i < N; ++i) {
    for (size_t j = 0; j < M; ++j) {
      double val = 0;
      for (size_t k = 0; k < L; ++k) {
        val += A(i, k) * B(k, j);
      }
      C(i, j) = val;
    }
  }
  return C;
}

// Inactive matrix-vector multiplication prototype.
// template<class dtype>
// pele::Array<dtype> hacky_mat_mul(MatrixAdapter<dtype> const & A,
// pele::Array<dtype> const & v)
//{
//    assert(A.shape().second == v.size());
//    size_t const L = A.shape().second;
//    size_t const n = A.shape().first;
//
//    pele::Array<dtype> C(n, 0);
//    for (size_t i = 0; i<n; ++i){
//        dtype val = 0;
//        for (size_t k = 0; k<L; ++k){
//            val += A(i,k) * v[k];
//        }
//        C(i) = val;
//    }
//    return C;
//}

template <class dtype>
std::ostream &operator<<(std::ostream &out,
                         const pele::MatrixAdapter<dtype> &a) {
  out << "[ ";
  size_t const N = a.shape().first;
  size_t const M = a.shape().second;
  for (size_t n = 0; n < N; ++n) {
    for (size_t m = 0; m < M; ++m) {
      if (m > 0) out << ", ";
      out << a(n, m);
    }
    if (n < N - 1) out << ",\n  ";
  }
  out << " ]";
  return out;
}

}  // namespace pele
#endif