CUDA Lab Course Reference Manual 2011

/Users/sck/cuda/Praktikum_2011/include/lac/cublas_Vector.h

Go to the documentation of this file.

#ifndef cublas_Vector_H
#define cublas_Vector_H

#include <lac/expression_template.h>

// Vorwaertsdeklaration fuer expression templates
struct vmu;

#include <lac/cublas_Array.h>

template<typename>           class FullMatrixAccessor;

template<typename T> struct PrecisionTraits;

template <>
struct PrecisionTraits<cuComplex> {

    typedef float     NumberType;

    typedef cuComplex CudaComplex;
};

template <>
struct PrecisionTraits<cuDoubleComplex> {

    typedef double NumberType;

    typedef cuDoubleComplex CudaComplex;
};


template <>
struct PrecisionTraits<float> {

    typedef float     NumberType;

    typedef cuComplex CudaComplex;
};


template <>
struct PrecisionTraits<double> {

    typedef double          NumberType;

    typedef cuDoubleComplex CudaComplex;
};



template <typename T>
struct PrecisionTraits {

    typedef T NumberType;

};

template<typename T> struct One
{
    typedef T Type;
public:
    T operator()(bool plus=true);
};


template<>
inline cuComplex One<cuComplex>::operator ()(bool plus)
{
    Type result; result.x = (plus ? 1. : -1); result.y = 0.; return result;
}

template<>
inline cuDoubleComplex One<cuDoubleComplex>::operator ()(bool plus)
{
    Type result; result.x = (plus ? 1. : -1); result.y = 0.; return result;
}

template<typename T>
inline T One<T>::operator ()(bool plus) {  return (plus ? 1. : -1); }

// @sect3{Operator << for cuComplex Types}
//
std::ostream& operator << (std::ostream& out, const cuComplex& c)
{
    out << "(" << c.x << ", " << c.y << ")";
    return out;
}


std::ostream& operator << (std::ostream& out, const cuDoubleComplex& c)
{
    out << "(" << c.x << ", " << c.y << ")";
    return out;
}


namespace bw_types {

    template<typename, typename> class Array;
    template<typename, typename> class Matrix;

    template<typename, typename> class VectorView;
    template<typename, typename> class ColVectorView;

    template<typename, typename> class SubMatrixView;
}

namespace bw_types {



// @sect3{Klasse: Vector}
//
// Mit irgendwas muss man ja rechnen.
template<typename T, typename BW>
class Vector : public ::Subscriptor, protected Array<T, BW> {



    friend class bw_types::Matrix<T, BW>;

    friend class bw_types::SubMatrixView<T, BW>;

    template<typename, typename> friend class VectorView;

    template<typename, typename> friend class ColVectorView;


    Vector(const Vector<T, BW> & ) {}

public:

    typedef BW blas_wrapper_type;

        // Fuer Optimierungen der Matrix-Vektor-Operationen zwischen
        // (Teil-)Matrizen und (Teil-)Vektoren muss bekannt sein, was fuer ein
        // Typ die jeweiligen Objekte sind. Da statische Konstanten zur Compile-Zeit ausgewertet werden
        // koennen dadurch bedingte Ausdruecke fuer die Komilierung geschrieben werden.
    // Siehe beispielsweise SubMatrixView::vmult().
    static const bool is_vector_view = false;



    Vector();

    Vector(int n_elements);

    Vector(int n_elements, const Array<T, BW> & raw_data);

    Vector(const FullMatrixAccessor<T> & src,
                 int r_begin, int c);

    Vector(const Matrix<T, BW> & src,
                 int r_begin, int c);

    template<typename M, typename Op>
    Vector(const X_read_read<M, Op, Vector<T, BW> > & Ax);

    template<typename M,  typename Op,
    typename T_src>
    Vector(const // X_read_read<M, Op, Vector<T, BW> >
                                    X_read_read<M, Op, bw_types::ColVectorView<T, T_src> >
                                    & Ax);

        // Lese- und Schreibzugriff auf das unterliegende array - falls man es doch mal braucht.
        // Der Zwang explizit diese Funktion aufzurufen sollte Schutz genug sein gegen Missbrauch,
        // d.h., ist man der Meinung diese Funktion nehmen zu muessen, hat man im Allgemeinen was
        // falsch gemacht und sich nicht an die Philosophie der Objekt-Orientierung gehalten..
    inline Array<T, BW> & array() { return *this; }

    inline const Array<T, BW> & array() const { return *this; }

    Vector<T, BW> & operator = (const Vector<T, BW> & other);
    Vector<T, BW> & operator = (const std::vector<T> & other);


