/var/lib/jenkins/jobs/Skylark/workspace/sketch/RFT_Elemental.hpp

Go to the documentation of this file.

#ifndef SKYLARK_RFT_ELEMENTAL_HPP
#define SKYLARK_RFT_ELEMENTAL_HPP

namespace skylark {
namespace sketch {

template <typename ValueType,
          template <typename> class InputType,
          template <typename> class KernelDistribution>
struct RFT_t <
    InputType<ValueType>,
    elem::Matrix<ValueType>,
    KernelDistribution> :
        public RFT_data_t<KernelDistribution> {
    // Typedef value, matrix, transform, distribution and transform data types
    // so that we can use them regularly and consistently.
    typedef ValueType value_type;
    typedef InputType<value_type> matrix_type;
    typedef elem::Matrix<value_type> output_matrix_type;
    typedef RFT_data_t<KernelDistribution> data_type;
private:
    typedef skylark::sketch::dense_transform_t <matrix_type,
                                                output_matrix_type,
                                                KernelDistribution>
    underlying_t;


protected:
    RFT_t (int N, int S, base::context_t& context)
        : data_type (N, S, context) {

    }

public:
    RFT_t(const RFT_t<matrix_type,
                      output_matrix_type,
                      KernelDistribution>& other)
        : data_type(other) {

    }

    RFT_t(const data_type& other_data)
        : data_type(other_data) {

    }

    template <typename Dimension>
    void apply (const matrix_type& A,
                output_matrix_type& sketch_of_A,
                Dimension dimension) const {
        try {
            apply_impl(A, sketch_of_A, dimension);
        } catch (std::logic_error e) {
            SKYLARK_THROW_EXCEPTION (
                base::elemental_exception()
                    << base::error_msg(e.what()) );
        } catch(boost::mpi::exception e) {
            SKYLARK_THROW_EXCEPTION (
                base::mpi_exception()
                    << base::error_msg(e.what()) );
        }
    }

private:
    void apply_impl(const matrix_type& A,
        output_matrix_type& sketch_of_A,
        skylark::sketch::columnwise_tag tag) const {

        // TODO verify sizes etc.
        underlying_t underlying(*data_type::_underlying_data);
        underlying.apply(A, sketch_of_A, tag);

#       if SKYLARK_HAVE_OPENMP
#       pragma omp parallel for collapse(2)
#       endif
        for(int j = 0; j < base::Width(A); j++)
            for(int i = 0; i < data_type::_S; i++) {
                value_type x = sketch_of_A.Get(i, j);
                x += data_type::_shifts[i];

#               ifdef SKYLARK_EXACT_COSINE
                x = std::cos(x);
#               else
                // x = std::cos(x) is slow
                // Instead use low-accuracy approximation
                if (x < -3.14159265) x += 6.28318531;
                else if (x >  3.14159265) x -= 6.28318531;
                x += 1.57079632;
                if (x >  3.14159265)
                    x -= 6.28318531;
                x = (x < 0) ?
                    1.27323954 * x + 0.405284735 * x * x :
                    1.27323954 * x - 0.405284735 * x * x;
#               endif

                x = data_type::_outscale * x;
                sketch_of_A.Set(i, j, x);
            }
    }

    void apply_impl(const matrix_type& A,
        output_matrix_type& sketch_of_A,
        skylark::sketch::rowwise_tag tag) const {

        // TODO verify sizes etc.
        underlying_t underlying(*data_type::_underlying_data);
        underlying.apply(A, sketch_of_A, tag);

#       if SKYLARK_HAVE_OPENMP
#       pragma omp parallel for collapse(2)
#       endif
        for(int j = 0; j < data_type::_S; j++)
            for(int i = 0; i < base::Height(A); i++) {
                value_type x = sketch_of_A.Get(i, j);
                x += data_type::_shifts[j];

#               ifdef SKYLARK_EXACT_COSINE
                x = std::cos(x);
#               else
                // x = std::cos(x) is slow
                // Instead use low-accuracy approximation
                if (x < -3.14159265) x += 6.28318531;
                else if (x >  3.14159265) x -= 6.28318531;
                x += 1.57079632;
                if (x >  3.14159265)
                    x -= 6.28318531;
                x = (x < 0) ?
                    1.27323954 * x + 0.405284735 * x * x :
                    1.27323954 * x - 0.405284735 * x * x;
#               endif

                x = data_type::_outscale * x;
                sketch_of_A.Set(i, j, x);
            }
    }
};

template <typename ValueType,
          elem::Distribution ColDist,
          template <typename> class KernelDistribution>
struct RFT_t <
    elem::DistMatrix<ValueType, ColDist, elem::STAR>,
    elem::DistMatrix<ValueType, ColDist, elem::STAR>,
    KernelDistribution> :
        public RFT_data_t<KernelDistribution> {
    // Typedef value, matrix, transform, distribution and transform data types
    // so that we can use them regularly and consistently.
    typedef ValueType value_type;
    typedef elem::DistMatrix<value_type, ColDist, elem::STAR> matrix_type;
    typedef elem::DistMatrix<value_type,
                             ColDist, elem::STAR> output_matrix_type;
    typedef RFT_data_t<KernelDistribution> data_type;
private:
    typedef skylark::sketch::dense_transform_t <matrix_type,
                                                output_matrix_type,
                                                KernelDistribution>
    underlying_t;

protected:

    // Allow only creation by subclasses.

    RFT_t (int N, int S, base::context_t& context)
        : data_type (N, S, context) {

    }

public:

    RFT_t(const RFT_t<matrix_type,
                      output_matrix_type,
                      KernelDistribution>& other)
        : data_type(other) {

    }

    RFT_t(const data_type& other_data)
        : data_type(other_data) {

    }

    template <typename Dimension>
    void apply (const matrix_type& A,
                output_matrix_type& sketch_of_A,
                Dimension dimension) const {

        switch(ColDist) {
        case elem::VR:
        case elem::VC:
            try {
            apply_impl_vdist (A, sketch_of_A, dimension);
            } catch (std::logic_error e) {
                SKYLARK_THROW_EXCEPTION (
                    base::elemental_exception()
                        << base::error_msg(e.what()) );
            } catch(boost::mpi::exception e) {
                SKYLARK_THROW_EXCEPTION (
                    base::mpi_exception()
                        << base::error_msg(e.what()) );
            }

            break;

        default:
            SKYLARK_THROW_EXCEPTION (
                base::unsupported_matrix_distribution() );
        }
    }

private:
    void apply_impl_vdist (const matrix_type& A,
                     output_matrix_type& sketch_of_A,
                     skylark::sketch::columnwise_tag tag) const {
        underlying_t underlying(*data_type::_underlying_data);
        underlying.apply(A, sketch_of_A, tag);
        elem::Matrix<value_type> &Al = sketch_of_A.Matrix();
        for(int j = 0; j < base::Width(Al); j++)
            for(int i = 0; i < data_type::_S; i++) {
                value_type val = Al.Get(i, j);
                value_type trans =
                    data_type::_outscale * std::cos(val + data_type::_shifts[i]);
                Al.Set(i, j, trans);
            }
    }

    void apply_impl_vdist(const matrix_type& A,
        output_matrix_type& sketch_of_A,
        skylark::sketch::rowwise_tag tag) const {

        // TODO verify sizes etc.
        underlying_t underlying(*data_type::_underlying_data);
        underlying.apply(A, sketch_of_A, tag);
        elem::Matrix<value_type> &Al = sketch_of_A.Matrix();
        for(int j = 0; j < data_type::_S; j++)
            for(int i = 0; i < base::Height(Al); i++) {
                value_type val = Al.Get(i, j);
                value_type trans =
                    data_type::_outscale * std::cos(val + data_type::_shifts[j]);
                Al.Set(i, j, trans);
            }
    }

};

} } 
#endif // SKYLARK_RFT_ELEMENTAL_HPP