CUDA Lab Course Reference Manual 2011

/Users/sck/cuda/Praktikum_2011/include/lac/blas_wrapper.hh

Go to the documentation of this file.

#ifndef BLAS_WRAPPER_HH
#define BLAS_WRAPPER_HH

#ifdef __APPLE__
// On Apple systems BLAS is part of the Accelerate framework.
#include <Accelerate/Accelerate.h>
#else
#include <cblas.h>
#endif

#include <vector>

// @sect3{class CBW: cblas S,D,C,Z-Wrapper-Funktionen}
//
// Dasd Hauptproblem bei der Entwicklung eines modernen
// Interfaces fuer cublas ist die explizite Zahlentyp-Abhaengigkeit
// der einzelnen Funktionen. Heutzutage wuerde so etwas per template-Parameter
// erledigt.
// Um nicht fuer jede Klasse eine Spezialisierung fuer float, double
// und komplexe Zahlen in einfacher
// und doppelter Genauigkeit schreiben zu muessen, werden hier zuerst einige
// inline-Wrapper-Funktionen
// definiert, die den Zahlentyp als template-Parameter bekommen und fuer die
// die einzelnen
// Typen die entsprechende Spezialisierung anbieten.
//
// Namenskonvention: Aus dem cblas-Funktionsnamen das cblas_-Praefix
// und das Genauigkeitskuerzel weglassen Z.B.: cblas_sgemv wird zu gemv.
// Ansonsten siehe die CBLAS-Doku.
struct blas {

    static inline std::string name () { return "cpu-blas"; }

 #ifdef dfkgskdd

    // @sect4{Funktion: check_status}
    //
    // Wertet das cublasStatus Argument aus und wirft im Falle
    // eines Fehlers eine Exception aus.
    static void  check_status(const cublasStatus & status)
    {

        std::string cublas_errors(" ");

        if (status != CUBLAS_STATUS_SUCCESS)
        {
            if (status == CUBLAS_STATUS_NOT_INITIALIZED)
                cublas_errors += "cublas not initialized ";
            if (status == CUBLAS_STATUS_MAPPING_ERROR)
                cublas_errors +="mapping error ";
            if (status == CUBLAS_STATUS_INVALID_VALUE)
                cublas_errors +="invalid value ";
            if (status == CUBLAS_STATUS_ALLOC_FAILED)
                cublas_errors +="allocation failed ";
            if (status == CUBLAS_STATUS_ARCH_MISMATCH)
                cublas_errors +="architecture mismatch ";
            if (status == CUBLAS_STATUS_EXECUTION_FAILED)
                cublas_errors +="execution failed ";
            if (status == CUBLAS_STATUS_INTERNAL_ERROR)
                cublas_errors +="cublas internal error ";

            if (cublas_errors == " ")
                cublas_errors = "unknown cublas error state";

#ifdef QT_NO_DEBUG
       AssertThrow(false, ::ExcMessage(cublas_errors.c_str() ) );
#else
       AssertThrow(false, ::ExcMessage(cublas_errors.c_str() ) );
#endif
        }

    }
#endif

    // @sect4{Funktion: Init}
    //
    // Initialisiert cublas. Im Falle eines Fehlschlags
    // wird eine Exception ausgeworfen. Dazu wird auf die Asset und AssertThrow Makros
    // von deal.II zurueckgeworfen. Wenn man es ganz besonders genau machen wollte, wuerde
    // man spezielle Exceptions implementieren die man hinterher mit dem catch-Mechanismus von C++
    // abfangen koennte.
    static void  Init() {
/*
        cublasStatus s = cublasInit();

        if (s == CUBLAS_STATUS_SUCCESS)
            std::cout << "cublas init succeeded" << std::endl;

#ifdef QT_NO_DEBUG
        AssertThrow(s == CUBLAS_STATUS_SUCCESS,
                    ::ExcMessage("cublas init failed"));
#else
        Assert(s == CUBLAS_STATUS_SUCCESS,
               ::ExcMessage("cublas init failed"));
#endif
               */

        std::cout << "ATLAS Init" << std::endl;
    }