    template<typename T2>
    void push_to(std::vector<T2> & dst) const;

    Vector<T, BW> & operator = (const T value);

    template<typename T_src>
    Vector<T, BW> & operator = (const VectorView<T, T_src > & other);

    template<typename M, typename Op>
    Vector<T, BW> & operator = (const X_read_read<M, Op, Vector<T, BW> > & Ax);

    template<typename M, typename T_src>
    Vector<T, BW> & operator =(const
                                    X_read_read<M, vmu, ColVectorView<T, T_src> >
                                    & Ax);


    Vector<T, BW> & operator += (const Vector<T, BW> & other);

    template<typename T_src>
    Vector<T, BW> & operator += (const VectorView<T, T_src > & other);

    Vector<T, BW> & operator -= (const Vector<T, BW> & other);

    template<typename T_src>
    Vector<T, BW> & operator -= (const VectorView<T, T_src > & other);

    Vector<T, BW> & operator *= (const T scale);

    T operator * (const Vector<T, BW> & other);

    Vector<T, BW> & operator /= (const T scale);

    template<typename VECTOR>
    T dot(const VECTOR & other) const;

    void sadd (T alpha, const Vector<T, BW> & other);

    void reinit(int new_size);

    void print() const;

    T l2_norm() const;

    typename PrecisionTraits<T>::NumberType sum() const;

    int size()   const { return this->__n; }

    int n_rows() const { return this->__n; }

    int n_cols() const { return this->__n; }

    T operator () (int k) const;

    void set(int k,const T value);

    void add(int k,const T value);

protected:
    // A Vector is just the hull for adding array math operations to
    // a 1D array. Therefore, the stride is always 1.
   static const int _stride = 1;

};
}

// @sect3{bw_types::Vector Methoden}


// @sect4{Konstruktor: Vector}
//
// Standardkontruktor fuer Vector

template<typename T, typename BW>
bw_types::Vector<T, BW>::Vector()
    :
    Array<T, BW>()
{}

// @sect4{Konstruktor: Vector(n)}
//
// Konstruktor fuer n-elementigen Vector
// @param n_elements : Groesse des Vectors

template<typename T, typename BW>
bw_types::Vector<T, BW>::Vector(int n_elements)
    :
    Array<T, BW>(n_elements)
{
    this->reinit(n_elements);
}

// @sect4{Konstruktor: Vector(n,raw_data)}
//
// Konstruktor fuer n-elementigen Vector, dessen Inhalt aus einem Array-Objekt kommen
// @param n_elements : Groesse des Vectors
// @param raw_data : Pointer auf ein Array Objekt

template<typename T, typename BW>
bw_types::Vector<T, BW>::Vector(int n_elements,
                              const Array<T, BW> & raw_data)
    :
    Array<T, BW>(n_elements)
{
    this->reinit(n_elements);
    *this = raw_data;
}


// @sect4{Konstruktor: Vector(FullMatrixAccessor,r_begin,c)}
//
// Konstruiert einen Spalten-Vector aus einer Matrix, indem das Startelement angegeben wird
// @param src : Pointer auf einen FullMatrixAccessor
// @param r_begin : Zeilenindex des der Matrix
// @param c : Spaltenindex der Matrix

template<typename T, typename BW>
bw_types::Vector<T, BW>::Vector(const FullMatrixAccessor<T> & src,
                              int r_begin, int c)
    :
    Array<T, BW>(src.n_rows() - r_begin)
{

    int n_el = src.n_rows() - r_begin;

    this->reinit(n_el);

    int src_begin = c*src.n_rows() + r_begin;
    const T * src_val = src.val();

    int inc_src = 1;
    if (!src.is_column_major())
    {
        inc_src = src.n_cols();
        src_begin = r_begin*inc_src + c;
    }

    const T * src_ptr = &(src_val[src_begin]);

    BW::SetVector(n_el, src_ptr, inc_src, this->val(), 1);
}


// @sect4{Konstruktor: Vector(Matrix,r_begin,c)}
//
// Konstruiert einen Spalten-Vector aus einer Matrix, indem das Startelement angegeben wird
// @param src : Pointer auf eine Matrix
// @param r_begin : Zeilenindex des der Matrix
// @param c : Spaltenindex der Matrix

template<typename T, typename BW>
bw_types::Vector<T, BW>::Vector(const Matrix<T, BW> & src,
                              int r_begin, int c)
    :
    Array<T, BW>(src.n_rows() - r_begin)
{
    int n_el = src.n_rows() - r_begin;

    this->reinit(n_el);

        // Elementweise Kopie des arrays.
    int inc_src = 1;
    int inc_this = 1;

    const T * col = &(src.val()[c*src.n_rows() + r_begin]);

    // cublasScopy
    BW::copy(this->__n, col, inc_src, this->val(), inc_this);

}



