CUDA Lab Course Reference Manual 2011

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

Go to the documentation of this file.

#ifndef cublas_Matrix_H
#define cublas_Matrix_H


#include <iomanip>
#include <cmath>
#include <QtGlobal>

#include <base/subscriptor.h>

#include <lac/cublas_Array.h>
#include <lac/expression_template.h>
#include <lac/cublas_Vector.h>

struct mmu;

template<typename> class FullMatrixAccessor;


namespace bw_types {

    template<typename, typename> class Vector;

    template<typename, typename> class VectorView;

    template<typename, typename> class SubMatrixView;
}


namespace bw_types {

    // @sect3{Klasse: Matrix}
    //
    template<typename T, typename BW>
    class Matrix
        :
        protected bw_types::Array<T, BW>,
        public ::Subscriptor {

        friend class Vector<T, BW>;

        friend class VectorView<T, Matrix<T, BW> >;

        friend class SubMatrixView<T, BW>;

        friend class FullMatrixAccessor<T>;

    public:

        typedef T Number;

        typedef T value_type;

        typedef BW blas_wrapper_type;

        Matrix();

        Matrix(int n_rows, int n_cols);

        Matrix(int n_rows, int n_cols, const Matrix<T, BW> & src_data);

        Matrix(const Matrix<T, BW> & other);

        template<typename L, typename R>
        Matrix(const X_read_read<L, mmu, R> & AB);


        Matrix(const ::IdentityMatrix & Id);

        void reinit(int n_rows, int n_cols);

        Matrix<T, BW> & operator = (const FullMatrixAccessor<T> & matrix);

        Matrix<T, BW> & operator = (const ::IdentityMatrix & Id);

        Matrix<T, BW> & operator = (const Matrix<T, BW> & array);


        template<typename L, typename R>
        Matrix<T, BW> & operator = (const X_read_read<L, mmu, R> & AB);


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

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


        template<typename VECTOR1, typename VECTOR2>
        void vmult(VECTOR1& dst, const VECTOR2& src) const;

        template<typename VECTOR1, typename VECTOR2>
        void Tvmult(VECTOR1& dst, const VECTOR2& src) const;

        void scaled_vmult(T beta, Vector<T, BW>& dst,
                          T alpha, const Vector<T, BW>& src) const;

        // @f$dst = this \cdot src@f$
        void mmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const;

        void mmult(SubMatrixView<T, BW>& dst, const Matrix<T, BW>& src) const;

        // @f$dst = this \cdot src^T@f$
        void mTmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const;

        // @f$dst = this^T \cdot src@f$
        void Tmmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const;

        // @f$dst = this^T \cdot src^T@f$
        void TmTmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const;

        template<typename VECTOR1, typename VECTOR2>
        void add_scaled_outer_product(T alpha, const VECTOR1& col, const VECTOR2& row);



        int n_rows() const { return __n_rows; }

        int n_cols() const { return __n_cols; }

        void print() const;

        T operator () (const unsigned int i, const unsigned int j) const;

        T l2_norm() const;

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

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

    private:

        Matrix<T, BW> & operator = (const Array<T, BW> & src);

        int __n_rows;
        int __n_cols;
    };

}

// @sect3{bw_types::Matrix Methoden}
//
// @sect4{Konstruktoren}
//
template<typename T, typename BW>
bw_types::Matrix<T, BW>::Matrix()
    :
    Array<T, BW>(), __n_rows(0), __n_cols(0)
{}


// Initialisiert GPU-seitig die Matrix mit NUllen.
template<typename T, typename BW>
bw_types::Matrix<T, BW>::Matrix(int n_rows, int n_cols)
    :
    Array<T, BW>(n_rows*n_cols),
    __n_rows(n_rows),
    __n_cols(n_cols)
{
    const std::vector<T> tmp(n_rows*n_cols, 0);

    const T * const tmp_ptr = &tmp[0];

    T * this_data = this->data();

    BW::SetMatrix(n_rows, n_cols, tmp_ptr, n_rows,
                  this_data, n_rows);
}


template<typename T, typename BW>
bw_types::Matrix<T, BW>::Matrix(int n_rows, int n_cols,
                                const Matrix<T, BW> & src_data)
                                    :
                                    Array<T, BW>(n_rows*n_cols),
                                    __n_rows(n_rows),
                                    __n_cols(n_cols)
{
    const std::vector<T> tmp(n_rows*n_cols, 0);

    BW::SetMatrix(n_rows, n_cols, &tmp[0], n_rows,
                  this->data(), n_rows);

    Array<T, BW> & self = *this;

    self = src_data;
}


