Code: Select all
int magma_dgeqrs_gpu(int *m, int *n, int *nrhs,
double *a, int *lda, double *tau, double *c, int *ldc,
double *work, int *lwork, double *td, int *info);
Code: Select all
int magma_dgeqrf_gpu2(int *, int *, double *, int *, double *,
double *, int *, double *, int *);
Code: Select all
#include <stdio.h>
#include <stdlib.h>
#include <cuda.h>
#include <cublas.h>
#include <magma.h>
/***************************************************************************
* Program to compute the inverse of a symmetric positive-definite matrix A
* using QR factorization.
* - magma_dgeqrs_gpu fails to compute the inverse correctly
* - LAPACK routines do it correctly
***************************************************************************/
void mprint(int m, int n, double *x, int ldx) {
int i, j;
for(i = 0; i < m; i++) {
for(j = 0; j < n; j++) {
printf("%f ", x[i + j * ldx]);
}
printf("\n");
}
}
int main(int argc, char** argv)
{
int m = 3, n = 3, nrhs = 3, lda = m, ldb = m,
nb = magma_get_dgeqrf_nb(m),
lhwork1 = (m - n + nb + 2 * nrhs) * nb,
lhwork2 = lhwork1,
ldwork = lhwork1,
info, i, j;
double A[9] = {1, 0.4, 0.2, 0.4, 1.0, 0.3, 0.2, 0.3, 1},
B[9] = {1, 0, 0, 0, 1, 0, 0, 0, 1},
Ai[9], AAi[9], QR[9], tau[3], *d_A, *d_B, *h_work1, *h_work2, *d_work,
alpha, beta;
printf("Matrix A to be inverted:\n");
mprint(m, n, A, lda);
cublasInit();
cublasAlloc(m * n, sizeof(double), (void**)&d_A);
cublasAlloc(m * nrhs, sizeof(double), (void**)&d_B);
cublasAlloc(ldwork, sizeof(double), (void**)&d_work);
cudaMallocHost((void**)&h_work1, lhwork1 * sizeof(double));
cudaMallocHost((void**)&h_work2, lhwork2 * sizeof(double));
if(d_A == NULL || d_B == NULL || d_work == NULL|| h_work2 == NULL
|| h_work1 == NULL ) {
fprintf(stderr, "Error: CUDA memory allocation failed\n");
exit(1);
}
/***************************************************************************
* MAGMA QR factorization
* - Correctly returns factorization matrix (d_A/QR) and scalars (tau)
***************************************************************************/
printf("\nMAGMA Results\n=============\n");
cublasSetVector(m * n, sizeof(double), A, 1, d_A, 1);
magma_dgeqrf_gpu2(&m, &n, d_A, &lda, tau, h_work1, &lhwork1, d_work, &info);
printf("\nmagma_dgeqrf_gpu2 info = %d\n", info);
// Save d_A factorization to QR for later use in LAPACK routines
cublasGetVector(m * n, sizeof(double), d_A, 1, QR, 1);
/***************************************************************************
* MAGMA QR solver
* - Returns incorrect solution
***************************************************************************/
cublasSetVector(m * nrhs, sizeof(double), B, 1, d_B, 1);
magma_dgeqrs_gpu(&m, &n, &nrhs, d_A, &lda, tau, d_B, &ldb, h_work2, &lhwork2,
d_work, &info);
cublasGetMatrix(n, nrhs, sizeof(double), d_B, ldb, Ai, lda);
printf("magma_dgeqrs_gpu info = %d\n", info);
printf("\nA^{-1} (magma_dgeqrs_gpu):\n");
mprint(n, nrhs, Ai, lda);
alpha = 1.0; beta = 0.0;
dgemm_("N", "N", &m, &n, &nrhs, &alpha, A, &lda, Ai, &lda, &beta, AAi, &lda);
printf("\nA %%*%% A^{-1}:\n");
mprint(m, nrhs, AAi, lda);
/***************************************************************************
* LAPACK QR solver
* - Returns correct solution
***************************************************************************/
printf("\nLAPACK Results\n==============\n");
alpha = 1.0;
dormqr_("L", "T", &m, &nrhs, &n, QR, &lda, tau, B, &ldb, h_work2, &lhwork2, &info);
dtrsm_("L", "U", "N", "N", &m, &nrhs, &alpha, QR, &lda, B, &ldb);
printf("\ndrmqr info = %d\n", info);
printf("\nA^{-1} (dormqr/dtrsm):\n");
mprint(n, nrhs, B, ldb);
alpha = 1.0; beta = 0.0;
dgemm_("N", "N", &m, &n, &nrhs, &alpha, A, &lda, B, &ldb, &beta, AAi, &lda);
printf("\nA %%*%% A^{-1}:\n");
mprint(m, nrhs, AAi, lda);
cublasFree(d_A);
cublasFree(d_B);
cublasFree(d_work);
cudaFreeHost(h_work1);
cudaFreeHost(h_work2);
cublasShutdown();
return(EXIT_SUCCESS);
}