---

linalg.h

Go to the documentation of this file.

#ifndef __LINALG_H_
#define __LINALG_H_

#define __LINALG_RANGE_CHECK    0

#include <mtl/mtl.h>
#include <mtl/matrix.h>
#include <mtl/scaled1D.h>
#include <mtl/utils.h>
#include <typeinfo>

template <class _Tp>
class matrix : public mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::row_major>::type
{
  typedef matrix<_Tp> _Self;
  typedef typename mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::row_major>::type _Base;
public:
  typedef typename _Base::OneD Row;

  matrix() : _Base() {}
  matrix(const _Self& __m) : _Base(__m) {}
  matrix(const size_t& n, const size_t& m) : _Base(n,m) {}
};

template <class _Tp>
class c_matrix : public mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::column_major>::type
{
  typedef matrix<_Tp> _Self;
  typedef typename mtl::matrix<_Tp, mtl::rectangle<>, mtl::compressed<>, mtl::column_major>::type _Base;
public:
  typedef typename _Base::OneD Col;

  c_matrix() : _Base() {}
  c_matrix(const _Self& __m) : _Base(__m) {}
  c_matrix(const size_t& n, const size_t& m) : _Base(n,m) {}
};

template <class _Tp>
class sparse_vector : public mtl::compressed1D<_Tp>
{
  typedef sparse_vector<_Tp> _Self;
  typedef mtl::compressed1D<_Tp> _Base;
};

template <class _Tp>
bool operator==(const sparse_vector<_Tp>& a, const sparse_vector<_Tp>& b)
{
  typename sparse_vector<_Tp>::const_iterator _b1, _b2, _e1, _e2;

  _b1 = a.begin();
  _b2 = b.begin();
  _e1 = a.end();
  _e2 = b.end();
  while(_b1 != _e1 || _b2 != _e2)
  {
    if(_b1 == _e1)
    {
      if(*_b2 == _Tp())
      {
        ++_b2;
        continue;
      }
      else
        return false;
    }
    else if(_b2 == _e2)
    {
      if(*_b1 == _Tp())
      {
        ++_b1;
        continue;
      }
      else
        return false;
    }
    else if(_b1.index() < _b2.index())
    {
      if(*_b1 == _Tp())
      {
        ++_b1;
        continue;
      }
      else
        return false;
    }
    else if(_b2.index() < _b1.index())
    {
      if(*_b2 == _Tp())
      {
        ++_b2;
        continue;
      }
      else
        return false;
    }
    else if(*_b1 != *_b2)
      return false;
    ++_b1;
    ++_b2;
  }
  return true;
}

#define linalg_scale(X, Y) mtl::scaled(X, Y)


template <class _S>
class __linalg_cvec : std::vector<_S>
{
private:
  typedef std::vector<_S> _Base;
  typedef __linalg_cvec<_S> _Self;
  const std::vector<_S>& vec;

public:
  enum { N = 0 };

  typedef typename std::vector<_S>::value_type value_type;
  typedef typename std::vector<_S>::reference reference;
  typedef typename std::vector<_S>::pointer pointer;
  typedef typename std::vector<_S>::const_reference const_reference;
  typedef typename std::vector<_S>::const_pointer const_pointer;
  typedef typename std::vector<_S>::size_type size_type;
  typedef typename std::vector<_S>::difference_type difference_type;
  typedef typename std::vector<_S>::allocator_type allocator_type;
  typedef typename std::vector<_S>::iterator iterator;
  typedef typename std::vector<_S>::const_iterator const_iterator;
  typedef typename std::vector<_S>::reverse_iterator reverse_iterator;
  typedef typename std::vector<_S>::const_reverse_iterator const_reverse_iterator;

  typedef mtl::dense_tag sparsity;
  typedef mtl::oned_tag dimension;

  typedef mtl::scaled1D<_Self> scaled_type;

  typedef mtl::light1D<_S> subrange_type;
  typedef _Self IndexArray;
  typedef _Self IndexArrayRef;

  __linalg_cvec(const std::vector<_S>& s) : vec(s) {}
  ~__linalg_cvec() {}

  iterator begin() { return vec.begin(); }
  iterator end() { return vec.end(); }

  const_iterator begin() const { return vec.begin(); }
  const_iterator end() const { return vec.end(); }
};