template<typename T, typename BW>
bw_types::Matrix<T, BW>::Matrix(const ::IdentityMatrix & Id)
    :
    Array<T, BW>(), __n_rows(0), __n_cols(0)
{
    *this = Id;
}




template<typename T, typename BW>
bw_types::Matrix<T, BW>::Matrix(const Matrix<T, BW> & other)
{
    *this = other;
}



template<typename T, typename BW>
void bw_types::Matrix<T, BW>::reinit(int n_rows, int n_cols)
{

    this->Array<T, BW>::reinit(n_rows*n_cols);
    __n_rows = n_rows;
    __n_cols = n_cols;

    /*
    std::vector<T> tmp(n_rows*n_cols, 0);

    BW::SetMatrix(n_rows, n_cols, &tmp[0], n_rows,
                  this->data(), n_rows);
    */
}


// @sect4{Operator: =}
//
// Elementweise Kopie des arrays.
// Die Quelle muss mindestens so viele Elemente wie das Ziel haben.
template<typename T, typename BW>
bw_types::Matrix<T, BW> &
        bw_types::Matrix<T, BW>::operator = (const Array<T, BW> & src)
                                            {
    Assert(this->n_elements() <= src.n_elements(),
           ::ExcMessage("n_element mismatch") );

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

    BW::copy(this->n_elements(), src.data(), inc_src,
             this->data(), inc_this);

    return *this;
}

template<typename T, typename BW>
bw_types::Matrix<T, BW> &
        bw_types::Matrix<T, BW>::operator = (const ::IdentityMatrix & Id)
                                            {
    this->__n_rows = Id.m();
    this->__n_cols = Id.n();

    // Zeilen- und Spaltenzahl einer Identitaetsmatrix sind gleich
    // und entsprechen der Anzahl der Freiheitsgrade.
    int n_dofs   = Id.m();
    int n_src_el = n_dofs*n_dofs;

    // Passe bei Bedarf die Array-Groesse an.
    this->Array<T,BW>::reinit(n_src_el);

    FullMatrixAccessor<T> tmp_id(Id);

    // Laut CUBLAS Doku wird die 'leading dimension' immer durch die Anzahl
    // der Zeilen angegeben.
    // Ausserdem muss man sich bei Einheitsmatrizen nicht darum kuemmern,
    // ob sie column-major oder row-major sind.
    const T *  id_val = tmp_id.val();
    T * dst_val = this->data();
    BW::SetMatrix(n_dofs, n_dofs, id_val,
                  n_dofs, dst_val, n_dofs );

    return *this;
}



template<typename T, typename BW>
bw_types::Matrix<T, BW> &
        bw_types::Matrix<T, BW>::operator = (const FullMatrixAccessor<T> & src_matrix)
                                            {


    Assert(src_matrix.is_column_major(),
           ::ExcMessage("bw_types:Matrix expects a matrix in column major"
                              " format as source.") );

    int n_src_el = src_matrix.n_elements();

    this->Array<T, BW>::reinit(n_src_el);

    int nr = src_matrix.n_rows();
    int nc = src_matrix.n_cols();

    const T * tmp_src = src_matrix.val();
    // Laut CUBLAS Doku wird die 'leading dimension' immer durch die Anzahl
    // der Zeilen angegeben.
    T * tmp_dst = this->data();
    BW::SetMatrix(nr, nc, tmp_src, nr, tmp_dst, nr);


    this->__n_rows = nr;
    this->__n_cols = nc;

    return *this;
}


// Beim Kopieren wird das Ziel neu angelegt. Jeglicher frueherer Inhalt geht verloren.
// @param other : zu kopierende Matrix
template<typename T, typename BW>
bw_types::Matrix<T, BW> &
bw_types::Matrix<T, BW>::operator = (const Matrix<T, BW> & other)
{
    this->__n_rows = other.__n_rows;
    this->__n_cols = other.__n_cols;

    this->Array<T, BW>::reinit(__n_rows*__n_cols);

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

    BW::copy(this->n_elements(), other.data(), inc_src,
             this->data(), inc_this);

    return *this;
}

