CUDA Lab Course Reference Manual 2011

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

Go to the documentation of this file.

#ifndef CUBLAS_WRAPPER_HH
#define CUBLAS_WRAPPER_HH

#include <cutil_inline.h>
#include <cublas.h>

// @sect3{namespace CW: cublas 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 cublas-Funktionsnamen das cublas-Praefix
// und das Genauigkeitskuerzel weglassen Z.B.: cublasSgemv wird zu gemv.
// Ansonsten siehe die CUBLAS-Doku.
struct cublas {

    // @sect4{Funktion: name}
    //
    // Gibt einen strgn zurueck, der die BLAS-Implementierug identifiziert.
    // Wird in den continuous-integration Tests fuer die Namen der Ausgabedateien
    // benoetigt.
    static inline std::string name () { return "cublas"; }

    // TODO: set und get-Routinen verpacken

    // @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)
    {
        // The following is for tracking cublas usage
        // std::cout << "Checking cublas status" << std::endl;
#ifdef DEBUG
        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
       Assert(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;

#ifndef DEBUG
        AssertThrow(s == CUBLAS_STATUS_SUCCESS,
                    ::ExcMessage("cublas init failed"));
#else
        Assert(s == CUBLAS_STATUS_SUCCESS,
               ::ExcMessage("cublas init failed"));
#endif
    }


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

        cublasStatus s = cublasShutdown();

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

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

    }


    template<typename T>
    struct Data {

        Data() : dev_ptr(0) {}

        Data(size_t n) { alloc(n); }

        ~Data() {  free_dev_ptr(); }
        T * data() { return dev_ptr; }

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

        void resize(size_t n)
        {
            free_dev_ptr();

            alloc(n);
        }

    private:
        void free_dev_ptr() {



             if (dev_ptr == 0) return;

            cublasStatus status;

                status = cublasFree(dev_ptr);

            check_status(status);

            dev_ptr = 0;

             std::cout << "GPU mem freed " << std::endl;

        }

        void alloc(size_t n)
        {
            cublasStatus status = cublasAlloc( n, sizeof(T), (void**)&dev_ptr);

            check_status(status);
        }


        T * dev_ptr;
    };


private:

    // @sect4{Funktion: Alloc}
    //
    // Ruft cublasAlloc auf und prueft anschliessend den status.
    template<typename T>
    static void Alloc(T *&dev_ptr, 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

        cublasStatus status = cublasAlloc( n, sizeof(T), (void**)&dev_ptr);

#ifdef QT_NO_DEBUG
       AssertThrow(status == CUBLAS_STATUS_SUCCESS,
                   ::ExcMessage("cublas allocation failed" ) );
#else
       AssertThrow(status == CUBLAS_STATUS_SUCCESS,
                   ::ExcMessage("cublas allocation failed" ) );
#endif

//        check_status(status);
    }

