mex_matrix_separate.hpp

// This file is a part of NiHu, a C++ BEM template library.
//
// Copyright (C) 2012-2014  Peter Fiala <fiala@hit.bme.hu>
// Copyright (C) 2012-2014  Peter Rucz <rucz@hit.bme.hu>
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
#ifndef MEX_MATRIX_SEPARATE_HPP_INCLUDED
#define MEX_MATRIX_SEPARATE_HPP_INCLUDED
 
#include "eigen_utils.hpp"
#include "../core/result_matrix.hpp"
 
#include <complex>
#include <cstddef>
#include <type_traits>
 
namespace NiHu
{
 
namespace mex {
 
class matrix_base
{
protected:
    matrix_base(size_t rows, size_t cols)
        : m_rows(rows), m_cols(cols)
    {
    }
 
    matrix_base(mxArray const *input)
        : m_rows(mxGetM(input))
        , m_cols(mxGetN(input))
    {
    }
 
public:
    size_t rows(void) const 
    {
        return m_rows;
    }
 
    size_t cols(void) const
    {
        return m_cols;
    }
 
protected:
    size_t m_rows;
    size_t m_cols;
};
 
 
template <class T>
struct classID;
 
template <>
struct classID<int>
{
    static mxClassID const value = mxINT32_CLASS;
};
 
template <>
struct classID<double>
{
    static mxClassID const value = mxDOUBLE_CLASS;
};
 
template <>
struct classID<float>
{
    static mxClassID const value = mxSINGLE_CLASS;
};
 
 
template <class T>
class real_matrix :
    public Eigen::Map<Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> >
{
public:
    typedef Eigen::Map< Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> > base_t;
 
    typedef T scalar_t;
 
    real_matrix(size_t rows, size_t cols, mxArray *&output)
        : base_t(
            static_cast<scalar_t *>(mxGetData(output = mxCreateNumericMatrix(rows, cols, classID<scalar_t>::value, mxREAL))),
            rows, cols)
    {
    }
 
    real_matrix(mxArray const *input)
        : base_t(static_cast<scalar_t *>(mxGetData(input)), mxGetM(input), mxGetN(input))
    {
    }
};
 
template <class Parent>
class index_proxy
{
public:
    typedef typename Parent::scalar_t scalar_t;
 
    index_proxy(Parent &matrix, size_t row, size_t col)
        : m_parent(matrix), m_row(row), m_col(col)
    {
    }
 
    operator std::complex<scalar_t>() const
    {
        return std::complex<scalar_t>(real(), imag());
    }
 
    scalar_t &real(void) const
    {
        return m_parent.m_real[m_row + m_parent.rows()*m_col];
    }
 
    scalar_t &imag(void) const
    {
        return m_parent.m_imag[m_row + m_parent.rows()*m_col];
    }
 
    template <class complex_rhs_t>
    void operator =(complex_rhs_t const &data) const
    {
        real() = data.real();
        imag() = data.imag();
    }
 
    template <class complex_rhs_t>
    void operator +=(complex_rhs_t const &data) const
    {
        real() += data.real();
        imag() += data.imag();
    }
 
    void operator =(scalar_t const &data)
    {
        real() = data;
        imag() = 0.0;
    }
 
    void operator +=(scalar_t const &data)
    {
        real() += data;
    }
 
private:
    Parent &m_parent;
    size_t const m_row;
    size_t const m_col;
};
 
 
template <class T>
class complex_matrix :
    public matrix_base
{
public:
    friend class index_proxy<complex_matrix<T> >;
    typedef T scalar_t;
 
    complex_matrix(size_t rows, size_t cols, mxArray *&output)
        : matrix_base(rows, cols)
    {
        output = mxCreateNumericMatrix(
            m_rows, m_cols, classID<scalar_t>::value, mxCOMPLEX);
        m_real = static_cast<scalar_t *>(mxGetData(output));
        m_imag = static_cast<scalar_t *>(mxGetImagData(output));
    }
 
    complex_matrix(mxArray const *input)
        : matrix_base(input),
        m_real(static_cast<scalar_t *>(mxGetData(input))),
        m_imag(static_cast<scalar_t *>(mxGetImagData(input)))
    {
    }
 
    std::complex<scalar_t> operator() (size_t row, size_t col) const
    {
        return std::complex<scalar_t>(
            m_real[row+m_rows*col], m_imag[row+m_rows*col]);
    }
 
    index_proxy<complex_matrix> operator() (size_t row, size_t col)
    {
        return index_proxy<complex_matrix>(*this, row, col);
    }
 
protected:
    scalar_t *m_real;
    scalar_t *m_imag;
};
 
template <class Scalar>
Scalar const &get_scalar(mxArray const *pa)
{
    return *static_cast<Scalar *>(mxGetPr(pa));
}
 
 
} // end of namespace mex
 
 
template <class T>
struct is_result_matrix_impl<mex::real_matrix<T> > : std::true_type {};
 
template <class T>
struct is_result_matrix_impl<mex::complex_matrix<T> > : std::true_type {};
 
} // end of namespace NiHu
 
#endif /* MEX_MATRIX_SEPARATE_HPP_INCLUDED */