// @sect4{Operator: +=}
//
// Elementweise Differenz von zwei Matrizen.
// @param other : abzuziehende Matrix. Muss die gleichen Dimensionen haben.
template<typename T, typename BW>
bw_types::Matrix<T, BW> & bw_types::Matrix<T, BW>::operator += (const Matrix<T, BW> & other)
                                                               {

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

    int n = this->n_elements();
    T alpha = 1.;

    const T * const x = other.array().val();
    int incx = 1;

    T * y = this->val();
    int incy = 1;

    BW::axpy(n, alpha, x, incx, y, incy);

    return *this;
}


// @sect4{Operator: -=}
//
// Elementweise Differenz von zwei Matrizen.
// @param other : abzuziehende Matrix. Muss die gleichen Dimensionen haben.
template<typename T, typename BW>
bw_types::Matrix<T, BW> & bw_types::Matrix<T, BW>::operator -= (const Matrix<T, BW> & other)
                                                               {

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

    int n = this->n_elements();
    T alpha = -1.;

    const T * const x = other.array().val();
    int incx = 1;

    T * y = this->val();
    int incy = 1;

    BW::axpy(n, alpha, x, incx, y, incy);

    return *this;
}



// @sect4{Function: l2_norm}
//
template<typename T, typename BW>
T bw_types::Matrix<T, BW>::l2_norm() const
{

    T result = BW::nrm2(this->n_elements(), this->val(), 1/*incx*/);

    return result;

}


// @sect4{Function: vmult}
//
// Matrix-Vektor-Multiplikation.
// @f$dst = A \cdot src@f$
template<typename T, typename BW>
template<typename VECTOR1, typename VECTOR2>
inline
void bw_types::Matrix<T, BW>::vmult(VECTOR1& dst, const VECTOR2& src) const
{
    // y = alpha*A*x + beta*y
    T alpha = 1.;
    T beta = 0.;
    int leading_dim_A = this->__n_rows;

    T * dst_ptr = dst.val();
    const T * src_ptr = src.val();

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

    BW::gemv('n', this->__n_rows, this->__n_cols, alpha, this->data(),
             leading_dim_A, src_ptr, 1, beta, dst_ptr, 1);
}



// @sect4{Function: Tvmult}
//
//
// Matrix-Vektor-Multiplikation mit transponierter Matrix.
// @f$dst = A^T \cdot src@f$
template<typename T, typename BW>
template<typename VECTOR1, typename VECTOR2>
inline
void bw_types::Matrix<T, BW>::Tvmult(VECTOR1& dst, const VECTOR2& src) const
{
    // y = alpha*A*x + beta*y
    T alpha = 1.;
    T beta = 0.;
    int leading_dim_A = this->__n_rows;

    BW::gemv('t', this->__n_rows, this->__n_cols, alpha, this->data(),
             leading_dim_A, src.val(), 1, beta, dst.val(), 1);
}



// @sect4{Function: scaled_vmult}
//
// Allgemeine Matrix-Vektor-Multiplikation.
// @f$dst = \alpha A \cdot src + \beta dst@f$
template<typename T, typename BW>
void
        bw_types::Matrix<T, BW>::scaled_vmult(T beta,        Vector<T, BW>& dst /*y*/,
                                              T alpha, const Vector<T, BW>& src /*x*/)
        const
{
    BW::gemv('n', this->__n_rows, this->__n_cols, alpha, this->data(),
             this->__n_rows, src.val(), 1, beta, dst.val(), 1);
}




// @sect4{Function: mmult}
//
// matrix-matrix multiplication
// @f$dst = A \cdot src@f$
template<typename T, typename BW>
void
bw_types::Matrix<T, BW>::mmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const
{
    T alpha = 1;
    T beta = 0;

    int lda = this->__n_rows;
    int ldb = src.__n_rows; /* == this->n_cols !!! */
    int ldc = dst.__n_rows;

    BW::gemm('n', 'n',
             this->__n_rows /* cublas doc : m == n_rows of op(A)*/,
             dst.__n_cols   /* cublas doc : n == n_cols of op(B), i.e. n_cols of C */,
             this->__n_cols /* cublas doc : k == n_cols of op(A), i.e. n_rows of op(B)*/,
             alpha,
             this->data(), lda,
             src.data(), ldb,
             beta,
             dst.data(), ldc);
}