    // @sect4{Funktion: Shutdown}
    //
    // Faehrt cublas wieder herunter. Im Falle eines Fehlschlags
    // wird eine Exception ausgeworfen.
    static void Shutdown() {
/*
        cublasStatus s = cublasShutdown();

#ifdef QT_NO_DEBUG
        AssertThrow(s == CUBLAS_STATUS_SUCCESS,
                    ::ExcMessage("cublas shutdown failed"));
#else
        Assert(s == CUBLAS_STATUS_SUCCESS,
                    ::ExcMessage("cublas shutdown failed"));
#endif
*/
        std::cout << "ATLAS Down" << std::endl;

    }



    template<typename T>
    class Data {

        T * __data;
            size_t __n_el;
    public:

        Data() : __data(0)
        {}

        Data(size_t n)
           : __data(0), __n_el(0)
        { resize(n); }

        void resize(size_t n)
        {
#ifdef QT_NO_DEBUG
        AssertThrow(n > 0,
                    ::ExcMessage("allocation of 0 elements not allowed"));
#else
        Assert(n > 0,
                    ::ExcMessage("allocation of 0 elements not allowed"));
#endif

            if (__n_el == 0)
            {
                __data = new T[n];
                __n_el = n;
            }
            else {
            if (__n_el != n)
                delete __data;
                __data = new T[n];
                __n_el = n;
            }
        };

        ~Data()
        {
            if (__n_el > 0)
            {
                delete __data;
                __data = 0;
            }
        }

        T * data() { return __data; }

        const T * data() const { return __data; }
    };


public:
    // @sect4{Funktion: SetMatrix}
    //
    // @param rows : Anzahl Zeilen.
    // @param cols : Anzahl Spalten.
    // @param A : Quellmatrix A.
    // @param lda : leading dimension von A.
    // @param B : Zielmatrix B.
    // @param ldb : leading dimension von B.
    template<typename T>
    static void SetMatrix(int rows, int cols, const T *const &A,
                   int /*lda*/, T *&B, int /*ldb*/)
    {
        // cublasStatus status = cublasSetMatrix(rows, cols, sizeof(T),
        //                                      A, lda, B, ldb);

        copy(rows*cols, A, 1, B, 1);

        // check_status(status);
    }

    // @sect4{Funktion: GetMatrix}
    //
    // @param rows : Anzahl Zeilen.
    // @param cols : Anzahl Spalten.
    // @param A : Quellmatrix A.
    // @param lda : leading dimension von A.
    // @param B : Zielmatrix B.
    // @param ldb : leading dimension von B.
    template<typename T>
    static void GetMatrix(int rows, int cols, const T * const &A,
                   int lda, T *&B, int ldb)
    {
        // cublasStatus status = cublasGetMatrix(rows, cols, sizeof(T),
        //                                      A, lda, B, ldb);

        //

        for (int c = 0; c < cols; c++)
        {
            copy(rows, A+c*lda, 1, B+c*ldb, 1);


        }

        // check_status(status);
    }

    // @sect4{Funktion: SetVector}
    //
    // @param n_el : Anzahl Elemente.
    // @param src : Quellvektor src.
    // @param inc_src : Speicher Abstand zwischen Elemente in Vector src.
    // @param dst : Zielvektor dst.
    // @param inc_dst : Speicher Abstand zwischen Elemente in Vector dst.
    template<typename T>
    static void SetVector(int n_el, const T * const src, int inc_src, T *dst, int inc_dst)
    {
        // cublasStatus status = cublasSetVector(n_el, sizeof(T),
        //                                      src, inc_src, dst, inc_dst);

        // check_status(status);

        copy(n_el, src, inc_src, dst, inc_dst);

    }

    // @sect4{Funktion: GetVector}
    //
    // @param n_el : Anzahl Elemente.
    // @param A: Quellvektor A.
    // @param inc_src : Speicher Abstand zwischen Elemente in Vector A.
    // @param B : Zielvektor B.
    // @param inc_dst : Speicher Abstand zwischen Elemente in Vector B.
    template<typename T>
    static void GetVector(int n_el, const T * const &A, int inc_src, T *&B, int inc_dst)
    {
        // cublasStatus status = cublasGetVector(n_el, sizeof(T),
        //                                      A, inc_src, B, inc_dst);

        copy (n_el, A, inc_src, B, inc_dst);

        // check_status(status);
    }