    // @sect4{Funktion: Free}
    //
    template<typename T>
    static void Free(T *&dev_ptr) {

        cublasStatus status = cublasFree(dev_ptr);

#ifdef QT_NO_DEBUG
       AssertThrow(status == CUBLAS_STATUS_SUCCESS,
                   ::ExcMessage("cublas deallocation failed" ) );
#else
       AssertThrow(status == CUBLAS_STATUS_SUCCESS,
                   ::ExcMessage("cublas deallocation failed" ) );
#endif
        //check_status(status);
    }

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);

        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);

        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);

    }

    // @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);

        check_status(status);
    }
    // @sect4{Funktion: asum}
    //
    // @param n : Anzahl Elemente.
    // @param x : Quellvektor x.
    // @param incx : Speicher Abstand zwischen Elemente in Vector x.

    static float
    asum (int n, const float *x,
          int incx)
    {
       float sum = cublasSasum(n, x, incx);

       cublasStatus status = cublasGetError();

       check_status(status);
       return sum;
    }

    static double
    asum (int n, const double *x,
          int incx)
    {
       double sum = cublasDasum(n, x, incx);

       cublasStatus status = cublasGetError();

       check_status(status);
       return sum;
    }

    static float
    asum (int n, const cuComplex *x,
          int incx)
    {
       float sum = cublasScasum(n, x, incx);

       cublasStatus status = cublasGetError();

       check_status(status);
       return sum;
    }

    static double
    asum (int n, const cuDoubleComplex *x,
          int incx)
    {
       double sum = cublasDzasum(n, x, incx);

       cublasStatus status = cublasGetError();

       check_status(status);
       return sum;
    }

    // @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)
    {
       cublasSaxpy(n, alpha, x, incx, y, incy);

       cublasStatus status = cublasGetError();

       check_status(status);
    }

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

       cublasStatus status = cublasGetError();

       check_status(status);
    }

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

       cublasStatus status = cublasGetError();

       check_status(status);
    }

    static void
    axpy (int n, cuDoubleComplex alpha, const cuDoubleComplex *x,
          int incx, cuDoubleComplex *y, int incy)
    {
       cublasZaxpy(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)
    {
        cublasScopy(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)
    {
        cublasDcopy(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);
    }
    static void
    copy(int n, const cuComplex *x, int incx, cuComplex *y, int incy)
    {
        cublasCcopy(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);
    }

    static void
    copy(int n, const cuDoubleComplex *x, int incx, cuDoubleComplex *y, int incy)
    {
        cublasZcopy(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)
    {
        cublasSscal (n, alpha, x, incx);

        cublasStatus status = cublasGetError();

        check_status(status);
    }

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

        cublasStatus status = cublasGetError();

        check_status(status);
    }



    static void
    scal (int n, cuComplex alpha, cuComplex *x, int incx)
    {
        cublasCscal (n, alpha, x, incx);

        cublasStatus status = cublasGetError();

        check_status(status);
    }



    static void
    scal (int n, cuDoubleComplex alpha, cuDoubleComplex *x, int incx)
    {
        cublasZscal (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)
    {
        cublasSgemv (trans, 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)
    {
        cublasDgemv (trans, 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)
    {
         cublasSger(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)
    {
         cublasDger(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)
    {
        cublasSgemm(transa, transb, 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)
    {
        cublasDgemm(transa, transb, 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 = cublasSnrm2 (n, x, incx);

        cublasStatus status = cublasGetError();

        check_status(status);

        return result;
    }

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

        cublasStatus status = cublasGetError();

        check_status(status);

        return result;
    }

    static cuComplex
    nrm2(int n, const cuComplex *x, int incx)
    {
        cuComplex result = dotc(n, x, incx, x, incx);
        result.x = sqrt(result.x);
        result.y = 0.;

        return result;

    }



    static cuDoubleComplex
    nrm2(int n, const cuDoubleComplex *x, int incx)
    {
        cuDoubleComplex result = dotc(n, x, incx, x, incx);
        result.x = sqrt(result.x);
        result.y = 0.;

        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 = cublasSdot(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 = cublasDdot(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);

        return result;
    }

    static cuComplex dot(int n, const cuComplex *x,
                      int incx, const cuComplex *y, int incy)
    {
        return dotc(n, x, incx, y, incy);
    }

    static cuDoubleComplex dot(int n, const cuDoubleComplex *x,
                      int incx, const cuDoubleComplex *y, int incy)
    {
        return dotc(n, x, incx, y, incy);
    }

    // @sect4{Funktion: dotu}
    //
    // @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.
    // Funktioniert auch fuer cufftComplex, da dieses per typedef
    // cuComlex cufftComplex definiert ist.
    static cuComplex dotu(int n, const cuComplex *x,
                      int incx, const cuComplex *y, int incy)
    {
        cuComplex result = cublasCdotu(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);

        std::cout << __FUNCTION__ << " : " << " c*c^* = " << result.x << ", " << result.y << std::endl;

        return result;
    }


    static cuDoubleComplex dotu(int n, const cuDoubleComplex *x,
                      int incx, const cuDoubleComplex *y, int incy)
    {
        cuDoubleComplex result = cublasZdotu(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);

        return result;
    }

    // @sect4{Funktion: dotc}
    //
    // @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.
    // Funktioniert auch fuer cufftComplex, da dieses per typedef
    // cuComlex cufftComplex definiert ist.
    static cuComplex dotc(int n, const cuComplex *x,
                      int incx, const cuComplex *y, int incy)
    {
        cuComplex result = cublasCdotc(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);

        //std::cout << __FUNCTION__ << " : " << " c*c^* = " << result.x << ", " << result.y << std::endl;

        return result;
    }


    static cuDoubleComplex dotc(int n, const cuDoubleComplex *x,
                      int incx, const cuDoubleComplex *y, int incy)
    {
        cuDoubleComplex result = cublasZdotc(n, x, incx, y, incy);

        cublasStatus status = cublasGetError();

        check_status(status);

        return result;
    }


    // @sect4{Funktion: trsm}
    //
    inline void trsm(char side, char uplo, char transa, char diag,  int m, int n,
                     float alpha,
                     const float * A, int lda, float * B, int ldb)
    {
        cublasStrsm(side, uplo, transa, diag, m, n, alpha,
                    A, lda, B, ldb);

        cublasStatus status = cublasGetError();

        check_status(status);
    }

    inline void trsm(char side, char uplo, char transa, char diag,  int m, int n,
                     double alpha,
                     const double * A, int lda, double * B, int ldb)
    {
        cublasDtrsm(side, uplo, transa, diag, m, n, alpha,
                    A, lda, B, ldb);

        cublasStatus status = cublasGetError();

        check_status(status);
    }

}; // struct CW END




#endif // CUBLAS_WRAPPER_HH

---