template<typename T, typename BW>
void
bw_types::Matrix<T, BW>::mmult(SubMatrixView<T, BW>& dst, const Matrix<T, BW>& src) const
{
    T alpha = 1;
    T beta = 0;


    int lda = this->__n_rows;
    int ldb = src.__n_rows; /* == this->n_cols !!! */
    int ldc = dst.__n_rows;

    BW::gemm('n', 'n',
             this->__n_rows /* cublas doc : m == n_rows of op(A)*/,
             dst.__n_cols   /* cublas doc : n == n_cols of op(B), i.e. n_cols of C */,
             this->__n_cols /* cublas doc : k == n_cols of op(A), i.e. n_rows of op(B)*/,
             alpha,
             this->data(), lda,
             src.data(), ldb,
             beta,
             dst.data(), ldc);
}



// @sect4{Function: mTmult}
//
// Matrix-Matrix-Multiplikation
// @f$dst = this \cdot src^T@f$
template<typename T, typename BW>
void
bw_types::Matrix<T, BW>::mTmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const
{
    T alpha = 1;
    T beta = 0;

    int lda = this->__n_rows; // this == A
    int ldb = src.__n_rows; /* == this->n_cols !!! */ // src == B
    int ldc = dst.__n_rows; // dst == C

    BW::gemm('n', 't',
             this->__n_rows /* cublas doc : m == n_rows of op(A), i.e. n_rows for A */,
             dst.__n_cols /* cublas doc : n == n_cols of op(B), i.e. n_cols of C */,
             this->__n_cols /* cublas doc : k == n_cols of op(A), i.e. n_rows of op(B) or n_cols for A */,
             alpha,
             this->data(), lda,
             src.data(), ldb,
             beta,
             dst.data(), ldc);
}


// @sect4{Function: Tmmult}
//
// Matrix-Matrix-Multiplikation.
// @f$dst = this^T \cdot src@f$
template<typename T, typename BW>
void
        bw_types::Matrix<T, BW>::Tmmult(Matrix<T, BW>& dst, const Matrix<T, BW>& src) const
{
    T alpha = 1;
    T beta = 0;

    int lda = this->__n_rows; // this == A
    int ldb = src.__n_rows; /* == this->n_cols !!! */ // src == B
    int ldc = dst.__n_rows; // dst == C

    BW::gemm('t', 'n',
             this->__n_cols /* cublas doc : m == n_rows of op(A), i.e. n_cols for A^T*/,
             dst.__n_cols /* cublas doc : n == n_cols of op(B), i.e. n_cols of C */,
             this->__n_rows /* cublas doc : k == n_cols of op(A), i.e. n_rows of op(B) or n_rows for A^T */,
             alpha,
             this->data(), lda,
             src.data(), ldb,
             beta,
             dst.data(), ldc);
}


// @sect4{Function: TmTmult}
//
// Matrix-Matrix-Multiplikation.
// @f$dst = this^T \cdot src^T@f$
template<typename T, typename BW>
void
        bw_types::Matrix<T, BW>::TmTmult(Matrix<T, BW>& , const Matrix<T, BW>& ) const
{
    AssertThrow(false, ::ExcNotImplemented() );
}



// @sect4{Function: add_scaled_outer_product}
//
// Wird fuer Householder-Transformationen gebraucht.
// @f$ this += \alpha x \cdot y^T @f$
template<typename T, typename BW>
template<typename VECTOR1, typename VECTOR2>
void
bw_types::Matrix<T, BW>::add_scaled_outer_product(T alpha,
                                                  const VECTOR1& x,
                                                  const VECTOR2& y)
{

    int m = this->__n_rows;
    int n = this->__n_cols;
    int lda = this->__n_rows;

#ifdef DEBUG
    if (x.size() != m) std::cout << "Dimension mismatch. "
            "Vector x.size() should be " << m << " but is " << x.size() << std::endl;
    if (y.size() != n) std::cout << "Dimension mismatch. "
            "Vector y.size() should be " << n << " but is " << y.size() << std::endl;
#endif
    Assert(x.size() == m, ::ExcMessage("Dimension mismatch"));
    Assert(y.size() == n, ::ExcMessage("Dimension mismatch"));

    int incx = x._stride;
    int incy = y._stride;

    BW::ger(m, n, alpha,
            x.val(), incx,
            y.val(), incy,
            this->data(), lda);
}