template <class _S>
class __linalg_vec : std::vector<_S>
{
private:
  typedef std::vector<_S> _Base;
  typedef __linalg_vec<_S> _Self;
  std::vector<_S>& vec;

public:
  enum { N = 0 };

  typedef typename std::vector<_S>::value_type value_type;
  typedef typename std::vector<_S>::reference reference;
  typedef typename std::vector<_S>::pointer pointer;
  typedef typename std::vector<_S>::const_reference const_reference;
  typedef typename std::vector<_S>::const_pointer const_pointer;
  typedef typename std::vector<_S>::size_type size_type;
  typedef typename std::vector<_S>::difference_type difference_type;
  typedef typename std::vector<_S>::allocator_type allocator_type;
  typedef typename std::vector<_S>::iterator iterator;
  typedef typename std::vector<_S>::const_iterator const_iterator;
  typedef typename std::vector<_S>::reverse_iterator reverse_iterator;
  typedef typename std::vector<_S>::const_reverse_iterator const_reverse_iterator;

  typedef mtl::dense_tag sparsity;
  typedef mtl::oned_tag dimension;
  
  typedef mtl::scaled1D<_Self> scaled_type;

  typedef mtl::light1D<_S> subrange_type;
  typedef _Self IndexArray;
  typedef _Self IndexArrayRef;


  __linalg_vec(std::vector<_S>& s) : vec(s) {}
  ~__linalg_vec() {}

  iterator begin() { return vec.begin(); }
  iterator end() { return vec.end(); }

  const_iterator begin() const { return vec.begin(); }
  const_iterator end() const { return vec.end(); }
};

template <class _S>
inline _S linalg_dot(const std::vector<_S>& a, const std::vector<_S>& b, _S s)
{
  return mtl::dot(__linalg_cvec<_S>(a), __linalg_cvec<_S>(b), s);
}

template <class _Va, class _S>
inline _S linalg_dot(const _Va& a, const std::vector<_S>& b, _S s)
{
  __linalg_cvec<_S> _b(b);
  return mtl::dot(a, _b, s);
}

template <class _Vb, class _S>
inline _S linalg_dot(const std::vector<_S>& a, const _Vb& b, _S s)
{
  return mtl::dot(__linalg_cvec<_S>(a), b, s);
}

template <class _Va, class _Vb, class _S>
inline _S linalg_dot(const _Va& a, const _Vb& b, _S s)
{
  return mtl::dot(a, b, s);
}

template <class _S>
inline void linalg_add(const std::vector<_S>& a, const std::vector<_S>& b, std::vector<_S>& c)
{
  mtl::add(__linalg_cvec<_S>(a), __linalg_cvec<_S>(b), __linalg_vec<_S>(c));
}

template <class _Va, class _S>
inline void linalg_add(const _Va& a, const std::vector<_S>& b, std::vector<_S>& c)
{
  mtl::add((a), __linalg_cvec<_S>(b), __linalg_vec<_S>(c));
}

template <class _Vb, class _S>
inline void linalg_add(const std::vector<_S>& a, const _Vb& b, std::vector<_S>& c)
{
  mtl::add(__linalg_cvec<_S>(a), (b), __linalg_vec<_S>(c));
}

template <class _S, class _Vc>
inline void linalg_add(const std::vector<_S>& a, const std::vector<_S>& b, _Vc& c)
{
  mtl::add(__linalg_cvec<_S>(a), __linalg_cvec<_S>(b), (c));
}

template <class _Va, class _Vb, class _S>
inline void linalg_add(const _Va& a, const _Vb& b, std::vector<_S>& c)
{
  mtl::add((a), (b), __linalg_vec<_S>(c));
}

template <class _Va, class _S, class _Vc>
inline void linalg_add(const _Va& a, const std::vector<_S>& b, _Vc& c)
{
  mtl::add((a), __linalg_cvec<_S>(b), (c));
}

template <class _S, class _Vb, class _Vc>
inline void linalg_add(const std::vector<_S>& a, const _Vb& b, _Vc& c)
{
  mtl::add(__linalg_cvec<_S>(a), (b), (c));
}

template <class _Va, class _Vb, class _Vc>
inline void linalg_add(const _Va& a, const _Vb& b, _Vc& c)
{
  mtl::add((a), (b), (c));
}

