include/util.h

Go to the documentation of this file.

#ifndef UTIL_H
#define UTIL_H

//== INCLUDES =================================================================

#include <cassert>
#include <iostream>

#ifdef  __CUDA_ARCH__
#   ifdef assert
#       undef assert
#   endif
#   define assert (void)
#endif

//== NAMESPACES ===============================================================

namespace gpufilter {

//== CLASS DEFINITION =========================================================

template <class T, int N>
class Vector {
    
public:

    std::vector<T> to_vector() const {
        return std::vector<T>(&m_data[0], &m_data[0]+size());
    }

    __host__ __device__ const T& operator [] ( int i ) const {
        assert(i >= 0 && i < size());
        return m_data[i];
    }

    __host__ __device__ T& operator [] ( int i ) {
        assert(i >= 0 && i < size());
        return m_data[i];
    }

    __host__ __device__ int size() const { return N; }

    friend std::ostream& operator << ( std::ostream& out,
                                       const Vector& v ) {
        out << '[';
        for (int i=0; i<v.size(); ++i) {
            out << v[i];
            if (i < v.size()-1) out << ',';
        }
        return out << ']';
    }

    __host__ __device__ Vector& operator += ( const Vector& v ) {
        assert(size() == v.size());
#pragma unroll
        for (int j=0; j<size(); ++j)
            m_data[j] += v[j];
        return *this;
    }

    __host__ __device__ Vector operator + ( const Vector& v ) const {
        return Vector(*this) += v;
    }

    __host__ __device__ Vector &operator *= ( const T& v ) {
#pragma unroll
        for (int j=0; j<size(); ++j)
            m_data[j] *= v;
        return *this;
    }

    __host__ __device__ Vector operator * ( const T& v ) const {
        return Vector(*this) *= v;
    }

    __host__ __device__ Vector& operator /= ( const T& v ) {
#pragma unroll
        for (int j=0; j<size(); ++j)
            m_data[j] /= v;
        return *this;
    }

    __host__ __device__ Vector operator / ( const T& v ) const {
        return Vector(*this) /= v;
    }

    __host__ __device__ operator const T* () const { return &m_data[0]; }

    __host__ __device__ operator T* () { return &m_data[0]; }

private:

    T m_data[N]; 

};

//== CLASS DEFINITION =========================================================

template <class T, int M, int N=M>
class Matrix {

public:

    __host__ __device__ int rows() const { return M; }

    __host__ __device__ int cols() const { return N; }

    __host__ __device__ const Vector<T,N>& operator [] ( int i ) const {
        assert(i >= 0 && i < rows());
        return m_data[i];
    }

    __host__ __device__ Vector<T,N>& operator [] ( int i ) {
        assert(i >= 0 && i < rows());
        return m_data[i];
    }

    friend std::ostream& operator << ( std::ostream& out,
                                       const Matrix& m ) {
        out << '[';
        for (int i=0; i<m.rows(); ++i) {
            for (int j=0; j<m.cols(); ++j) {
                out << m[i][j];
                if (j < m.cols()-1) out << ',';
            }
            if (i < m.rows()-1) out << ";\n";
        }
        return out << ']';
    }

    template <int P, int Q>
    __host__ __device__ Matrix<T,M,Q> operator * ( const Matrix<T,P,Q>& rhs ) const {
        assert(cols()==rhs.rows());
        Matrix<T,M,Q> r;
        for (int i=0; i<r.rows(); ++i) {
            for (int j=0; j<r.cols(); ++j) {
                r[i][j] = m_data[i][0]*rhs[0][j];
                for (int k=1; k<cols(); ++k)
                    r[i][j] += m_data[i][k]*rhs[k][j];
            }
        }
        return r;
    }

    __host__ __device__ Matrix& operator *= ( T val ) {
#pragma unroll
        for (int i=0; i<rows(); ++i)
#pragma unroll
            for (int j=0; j<cols(); ++j)
                m_data[i][j] *= val;
        return *this;
    }

    __host__ __device__ Vector<T,M> col( int j ) const {
        Vector<T,M> c;
#pragma unroll
        for (int i=0; i<rows(); ++i)
            c[i] = m_data[i][j];
        return c;
    }