// @sect4{Operator: ()}
//
// Die elementweisen Zugriffe
// @param i : Zeilenindex in C-Zaehlung, d.h. bei 0 beginnend.
// @param j : Spaltenindex in C-Zaehlung, d.h. bei 0 beginnend.
template <typename T, typename BW>
        inline
        T
        bw_types::Matrix<T, BW>::operator () (const unsigned int r,
                                              const unsigned int c) const
{
    int lead_dim = this->n_rows();
    const T * tmp_d =  & this->data()[c*this->__n_rows+r];
    T entry;
    T * p_e = &entry;
    BW::GetMatrix(1, 1,
                  tmp_d, lead_dim, p_e, 1);

    return entry;
}


// @sect4{Function: Matrix::print}
//
// Kopiert die Daten vom device zurueck zum host
// und gibt sie auf der Standardausgabe aus.
template<typename T, typename BW>
void
        bw_types::Matrix<T, BW>::print() const
{
    // T numerical_zero = 2.5e-16;

    std::cout << "Matrix dims : " << this->__n_rows << " "
            << this->__n_cols << std::endl;

    int n_el = this->__n_rows * this->__n_cols;
    T * tmp = new T[n_el];

    BW::GetMatrix(this->__n_rows, this->__n_cols,
                  this->data(), this->__n_rows, tmp, this->__n_rows);

    for (int r = 0; r < this->__n_rows; ++r)
    {
        for (int c = 0; c < this->__n_cols; ++c)
            std::cout << std::setprecision(4) << std::fixed
                    << std::setw(15) <<
                    //(std::abs(tmp[c*this->__n_rows + r])> numerical_zero
                    //           ?
                    tmp[c*this->__n_rows + r]
                    //                 : 0.)
                    << " ";
        std::cout <<";" << std::endl;
    }

    delete [] tmp;

}




// @sect4{struct: mmu}
//
// Structure implementing matrix-matrix multiplication.
// @param C : destination
// @param A : left operand
// @param B : right operand
struct mmu
{
    template<typename T, typename BW>
    static void apply(      bw_types::Matrix<T, BW> & C,
                      const bw_types::Matrix<T, BW> & A,
                      const bw_types::Matrix<T, BW> & B)
    {
        C.reinit(A.n_rows(), B.n_cols());   A.mmult(C,B);
    }

    template<typename T, typename BW>
    static void apply(      bw_types::Matrix<T, BW> & C,
                      const bw_types::SubMatrixView<T, BW> & A,
                      const bw_types::Matrix<T, BW> & B)
    {
        AssertThrow(false, ::ExcMessage("Not implemented"));
       // C.reinit(A.n_rows(), B.n_cols());   A.mmult(C,B);
    }


    template<typename T, typename BW>
    static void apply(      bw_types::Matrix<T, BW> & C,
                      const bw_types::Matrix<T, BW> & A,
                      const transpose<bw_types::Matrix<T, BW> > & B_t)
    {
        C.reinit(A.n_rows(), B_t.A.n_rows());   A.mTmult(C, B_t.A);
    }


    template<typename T, typename BW>
    static void apply(      bw_types::Matrix<T, BW> & C,
                      const transpose<bw_types::Matrix<T, BW> > & A_t,
                      const bw_types::Matrix<T, BW>  & B)
    {
        C.reinit(A_t.A.n_cols(), B.n_cols());   A_t.A.Tmmult(C, B);
    }


    template<typename T, typename BW>
    static void apply(      bw_types::Matrix<T, BW> & C,
                      const transpose<bw_types::Matrix<T, BW> > & A_t,
                      const transpose<bw_types::Matrix<T, BW> > & B_t)
    {
        C.reinit(A_t.A.n_cols(), B_t.A.n_rows());   A_t.A.TmTmult(C, B_t.A);
    }

};

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

template<typename L, typename T, typename BW> // , typename R>
inline
X_read_read<L, mmu,  transpose<bw_types::Matrix<T, BW> > /*R*/>
operator * (const L & _l, const transpose<bw_types::Matrix<T, BW> > & _r)
{
    typedef transpose<bw_types::Matrix<T, BW> > R;
    return X_read_read<L, mmu, R> (_l,_r);
}



template<typename T, typename BW>
template<typename L, typename R>
bw_types::Matrix<T, BW>::Matrix(const X_read_read<L, mmu, R> & AB)
{
    AB.apply(*this);
}



template<typename T, typename BW>
template<typename L, typename R>
bw_types::Matrix<T, BW> &
bw_types::Matrix<T, BW>::operator = (const X_read_read<L, mmu, R > & AB)
{
    AB.apply(*this);
    return *this;
}




#endif

---