// @sect4{Funktion: reinit}
//
// Erzeugt Nullvektor der Lange @p new_size auf der GPU.
// @param newsize : Groesse des Vectors

template<typename T, typename BW>
void bw_types::Vector<T, BW>::reinit(int new_size)
{
    this->Array<T, BW>::reinit(new_size);

    // *this *= 0.;
}



// @sect4{Funktion: l2_norm}
//
// Berechnet die L2-Norm eines Vectors

template<typename T, typename BW>
T bw_types::Vector<T, BW>::l2_norm() const
{
    return BW::nrm2(this->__n, &(this->val()[0]), 1);
}

// @sect4{Funktion: sum}
//
// Berechnet die Summe über alle Vektorelemente

template<typename T, typename BW>
typename PrecisionTraits<T>::NumberType bw_types::Vector<T, BW>::sum() const
{
    return BW::asum(this->__n, &(this->val()[0]), 1);
}


// @sect4{Funktion: Vector::print}
//
// Kopiert die Daten vom device zurueck zum host
// und gibt sie auf der Standardausgabe aus.

template<typename T, typename BW>
void
bw_types::Vector<T, BW>::print() const
{

    std::vector<T> tmp(this->__n);
    T * dst_ptr = &tmp[0];

    const T * src_ptr = this->val();

    BW::GetVector(this->__n, src_ptr, 1, dst_ptr, 1);

    for (int i = 0; i < this->__n; ++i)
        std::cout  << tmp[i] << std::endl;
}

// @sect4{Funktion: Vector::()}
//
// Elementzugriff auf Vektor - lesend
//
// @param k : Index des Vectors

template<typename T, typename BW>
T
bw_types::Vector<T, BW>::operator () (int k) const
{

    std::vector<T> tmp(this->__n);
    T * dst_ptr = &tmp[0];

    BW::GetVector(1, this->val()+k, 1, dst_ptr, 1);
    return tmp[0];
}

// @sect4{Funktion: Vector::set}
//
// Elementzugriff auf Vektor - schreibend
// @param k : Index des Vectors
// @param value : zu schreibender Wert
template<typename T, typename BW>
void
bw_types::Vector<T, BW>::set(int k,const T value)
{
    int inc_src = 1;
    int inc_dst = 1;

    BW::SetVector(1, &value, inc_src, &(this->val()[k]), inc_dst);
}

// @sect4{Funktion: Vector::add}
//
// Elementzugriff auf Vektor - schreibend
// @param k : Index des Vectors
// @param value : zu addierender Wert
template<typename T, typename BW>
void
bw_types::Vector<T, BW>::add(int k,const T value)
{
    Assert((k >= 0) && (k < this->size()),
           ::ExcMessage("Index out of range") );


    Vector<T, BW> tmp_d(1);
    tmp_d.set(0,value);


    BW::axpy(1/*other.size()*/, 1, tmp_d.val()/*other.val()*/, 1,
             &(this->val()[k]), 1);
}



// @sect4{Operator: =}
//
// Elementweise Kopie eines Vektors. Das Ziel wird bei Bedarf
// in der Groesse der Quelle angepassst.
// @param other : rechte Seite des = ist ein Vector
template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator = (const Vector<T, BW> & other)
{

    if(this->__n != other.__n) this->reinit(other.__n);

    // Elementweise Kopie des arrays.
    int inc_src = 1;
    int inc_this = 1;

    BW::copy(this->__n, other.val(), inc_src, this->val(), inc_this);

    return *this;
}

// @sect4{Operator: =}
//
// Elementweise Kopie eines Vektors. Das Ziel wird bei Bedarf
// in der Groesse der Quelle angepassst.
// @param other : rechte Seite des = ist ein std::Vector

template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator = (const std::vector<T> & other)
{
    this->reinit(other.size());

    const T * tmp_ptr = &other[0];

    BW::SetVector(size(), tmp_ptr, 1, this->data(), 1);

    return *this;
}

template<typename T, typename BW>
template<typename T2>
void
bw_types::Vector<T, BW>::push_to(std::vector<T2> & dst) const
{
    dst.resize(this->size());

    const T * const src_ptr = this->val();

    T * dst_ptr = &dst[0];

    static int inc_src = 1;

    static int inc_dst = 1;

    BW::GetVector(this->size(), src_ptr, inc_src, dst_ptr, inc_dst);
}



