fmat.h

#pragma once
 
#include "fvec.h"
#include <cassert>
 
namespace cgv {
        namespace math {
 
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
class fmat : public fvec<T,N*M>
{
public:
        typedef fvec<T,N*M> base_type;
        typedef fmat<T,N,M> this_type;
        fmat() {}
        fmat(const T& c) : base_type(c) {}
        fmat(cgv::type::uint32_type n, cgv::type::uint32_type m, const T* a, bool column_major = true) {
                for (cgv::type::uint32_type j = 0; j < std::min(m, M); ++j)
                        for (cgv::type::uint32_type i = 0; i < std::min(n, N); ++i)
                                (*this)(i, j) = a[column_major ? j * n + i : i * m + j];
        }
        template <typename S>
        fmat(cgv::type::uint32_type n, cgv::type::uint32_type m, const S* a, bool column_major = true) {
                for (cgv::type::uint32_type j = 0; j < std::min(m, M); ++j)
                        for (cgv::type::uint32_type i = 0; i < std::min(n, N); ++i)
                                (*this)(i, j) = (T)a[column_major ? j * n + i : i * m + j];
        }
        template <typename S>
        fmat(const fmat<S,N,M>& m) : base_type(m) {}
        template <typename T1, typename T2>
        fmat(const fvec<T1,N>& v, const fvec<T2,M>& w) { 
                for(unsigned i = 0; i < N; i++)
                        for(unsigned j = 0; j < M; j++)
                                (*this)(i,j) = (T)(v(i)*w(j)); 
        }
        static unsigned nrows() { return N; }
        static unsigned ncols() { return M; }
        template <typename S> 
        fmat<T,N,M>& operator = (const fmat<S,N,M>& m) {
                base_type::operator = (m);
                return *this;
        }
        this_type& operator  = (const T& s) { 
                fill (s);
                return *this; 
        }
        bool is_square() const { return N == M; }
        T& operator () (unsigned i, unsigned j) {
                assert(i < N && j < M);
                return base_type::v[j*N+i];
        }
        const T& operator () (unsigned i, unsigned j) const {
                assert(i < N && j < M);
                return base_type::v[j*N+i];
        }
        //in place scalar multiplication
        this_type& operator *= (const T& s)     { base_type::operator *= (s); return *this; } 
        this_type operator * (const T& s) const { this_type r=(*this); r*=(T)s; return r; }
        fmat<T,N,M>& operator /= (const T& s) { base_type::operator /= (s); return *this; }
        const fmat<T,N,M> operator / (const T& s) const { this_type r=(*this); r/=(T)s; return r; }
        fmat<T,N,M>& operator += (const T& s) { base_type::operator += (s); return *this; }
        const fmat<T,N,M> operator + (const T& s) { this_type r=(*this); r+=(T)s; return r; }
        fmat<T,N,M>& operator -= (const T& s) { base_type::operator -= (s); return *this; }
        const fmat<T,N,M> operator - (const T& s) { this_type r=(*this); r-=(T)s; return r; }
        const fmat<T,N,M> operator-() const { return (*this)*(-1); }
        template <typename S> 
        fmat<T,N,M>& operator += (const fmat<S,N,M>& m) { base_type::operator += (m); return *this; } 
        template <typename S> 
        fmat<T,N,M>& operator -= (const fmat<S,N,M>& m) { base_type::operator -= (m); return *this; } 
        template <typename S>
        const fmat<T,N,M> operator+(const fmat<S,N,M> m2)const { fmat<T,N,M> r=(*this); r += m2; return r; }
        template <typename S>
        const fmat<T,N,M> operator-(const fmat<S,N,M> m2)const { fmat<T,N,M> r=(*this); r -= m2; return r; }
        template <typename S>
        const fmat<T,N,M> operator*=(const fmat<S,N,N>& m2) 
        {
                assert(N == M);
                fmat<T,N,N> r(T(0));    
                for(unsigned i = 0; i < N; i++)
                        for(unsigned j = 0; j < N;j++)
                                for(unsigned k = 0; k < N; k++)
                                        r(i,j) += operator()(i,k) * (T)(m2(k,j)); 
                (*this)=r;
                return *this;
        }
        template <typename S, cgv::type::uint32_type L>
        const fmat<T,N,L> operator*(const fmat<S,M,L>& m2) const
        {
                fmat<T,N,L> r; r.zeros();       
                for(unsigned i = 0; i < N; i++)
                        for(unsigned j = 0; j < L;j++)
                                for(unsigned k = 0; k < M; k++)
                                        r(i,j) += operator()(i,k) * (T)(m2(k,j)); 
                return r;
        }
 
