// tutorial2_gpu_interface.cc.cc
#include <exception>

#include "magma_v2.h"

// -----------------------------------------------------------------------------
// generates random n-by-n dA and n-by-nrhs dX in GPU device memory
void generate_problem( int n, int nrhs,
                       double* dA, int ldda,
                       double* dX, int lddx )
{
    int lda = ldda, ldx = lddx;
    double *A = new double[ lda*n ];
    double *X = new double[ ldx*nrhs ];

    for (int j = 0; j < n; ++j) {
        for (int i = 0; i < n; ++i) {
            A[ i + j*lda ] = rand() / double(RAND_MAX);
        }
    }
    for (int j = 0; j < nrhs; ++j) {
        for (int i = 0; i < n; ++i) {
            X[ i + j*ldx ] = rand() / double(RAND_MAX);
        }
    }

    magma_queue_t queue;
    int device;
    magma_getdevice( &device );
    magma_queue_create( device, &queue );
    magma_dsetmatrix( n, n,    A, lda, dA, ldda, queue );
    magma_dsetmatrix( n, nrhs, X, ldx, dX, lddx, queue );
    magma_queue_destroy( queue );

    delete[] A;
    delete[] X;
}

// -----------------------------------------------------------------------------
// GPU interface: input & output matrices dA and dX in GPU device memory
int main( int argc, char** argv )
{
    magma_init();

    int n = 100, nrhs = 10;
    int ldda = magma_roundup( n, 32 );
    int lddx = magma_roundup( n, 32 );
    int* ipiv = new int[ n ];

    double *dA, *dX;
    magma_dmalloc( &dA, ldda*n );
    magma_dmalloc( &dX, lddx*nrhs );
    assert( dA != nullptr );
    assert( dX != nullptr );

    // ... fill in dA and dX (on GPU)
    generate_problem( n, nrhs, dA, ldda, dX, lddx );

    // solve AX = B where B is in X
    int info;
    magma_dgesv_gpu( n, nrhs,
                     dA, ldda, ipiv,
                     dX, lddx, &info );
    if (info != 0) {
        throw std::exception();
    }

    // ... use result in dX

    magma_free( dA );
    magma_free( dX );
    delete[] ipiv;

    magma_finalize();
}