// @sect4{Operator: =}
//
// Elementweise Kopie eines Vektors. Das Ziel wird bei Bedarf
// in der Groesse der Quelle angepassst.
// @param other : rechte Seite des = ist ein VectorView
// Beim Kopieren wird das Koordinatensystem des Views uebernommen,
// d.h. Beginn und Ende des Views bestimmen den Bereich in den kopiert wird.
// Siehe Berechnung der Householder-Vektoren im QR-Beispiel.
template<typename T, typename BW>
template<typename T_src>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator = (const VectorView<T, T_src > & other)
{
    Assert(this->size() >= other.size(),
           ::ExcMessage("Dimension mismatch") );

    int incx = other._stride;
    int incy = 1;
    BW::copy(other.size(), other.val(), incx,
             &(this->val()[(other._is_col ?
                              other.r_begin() : other.c_begin())
                             ]), incy);

    return *this;
}



// @sect4{Operator: =}
//
template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator = (const T value)
{
    Vector<T, BW> tmp(1);
    tmp.set(0, value);
    int incx = 0; // copy the same element all the time
    int incy = 1;
    BW::copy(tmp.size(), tmp.val(), incx, this->val(), incy);

    return *this;
}



// @sect4{Operator: +=}
//
// Addiert zwei Vectoren
// @param other : rechte Seite des += ist ein Vector

template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator += (const Vector<T, BW> & other)
{
    One<T> one;
    BW::axpy(this->__n, one(), other.val(), 1, this->val(), 1);

    return *this;
}

// @sect4{Operator: +=}
//
// Addiert VectorView zu Vector
// @param other : rechte Seite des += ist ein VectorView

template<typename T, typename BW>
template<typename T_src>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator += (const VectorView<T, T_src > & other)
{
    Assert(this->size() >= other.size(),
           ::ExcMessage("Dimension mismatch") );

    BW::axpy(other.size(), 1, other.val(), 1,
             &(this->val()[other.r_begin()]), 1);

    return *this;
}



// @sect4{Operator: *=}
//
// Skalare Vectormultiplikation
// @param scale : Faktor

template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator *= (const T scale)
{
    int elem_dist = 1;

    BW::scal(this->__n, scale, &(this->val()[0]), elem_dist);

    return *this;
}

// @sect4{Operator: -=}
//
// Subtrahiert einen Vector von einem anderen
// @param other : Vector, der subtrahiert werden soll


template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator -= (const Vector<T, BW> & other)
{
    One<T> one;
    bool plus = true;

    BW::axpy(this->__n, one(!plus), other.val(), 1, this->val(), 1 );

    return *this;
}


// @sect4{Operator: -=}
//
// Subtrahiert einen VectorView von einem Vector
// @param other : VectorView, der subtrahiert werden soll

template<typename T, typename BW>
template<typename T_src>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator -= (const VectorView<T, T_src > & other)
{
    Assert(this->size() >= other.size(),
           ::ExcMessage("Dimension mismatch") );

    BW::axpy(other.size(), -1, other.val(), 1,
             &(this->val()[other.r_begin()]), 1);

    return *this;
}



 // @sect4{Operator: *}
 //
// Skalarprodukt zweier Vectoren
// @param other : zweiter Vektor

template<typename T, typename BW>
T
bw_types::Vector<T, BW>::operator * (const Vector<T, BW> & other)
{
     Assert(this->size() == other.size(),
           ::ExcMessage("Dimension mismatch") );
    return this->dot(other);
}



// @sect4{Operator: /=}
//
// Skalare Vectordivision
// @param scale : Faktor

template<typename T, typename BW>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator /= (const T scale)
{
    int elem_dist = 1;

    Assert(scale,
           ::ExcMessage("Division by Zero") );

    BW::scal(this->__n, 1/scale, this->val(), elem_dist);

    return *this;
}


// @sect4{Funktion: dot}
//
// Berechnet das Skalarprodukt von zwei Vektoren.
// @param other : zweiter Vector

template<typename T, typename BW>
template<typename VECTOR>
T
bw_types::Vector<T, BW>::dot(const VECTOR & other) const
{

    //AssertThrow(false,
     //           ::ExcMessage("Strange implementation. "
     //                              "Please check source code of the body of this function. SCK"));

     Assert(this->size() == other.size(),
     ::ExcMessage("Dimension mismatch"));

    int incx = 1; // this->_stride;
    int incy = 1; //other._stride;
    T result = BW::dot(this->__n,
                       this->val(), incx, other.val(), incy);

    return result;
}