    // @sect4{Funktion: axpy}
    //
    // @param n : Anzahl Elemente.
    // @param alpha: Skalar.
    // @param x : Quellvektor x.
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.
    // @param y : Zielvektor (y= alpha*x+y).
    // @param incy: Speicher Abstand zwischen Elemente in Vector y.

    static void
    axpy (int n, float alpha, const float *x,
          int incx, float *y, int incy)
    {
       cblas_saxpy(n, alpha, x, incx, y, incy);

       // cublas_status status = cublasGetError();

       // check_status(status);
    }

    static void
    axpy (int n, double alpha, const double *x,
          int incx, double *y, int incy)
    {
       cblas_daxpy(n, alpha, x, incx, y, incy);

       // cublasStatus status = cublasGetError();

       // check_status(status);
    }


    // @sect4{Funktion: copy}
    //
    // @param n : Anzahl Elemente.
    // @param x : Quellvektor x.
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.
    // @param y : Zielvektor (y= alpha*x+y).
    // @param incy: Speicher Abstand zwischen Elemente in Vector y.

    static void
    copy(int n, const float *x, int incx, float *y, int incy)
    {
        cblas_scopy(n, x, incx, y, incy);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }

    static void
    copy(int n, const double *x, int incx, double *y, int incy)
    {
        cblas_dcopy(n, x, incx, y, incy);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }

    // @sect4{Funktion: scal}
    //
    // @param n : Anzahl Elemente.
    // @param alpha: Skalar.
    // @param x : Quellvektor x (auch output).
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.


    static void
    scal (int n, float alpha, float *x, int incx)
    {
        cblas_sscal (n, alpha, x, incx);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }

    static void
    scal (int n, double alpha, double *x, int incx)
    {
        cblas_dscal (n, alpha, x, incx);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }

    // @sect4{Funktion: gemv}
    //
    // @param trans : gibt an ob A transponiert ist oder nicht. Sei trans = 'N' oder 'n' so ist op(A)= A, sei trans = 'T', 't','C' oder 'c' so ist op(A)= trans(A)
    // @param m : Anzahl Zeilen in Matrix A.
    // @param n : Anzahl Spalten in Matrix A.
    // @param alpha: Skalar fuer A.
    // @param A : Matrix A
    // @param lda : leading dimension von A.
    // @param x : Vektor mit der laenge von mindestens (1+(n-1)*abs(incx)) falls trans = 'N' oder 'n', sonst mindestens der laenge (1+(m-1)*abs(incx)).
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.
    // @param beta : Skalar fuer Vektor y.
    // @param y : Vektor mit der laenge von mindestens (1+(n-1)*abs(incy)) falls trans = 'N' oder 'n', sonst mindestens der laenge (1+(m-1)*abs(incy)).
    // @param incy: Speicher Abstand zwischen Elemente in Vector y.

    static void gemv (char trans, int m, int n, float alpha,
                      const float * const A, int lda,
                      const float * const x,  int incx, float beta,
                      float *y, int incy)
    {
        CBLAS_TRANSPOSE tr = ( ( (trans == 't') || (trans == 'T') ) ? CblasTrans : CblasNoTrans );
        cblas_sgemv (CblasColMajor, tr, m, n, alpha, A, lda, x, incx, beta, y, incy);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }

    static void gemv (char trans, int m, int n, double alpha,
                      const double * const A, int lda,
                      const double * const x,  int incx, double beta,
                      double *y, int incy)
    {
        CBLAS_TRANSPOSE tr = ( ( (trans == 't') || (trans == 'T') ) ? CblasTrans : CblasNoTrans );
        cblas_dgemv (CblasColMajor, tr, m, n, alpha, A, lda, x, incx, beta, y, incy);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }

    // @sect4{Funktion: ger}
    //
    // @param m : Anzahl Zeilen in Matrix A.
    // @param n : Anzahl Spalten in Matrix A.
    // @param alpha: Skalar fuer x*trans(y).
    // @param x : Vektor mit der laenge von mindestens der laenge (1+(m-1)*abs(incx)).
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.
    // @param y : Vektor mit der laenge von mindestens (1+(n-1)*abs(incy)).
    // @param incy: Speicher Abstand zwischen Elemente in Vector y.
    // @param A : Matrix A
    // @param lda : leading dimension von A.