    __host__ __device__ void set_col( int j,
                                      const Vector<T,M>& c ) {
#pragma unroll
        for (int i=0; i<rows(); ++i)
            m_data[i][j] = c[i];
    }

    __host__ __device__ friend Matrix operator * ( const Matrix& m,
                                                   T val ) {
        return Matrix(m) *= val;
    }

    __host__ __device__ friend Matrix operator * ( T val,
                                                   const Matrix& m ) {
        return operator*(m,val);
    }

    __host__ __device__ Matrix& operator += ( const Matrix& rhs ) {
#pragma unroll
        for (int i=0; i<rows(); ++i)
#pragma unroll
            for (int j=0; j<cols(); ++j)
                m_data[i][j] += rhs[i][j];
        return *this;
    }

    __host__ __device__ friend Matrix operator + ( const Matrix& lhs,
                                                   const Matrix& rhs ) {
        return Matrix(lhs) += rhs;
    }

    __host__ __device__ Matrix& operator -= ( const Matrix& rhs ) {
#pragma unroll
        for (int i=0; i<rows(); ++i)
#pragma unroll
            for (int j=0; j<cols(); ++j)
                m_data[i][j] -= rhs[i][j];
        return *this;
    }

    __host__ __device__ friend Matrix operator - ( const Matrix& lhs,
                                                   const Matrix& rhs ) {
        return Matrix(lhs) -= rhs;
    }

private:

    Vector<T,N> m_data[M]; 

};

//== IMPLEMENTATION ===========================================================

//-- Multiply -----------------------------------------------------------------

template <class T, int M, int N>
__host__ __device__ Vector<T,N> operator * ( const Vector<T,M>& v,
                                             const Matrix<T,M,N>& m ) {
    assert(v.size() == m.rows());
    Vector<T,N> r;
#pragma unroll
    for (int j=0; j<m.cols(); ++j) {
        r[j] = v[0]*m[0][j];
#pragma unroll
        for (int i=1; i<m.rows(); ++i)
            r[j] += v[i]*m[i][j];
    }
    return r;
}

//-- Basics -------------------------------------------------------------------

template <class T, int N>
__host__ __device__ Vector<T,N> zeros() {
    Vector<T,N> v;
    if (N>0)
    {
#if __CUDA_ARCH__
#pragma unroll
        for (int j=0; j<v.size(); ++j)
            v[j] = T();
#else
        std::fill(&v[0], &v[N-1]+1, T());
#endif
    }
    return v; // I'm hoping that RVO will kick in
}

template <class T, int M, int N>
__host__ __device__ Matrix<T,M,N> zeros() {
    Matrix<T,M,N> mat;
    if (M>0 && N>0)
    {
#if __CUDA_ARCH__
#pragma unroll
        for (int i=0; i<mat.rows(); ++i)
#pragma unroll
            for (int j=0; j<mat.cols(); ++j)
                mat[i][j] = T();
#else
        std::fill(&mat[0][0], &mat[M-1][N-1]+1, T());
#endif
    }
    return mat; // I'm hoping that RVO will kick in
}

template <class T, int M, int N>
Matrix<T,M,N> identity() {
    Matrix<T,M,N> mat;
    for (int i=0; i<M; ++i)
        for (int j=0; j<N; ++j)
            mat[i][j] = i==j ? 1 : 0;
    return mat;
}

template <class T, int M, int N>
__host__ __device__ Matrix<T,N,M> transp( const Matrix<T,M,N>& m ) {
    Matrix<T,N,M> tm;
#pragma unroll
    for (int i=0; i<m.rows(); ++i)
#pragma unroll
        for (int j=0; j<m.cols(); ++j)
            tm[j][i] = m[i][j];
    return tm;
}

//-- Forward ------------------------------------------------------------------

template <class T, int N, int R>
void fwd_inplace( const Vector<T,R>& p,
                  Vector<T,N>& b,
                  const Vector<T,R+1>& w ) {
    for (int j=0; j<b.size(); ++j) {
        b[j] *= w[0];
        for (int k=1; k<w.size(); ++k) {
            if (j-k < 0)
                b[j] -= p[p.size()+j-k]*w[k]; // use data from prologue
            else
                b[j] -= b[j-k]*w[k];
        }
    }
}

template <class T, int M, int N, int R>
void fwd_inplace( const Matrix<T,M,R>& p,
                  Matrix<T,M,N>& b,
                  const Vector<T,R+1>& w ) {
    for (int i=0; i<b.rows(); ++i)
        fwd_inplace(p[i], b[i], w);
}

template <class T, int M, int N, int R>
Matrix<T,M,N> fwd( const Matrix<T,M,R>& p,
                   const Matrix<T,M,N>& b, 
                   const Vector<T,R+1>& w) {
    Matrix<T,M,N> fb = b;
    fwd_inplace(p, fb, w);
    return fb;
}

template <class T, int M, int N, int R>
void fwd_inplace( const Matrix<T,R,N>& p,
                  Matrix<T,M,N>& b,
                  const Vector<T,R+1>& w ) {
    b = fwd(p, b, w);
}

template <class T, int M, int N, int R>
Matrix<T,M,N> fwdT( const Matrix<T,R,N>& pT,
                    const Matrix<T,M,N>& b, 
                    const Vector<T,R+1>& w ) {
    return transp(fwd(transp(pT), transp(b), w));
}

//-- Reverse ------------------------------------------------------------------

template <class T, int N, int R>
void rev_inplace( Vector<T,N>& b,
                  const Vector<T,R>& e,
                  const Vector<T,R+1>& w ) {
    for (int j=b.size()-1; j>=0; --j) {
        b[j] *= w[0];
        for (int k=1; k<w.size(); ++k) {
            if (j+k >= b.size())
                b[j] -= e[j+k-b.size()]*w[k]; // use data from epilogue
            else
                b[j] -= b[j+k]*w[k];
        }
    }
}

template <class T, int M, int N, int R>
void rev_inplace( Matrix<T,M,N>& b,
                  const Matrix<T,M,R>& e,
                  const Vector<T,R+1>& w ) {
    for (int i=0; i<b.rows(); ++i)
        rev_inplace(b[i], e[i], w);
}

template <class T, int M, int N, int R>
Matrix<T,M,N> rev( const Matrix<T,M,N>& b,
                   const Matrix<T,M,R>& e, 
                   const Vector<T,R+1>& w ) {
    Matrix<T,M,N> rb = b;
    rev_inplace(rb, e, w);
    return rb;
}

template <class T, int M, int N, int R>
void rev_inplace( Matrix<T,M,N>& b,
                  const Matrix<T,R,N>& e,
                  const Vector<T,R+1>& w ) {
    b = rev(b, e, w);
}

template <class T, int M, int N, int R>
Matrix<T,M,N> revT( const Matrix<T,M,N>& b,
                    const Matrix<T,R,N>& eT,
                    const Vector<T,R+1>& w ) {
    return transp(rev(transp(b), transp(eT), w));
}

//-- Head ---------------------------------------------------------------------

template <int R, int M, int N, class T>
Matrix<T,M,R> head( const Matrix<T,M,N>& mat ) {
    assert(mat.cols() >= R);
    Matrix<T,M,R> h;
    for (int j=0; j<R; ++j)
        for (int i=0; i<mat.rows(); ++i)
            h[i][j] = mat[i][j];
    return h;
}

template <int R, int M, int N, class T>
Matrix<T,R,N> headT( const Matrix<T,M,N>& mat ) {
    return transp(head<R>(transp(mat)));
}

//-- Tail ---------------------------------------------------------------------

template <int R, int M, int N, class T>
Matrix<T,M,R> tail( const Matrix<T,M,N>& mat ) {
    assert(mat.cols() >= R);
    Matrix<T,M,R> t;
    for (int j=0; j<R; ++j)
        for (int i=0; i<mat.rows(); ++i)
            t[i][j] = mat[i][mat.cols()-R+j];
    return t;
}

template <int R, int M, int N, class T>
Matrix<T,R,N> tailT( const Matrix<T,M,N>& mat ) {
    return transp(tail<R>(transp(mat)));
}

//=============================================================================
} // namespace gpufilter
//=============================================================================
#endif // UTIL_H
//=============================================================================