// matrix vector multiplication
template <class _Matrix, class _S>
inline void linalg_matvec(const _Matrix& A, const std::vector<_S>& a,
                          std::vector<_S>& c)
{
  mtl::mult((A), __linalg_cvec<_S>(a), __linalg_vec<_S>(c));
}

template <class _Matrix, class _S, class _Va>
inline void linalg_matvec(const _Matrix& A, const _Va& a,
                          std::vector<_S>& c)
{
  mtl::mult((A), (a), __linalg_vec<_S>(c));
}

template <class _Matrix, class _S, class _Vc>
inline void linalg_matvec(const _Matrix& A, const std::vector<_S>& a,
                          _Vc& c)
{
  mtl::mult((A), __linalg_cvec<_S>(a), __linalg_cvec<_S>(b), (c));
}

template <class _Matrix, class _Va, class _Vc>
inline void linalg_matvec(const _Matrix& A, const _Va& a, _Vc& c)
{
  mtl::mult((A), (a), (c));
}

template <class _Matrix, class _S>
inline void linalg_matvecv(const _Matrix& A, const std::vector<_S>& a,
                          const std::vector<_S>& b, std::vector<_S>& c)
{
  mtl::mult((A), __linalg_cvec<_S>(a), __linalg_cvec<_S>(b),
            __linalg_vec<_S>(c));
}

template <class _Matrix, class _S, class _Va>
inline void linalg_matvecv(const _Matrix& A, const _Va& a,
                          const std::vector<_S>& b, std::vector<_S>& c)
{
  mtl::mult((A), (a), __linalg_cvec<_S>(b), __linalg_vec<_S>(c));
}

template <class _Matrix, class _S, class _Vb>
inline void linalg_matvecv(const _Matrix& A, const std::vector<_S>& a,
                          const _Vb& b, std::vector<_S>& c)
{
  mtl::mult((A), __linalg_cvec<_S>(a), (b), __linalg_vec<_S>(c));
}

template <class _Matrix, class _S, class _Vc>
inline void linalg_matvecv(const _Matrix& A, const std::vector<_S>& a,
                          const std::vector<_S>& b, _Vc& c)
{
  mtl::mult((A), __linalg_cvec<_S>(a), __linalg_cvec<_S>(b), (c));
}

template <class _Matrix, class _Va, class _Vb, class _S>
inline void linalg_matvecv(const _Matrix& A, _Va& a,
                          const _Vb& b, const std::vector<_S>& c)
{
  mtl::mult((A), (a), (b), __linalg_vec<_S>(c));
}

template <class _Matrix, class _Va, class _S, class _Vc>
inline void linalg_matvecv(const _Matrix& A, const _Va& a,
                          const std::vector<_S>& b, _Vc& c)
{
  mtl::mult((A), (a), __linalg_cvec<_S>(b), (c));
}

template <class _Matrix, class _S, class _Vb, class _Vc>
inline void linalg_matvecv(const _Matrix& A, const std::vector<_S>& a,
                          const _Vb& b, _Vc& c)
{
  mtl::mult((A), __linalg_cvec<_S>(a), (b), (c));
}

template <class _Matrix, class _Va, class _Vb, class _Vc>
inline void linalg_matvecv(const _Matrix& A, const _Va& a, const _Vb& b,
                          _Vc& c)
{
  mtl::mult((A), (a), (b), (c));
}

#define linalg_trans(A) mtl::trans((A))

// other functions
template <class _Ve>
inline std::vector<_Ve>& linalg_ssum(std::vector<_Ve>& a, _Ve b, const std::vector<_Ve>& c)
{
#if __LINALG_RANGE_CHECK
  if(a.size() != b.size())
    throw "linalg: a and c have different size in linalg_ssum";
#endif
  typename std::vector<_Ve>::iterator _a(a.begin()), e_a(a.end());
  typename std::vector<_Ve>::const_iterator _c(c.begin());

  while(_a != e_a)
    *_a++ += b * *_c++;
  return a;
}

template <class _Ve>
inline std::vector<_Ve>& linalg_smult(std::vector<_Ve>& a, _Ve b)
{
  typename std::vector<_Ve>::iterator _a(a.begin()), e_a(a.end());

  while(_a != e_a)
    *_a++ *= b;
  return a;
}

#endif /* __LINALG_H_ */

---