    static void
    ger(int m, int n, float alpha, const float *x,
        int incx, const float *y, int incy, float *A,
        int lda)
    {
         cblas_sger(CblasColMajor, m, n, alpha, x, incx, y, incy, A, lda);

         // cublasStatus status = cublasGetError();

         // check_status(status);
    }

    static void
    ger(int m, int n, double alpha, const double *x,
        int incx, const double *y, int incy, double *A,
        int lda)
    {
         cblas_dger(CblasColMajor, m, n, alpha, x, incx, y, incy, A, lda);

         // cublasStatus status = cublasGetError();

         // check_status(status);
    }



    // @sect4{Funktion: gemm}
    //
    // @param transa : gibt an ob A transponiert ist oder nicht. Sei transa = 'N' oder 'n' so ist op(A)= A, sei trans = 'T', 't','C' oder 'c' so ist op(A)= trans(A)
    // @param transb : gibt an ob B transponiert ist oder nicht. Sei transb = 'N' oder 'n' so ist op(B)= A, sei trans = 'T', 't','C' oder 'c' so ist op(B)= trans(B)
    // @param m : Anzahl Zeilen in Matrix A und Matrix C.
    // @param n : Anzahl Spalten in Matrix B und Matrix C.
    // @param k : Anzahl Spalten in Matrix A und Zeilen in Matrix B.
    // @param alpha: Skalar fuer op(A)*op(B).
    // @param A : Matrix A
    // @param lda : leading dimension von A.
    // @param B : Matrix B.
    // @param ldb : leading dimension von B.
    // @param beta : Skalar fuer Matrix C.
    // @param C : Matrix C.
    // @param ldc : leading dimension von C.
    static void gemm(char transa, char transb, int m, int n, int k, float alpha,
                     const float * const A, int lda, const float * const B, int ldb,
                     float beta, float * C, int ldc)
    {
        CBLAS_TRANSPOSE tr_a = ( ( (transa == 't') || (transa == 'T') ) ? CblasTrans : CblasNoTrans );
        CBLAS_TRANSPOSE tr_b = ( ( (transb == 't') || (transb == 'T') ) ? CblasTrans : CblasNoTrans );

        cblas_sgemm(CblasColMajor,
                    tr_a, tr_b,
                    m, n, k,
                    alpha,
                    A, lda,
                    B, ldb,
                    beta,
                    C, ldc);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }


    static void gemm(char transa, char transb, int m, int n, int k, double alpha,
                     const double * const A, int lda, const double * const B, int ldb,
                     double beta, double * C, int ldc)
    {
        CBLAS_TRANSPOSE tr_a = ( ( (transa == 't') || (transa == 'T') ) ? CblasTrans : CblasNoTrans );
        CBLAS_TRANSPOSE tr_b = ( ( (transb == 't') || (transb == 'T') ) ? CblasTrans : CblasNoTrans );

        cblas_dgemm(CblasColMajor, tr_a, tr_b, m, n, k, alpha,
                    A, lda, B, ldb,
                    beta, C, ldc);

        // cublasStatus status = cublasGetError();

        // check_status(status);
    }


    // @sect4{Funktion: nrm2}
    //
    // @param n : Anzahl Elemente.
    // @param x : Quellvektor x.
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.

    static float
    nrm2(int n, const float *x, int incx)
    {
        float result = cblas_snrm2 (n, x, incx);

        // cublasStatus status = cublasGetError();

        // check_status(status);

        return result;
    }

    static double
    nrm2(int n, const double *x, int incx)
    {
        double result = cblas_dnrm2 (n, x, incx);

        // cublasStatus status = cublasGetError();

        // check_status(status);

        return result;
    }

    // @sect4{Funktion: dot}
    //
    // @param n : Anzahl Elemente.
    // @param x : Quellvektor x.
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.
    // @param y : Zielvektor.
    // @param incy: Speicher Abstand zwischen Elemente in Vector y.

    static float  dot(int n, const float *x, int incx, const float *y, int incy)
    {
        float result = cblas_sdot(n, x, incx, y, incy);

        // cublasStatus status = cublasGetError();

        // check_status(status);

        return result;
    }

    static double dot(int n, const double *x,
                      int incx, const double *y, int incy)
    {
        double result = cblas_ddot(n, x, incx, y, incy);

        // cublasStatus status = cublasGetError();

        // check_status(status);

        return result;
    }

}; // struct CBW END




#endif // BLAS_WRAPPER_HH

---