// tutorial3_blas.cc
#include <exception>

#include "magma_v2.h"

// -----------------------------------------------------------------------------
// generates random m-by-k dA, k-by-n dB, and m-by-n dC in GPU device memory
void generate_problem( int m, int n, int k,
                       double* dA, int ldda,
                       double* dB, int lddb,
                       double* dC, int lddc )
{
    int lda = ldda, ldb = lddb, ldc = lddc;
    double *A = new double[ lda*k ];
    double *B = new double[ ldb*n ];
    double *C = new double[ ldc*n ];
    for (int j = 0; j < k; ++j) {
        for (int i = 0; i < m; ++i) {
            A[ i + j*lda ] = rand() / double(RAND_MAX);
        }
    }
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < k; ++i) {
            B[ i + j*ldb ] = rand() / double(RAND_MAX);
        }
    }
    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < m; ++i) {
            C[ i + j*ldc ] = rand() / double(RAND_MAX);
        }
    }

    magma_queue_t queue;
    int device;
    magma_getdevice( &device );
    magma_queue_create( device, &queue );
    magma_dsetmatrix( m, k, A, lda, dA, ldda, queue );
    magma_dsetmatrix( k, n, B, ldb, dB, lddb, queue );
    magma_dsetmatrix( m, n, C, ldc, dC, lddc, queue );
    magma_queue_destroy( queue );

    delete[] A;
    delete[] B;
    delete[] C;
}

// -----------------------------------------------------------------------------
// Matrix-matrix multiply (gemm).
int main( int argc, char** argv )
{
    // ... setup matrices on GPU:
    // m-by-k matrix dA,
    // k-by-n matrix dB,
    // m-by-n matrix dC.
    int m = 100, k = 50, n = 200;
    int ldda = m, lddb = k, lddc = m;
    double *dA, *dB, *dC;
    magma_dmalloc( &dA, ldda*n );
    magma_dmalloc( &dB, lddb*n );
    magma_dmalloc( &dC, lddc*n );
    assert( dA != nullptr );
    assert( dB != nullptr );
    assert( dC != nullptr );
    generate_problem( m, n, k, dA, ldda, dB, lddb, dC, lddc );

    int device;
    magma_queue_t queue;
    magma_getdevice( &device );
    magma_queue_create( device, &queue );

    // C = -A B + C
    magma_dgemm( MagmaNoTrans,
                 MagmaNoTrans, m, n, k,
                 -1.0, dA, ldda,
                       dB, lddb,
                  1.0, dC, lddc, queue );

    // ... do concurrent work on CPU

    // wait for gemm to finish
    magma_queue_sync( queue );

    // ... use result in dC

    magma_queue_destroy( queue );

    // ... cleanup
    magma_free( dA );
    magma_free( dB );
    magma_free( dC );
}