// @sect4{Funktion: sadd}
//
// \f{eqnarray*} this & = & \alpha \cdot this \\ this & = & this + other \f}
// @param alpha : Faktor
// @param other : zweiter Vector

template<typename T, typename BW>
void
bw_types::Vector<T, BW>::sadd (T alpha, const Vector<T, BW> & other)
{
    int incx = 1;
    BW::scal(this->__n, alpha, this->val(), incx);

    // cublasSaxpy
    int incy = 1;
    BW::axpy(this->__n, 1., other.dev_ptr, 1, this->val(), 1);
}



// @sect3{Expression Templates}

template<typename M>
struct transpose {

    const M & A;

    transpose(const M & m) : A(m) {}

};


// @sect4{struct: vmu}
//
// Struktur zum Ausfuehren einer Matrix Vector Multiplikation
// @param b : Zielvector
// @param A : Matrix
// @param x : Vector
struct vmu
{
    template<typename T, typename BW>
    static void apply(      bw_types::Vector<T, BW> & b,
                      const bw_types::Matrix<T, BW> & A,
                      const bw_types::Vector<T, BW> & x)
    {
        b.reinit(A.n_rows());   A.vmult(b,x);
    }

    template<typename T, typename BW>
    static void apply(      bw_types::Vector<T, BW> & b,
                      const bw_types::SubMatrixView<T, BW> & A,
                      const bw_types::Vector<T, BW> & x)
    {
        b.reinit(A.r_end() - A.r_begin());   A.vmult(b,x);
    }


    template<typename T, typename BW>
    static void apply(      bw_types::Vector<T, BW> & b,
                      const transpose<bw_types::Matrix<T, BW> > & A_t,
                      const bw_types::Vector<T, BW> & x)
    {
        b.reinit(A_t.A.n_cols());    A_t.A.Tvmult(b,x);
    }



    template<typename T, typename BW>
    static void apply(      bw_types::Vector<T, BW> & b,
                      const transpose<bw_types::SubMatrixView<T, BW> > & A_t,
                      const bw_types::Vector<T, BW> & x)
    {
        b.reinit(A_t.A.matrix().n_rows());   A_t.A.Tvmult(b,x);
    }


    template<typename T, typename BW, typename T_src>
    static void apply(      bw_types::Vector<T, BW> & b,
                      const bw_types::Matrix<T, BW> & A,
                      const bw_types::ColVectorView<T, T_src>  & x)
    {
        b.reinit(A.n_rows());   A.vmult(b,x);
    }
};

// @sect4{Operator *}
//
// Operator fuer Expression Template vmu
// @param _l : Referenz auf rechte Seite (Matrix)
// @param _r : Referenz auf linke Seite (Vector)
template<typename L, typename T, typename BW>
inline
X_read_read<L, vmu, bw_types::Vector<T, BW> >
operator * (const L & _l, const bw_types::Vector<T, BW> & _r)
{
    return X_read_read<L, vmu, bw_types::Vector<T, BW> > (_l,_r);
}

template<typename L, typename T, //typename R> //
typename BW, typename T_src>
inline
X_read_read<L, vmu, bw_types::ColVectorView<T, T_src>
>
operator * (const L & _l, const // R
             bw_types::ColVectorView<T, T_src>
            & _r)
{
    return X_read_read<L, vmu, // R
            bw_types::ColVectorView<T, T_src>
            > (_l,_r);
}




template<typename T, typename BW>
template<typename M, typename Op>
bw_types::Vector<T, BW>::Vector(const X_read_read<M, Op, Vector<T, BW> > & Ax)
{
    Ax.apply(*this);
}


template<typename T, typename BW>
template<typename M, typename Op
          , typename T_src
>
bw_types::Vector<T, BW>::Vector(const // X_read_read<M, Op, Vector<T, BW> >
                                X_read_read<M, Op, bw_types::ColVectorView<T, T_src> >
                                & Ax)
{
    Ax.apply(*this);
}




template<typename T, typename BW>
template<typename M, typename Op>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator = (const X_read_read<M, Op, Vector<T, BW> > & Ax)
{
    Ax.apply(*this);

    return *this;
}



template<typename T, typename BW>
template<typename M, // typename Op,
typename T_src>
bw_types::Vector<T, BW> &
bw_types::Vector<T, BW>::operator =(const
                                X_read_read<M, vmu, bw_types::ColVectorView<T, T_src> >
                                & Ax)
{
    Ax.apply(*this);

     return *this;
}


#endif // cublas_Vector_H

---