        template < typename S>
        const fvec<T,N> operator * (const fvec<S,M>& v) const {
                fvec<T,N> r;
                for(unsigned i = 0; i < N; i++)
                        r(i) = dot(row(i),v);
                return r;
        }
        const fvec<T,M> row(unsigned i) const {
                fvec<T,M> r;
                for(unsigned j = 0; j < M; j++) 
                        r(j)=operator()(i,j);
                return r;
        }
        void set_row(unsigned i,const fvec<T,M>& v) {
                for(unsigned j = 0; j < M;j++) 
                        operator()(i,j)=v(j);           
        }
        const fvec<T,N>& col(unsigned j) const {        
                return *(const fvec<T,N>*)(&operator()(0,j));
        }
        void set_col(unsigned j,const fvec<T,N>& v) {
                for(unsigned i = 0; i < N;i++) 
                        operator()(i,j)=v(i);           
        }
        T trace() const {
                assert(N == M);
                T t = 0;
                for(unsigned i = 0; i < N;i++)
                        t+=operator()(i,i);
                return t;
        }
 
        void transpose() {
                assert(N == M);
                for(unsigned i = 1; i < N; i++)
                        for(unsigned j = 0; j < i; j++)
                                std::swap(operator()(i,j), operator()(j,i));
        }
        T frobenius_norm() const { return base_type::length(); }
        void identity()
        {
                base_type::zeros();
                for(unsigned i = 0; i < M && i < N; ++i)
                        operator()(i,i)=1;
        }
};
 
template <typename T, cgv::type::uint32_type N>
fmat<T,N,N> transpose(const fmat<T,N,N>& m)
{
        fmat<T,N,N> m_t(m);
        m_t.transpose();
        return m_t; 
}
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
fmat<T,N,M> operator * (const T& s, const fmat<T,N,M>& m)
{ 
        return m*s; 
}
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
fvec<T, M> operator * (const fvec<T, N>& v_row, const fmat<T, N, M>& m)
{
        fvec<T, M> r_row;
        for (unsigned i = 0; i < M; i++)
                r_row(i) = dot(m.col(i), v_row);
        return r_row;
}
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
std::ostream& operator<<(std::ostream& out, const fmat<T,N,M>& m)
{
        for (unsigned i=0;i<N;++i) {
                for(unsigned j =0;j < M-1;++j)
                        out << m(i,j) << " ";
                out << m(i,M-1);
                if(i != N-1)
                        out <<"\n";
        }
        return out;
}
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
std::istream& operator>>(std::istream& in, fmat<T,N,M>& m)
{
        for (unsigned i=0;i<m.nrows();++i)
                for(unsigned j =0;j < m.ncols();++j)
                        in >> m(i,j);
        return in;
}
 
 
 
template < typename T, cgv::type::uint32_type N, typename S, cgv::type::uint32_type M>
fmat<T, N, M> dyad(const fvec<T,N>& v, const fvec<S,M>& w)
{
        fmat<T, N, M> m;
        for (unsigned i = 0; i < N; i++)
                for (unsigned j = 0; j < M; j++)
                        m(i, j) = v(i)*(T)w(j);
        return m;
}
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
const fmat<T, N, M> lerp(const fmat<T, N, M>& m1, const fmat<T, N, M>& m2, T t)
{
        fmat<T, N, M> m;
        for(unsigned i = 0; i < N; i++)
                for(unsigned j = 0; j < M; j++)
                        m(i, j) = ((T)1 - t)*m1(i, j) + t * m2(i, j);
        return m;
}
 
template <typename T, cgv::type::uint32_type N, cgv::type::uint32_type M>
const fmat<T, N, M> lerp(const fmat<T, N, M>& m1, const fmat<T, N, M>& m2, fmat<T, N, M> t)
{
        fmat<T, N, M> m;
        for(unsigned i = 0; i < N; i++)
                for(unsigned j = 0; j < M; j++)
                        m(i, j) = ((T)1 - t(i, j))*m1(i, j) + t(i, j)*m2(i, j);
        return m;
}
 
        }
 
}