/data/development/ViennaCL/dev/viennacl/vandermonde_matrix.hpp

Go to the documentation of this file.

#ifndef VIENNACL_VANDERMONDE_MATRIX_HPP
#define VIENNACL_VANDERMONDE_MATRIX_HPP

/* =========================================================================
   Copyright (c) 2010-2011, Institute for Microelectronics,
                            Institute for Analysis and Scientific Computing,
                            TU Wien.

                            -----------------
                  ViennaCL - The Vienna Computing Library
                            -----------------

   Project Head:    Karl Rupp                   rupp@iue.tuwien.ac.at
               
   (A list of authors and contributors can be found in the PDF manual)

   License:         MIT (X11), see file LICENSE in the base directory
============================================================================= */

#include <cmath>

#include "viennacl/forwards.h"
#include "viennacl/vector.hpp"
#include "viennacl/ocl/context.hpp"

#include "viennacl/fft.hpp"

#include "viennacl/linalg/vandermonde_matrix_operations.hpp"

namespace viennacl {
    template<class SCALARTYPE, unsigned int ALIGNMENT>
    class vandermonde_matrix
    {
      public:
        explicit vandermonde_matrix()
        {
          viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
        }

        explicit vandermonde_matrix(std::size_t rows, std::size_t cols) : elements_(rows)
        {
          assert(rows == cols && "Vandermonde matrix must be square in this release!");
          viennacl::linalg::kernels::fft<SCALARTYPE, 1>::init();
        }

        void resize(std::size_t sz, bool preserve = true) {
            elements_.resize(sz, preserve);
        }

        viennacl::ocl::handle<cl_mem> handle() const { return elements_.handle(); }

        viennacl::vector<SCALARTYPE, ALIGNMENT> & elements() { return elements_; }
        viennacl::vector<SCALARTYPE, ALIGNMENT> const & elements() const { return elements_; }

        std::size_t size1() const { return elements_.size(); }
        
        std::size_t size2() const { return elements_.size(); }

        std::size_t internal_size() const { return elements_.internal_size(); }

        entry_proxy<SCALARTYPE> operator()(std::size_t row_index)
        {
            return elements_[row_index];
        }

        SCALARTYPE operator()(std::size_t row_index, std::size_t col_index) const
        {
            assert(row_index < size1() && col_index < size2() && "Invalid access");
            
            return pow(elements_[row_index], static_cast<int>(col_index));
        }

    private:
        vandermonde_matrix(vandermonde_matrix const & t) {}
        vandermonde_matrix & operator=(vandermonde_matrix const & t) {}
        
        viennacl::vector<SCALARTYPE, ALIGNMENT> elements_;
    };

    template <typename SCALARTYPE, unsigned int ALIGNMENT>
    void copy(std::vector<SCALARTYPE>& cpu_vec, vandermonde_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat)
    {
        assert(cpu_vec.size() == gpu_mat.size1()  && "Size mismatch");
        copy(cpu_vec, gpu_mat.elements());
    }

    template <typename SCALARTYPE, unsigned int ALIGNMENT>
    void copy(vandermonde_matrix<SCALARTYPE, ALIGNMENT>& gpu_mat, std::vector<SCALARTYPE>& cpu_vec)
    {
        assert(cpu_vec.size() == gpu_mat.size1() && "Size mismatch");
        copy(gpu_mat.elements(), cpu_vec);
    }

    template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
    void copy(vandermonde_matrix<SCALARTYPE, ALIGNMENT>& vander_src, MATRIXTYPE& com_dst)
    {
        std::size_t size = vander_src.size1();
        assert(size == com_dst.size1() && "Size mismatch");
        assert(size == com_dst.size2() && "Size mismatch");
        std::vector<SCALARTYPE> tmp(size);
        copy(vander_src, tmp);

        for(std::size_t i = 0; i < size; i++) {
            for(std::size_t j = 0; j < size; j++) {
                com_dst(i, j) = pow(tmp[i], static_cast<int>(j));
            }
        }
    }
    
    template <typename SCALARTYPE, unsigned int ALIGNMENT, typename MATRIXTYPE>
    void copy(MATRIXTYPE& com_src, vandermonde_matrix<SCALARTYPE, ALIGNMENT>& vander_dst) 
    {
        std::size_t size = vander_dst.size1();
        assert(size == com_src.size1() && "Size mismatch");
        assert(size == com_src.size2() && "Size mismatch");
        std::vector<SCALARTYPE> tmp(size);

        for(std::size_t i = 0; i < size; i++)
            tmp[i] = com_src(i, 1);

        copy(tmp, vander_dst);
    }

    /*template <typename SCALARTYPE, unsigned int ALIGNMENT, unsigned int VECTOR_ALIGNMENT>
    void prod_impl(vandermonde_matrix<SCALARTYPE, ALIGNMENT>& mat,
                   vector<SCALARTYPE, VECTOR_ALIGNMENT>& vec,
                   vector<SCALARTYPE, VECTOR_ALIGNMENT>& result) {
        assert(mat.size1() == vec.size());

        fft::vandermonde_prod<SCALARTYPE>(mat.handle(), vec.handle(), result.handle(), mat.size1());
    } */

    template<class SCALARTYPE, unsigned int ALIGNMENT>
    std::ostream & operator<<(std::ostream& s, vandermonde_matrix<SCALARTYPE, ALIGNMENT>& gpu_matrix)
    {
        std::size_t size = gpu_matrix.size1();
        std::vector<SCALARTYPE> tmp(size);
        copy(gpu_matrix, tmp);
        s << "[" << size << "," << size << "](\n";

        for(std::size_t i = 0; i < size; i++) {
            s << "(";
            for(std::size_t j = 0; j < size; j++) {
                s << pow(tmp[i], j);
                if(j < (size - 1)) s << ",";
            }
            s << ")";
        }
        s << ")";
        return s;
    }

}

#endif // _VIENNACL_VANDERMONDE_MATRIX_HPP