28 #define max(a, b) ((a) > (b) ? (a) : (b))
31 #define min(a, b) ((a) < (b) ? (a) : (b))
34 int check_solution(
int,
int,
int,
double*,
int,
double*,
double*,
int);
53 int LDBxNRHS = LDB*NRHS;
55 double *A1 = (
double *)malloc(LDA*N*
sizeof(
double));
56 double *A2 = (
double *)malloc(LDA*N*
sizeof(
double));
57 double *B1 = (
double *)malloc(LDB*NRHS*
sizeof(
double));
58 double *B2 = (
double *)malloc(LDB*NRHS*
sizeof(
double));
62 if ((!A1)||(!A2)||(!B1)||(!B2)){
63 printf(
"Out of Memory \n ");
69 printf(
"-- PLASMA is initialized to run on %d cores. \n",cores);
75 LAPACKE_dlarnv_work(
IONE,
ISEED, LDAxN, A1);
76 for (i = 0; i < M; i++)
77 for (j = 0; j < N; j++)
78 A2[LDA*j+i] = A1[LDA*j+i] ;
81 LAPACKE_dlarnv_work(
IONE,
ISEED, LDBxNRHS, B1);
82 for (i = 0; i < M; i++)
83 for (j = 0; j < NRHS; j++)
84 B2[LDB*j+i] = B1[LDB*j+i] ;
92 N, NRHS, (
double)1.0, A2, LDA, B2, LDB);
97 if ((info_solution != 0)|(info != 0))
98 printf(
"-- Error in DORMQR example ! \n");
100 printf(
"-- Run of DORMQR example successful ! \n");
102 free(A1); free(A2); free(B1); free(B2); free(T);
113 int check_solution(
int M,
int N,
int NRHS,
double *A1,
int LDA,
double *B1,
double *B2,
int LDB)
116 double Rnorm, Anorm, Xnorm, Bnorm;
118 double *work = (
double *)malloc(
max(M, N)*
sizeof(double));
121 eps = LAPACKE_dlamch_work(
'e');
130 cblas_dgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, NRHS, N, (alpha), A1, LDA, B2, LDB, (beta), B1, LDB);
133 double *Residual = (
double *)malloc(M*NRHS*
sizeof(
double));
134 memset((
void*)Residual, 0, M*NRHS*
sizeof(
double));
135 cblas_dgemm(
CblasColMajor,
CblasTrans,
CblasNoTrans, N, NRHS, M, (alpha), A1, LDA, B1, LDB, (beta), Residual, M);
140 double *Residual = (
double *)malloc(N*NRHS*
sizeof(
double));
141 memset((
void*)Residual, 0, N*NRHS*
sizeof(
double));
142 cblas_dgemm(
CblasColMajor,
CblasTrans,
CblasNoTrans, N, NRHS, M, (alpha), A1, LDA, B1, LDB, (beta), Residual, N);
147 printf(
"============\n");
148 printf(
"Checking the Residual of the solution \n");
149 printf(
"-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||)_oo.N.eps) = %e \n",Rnorm/((Anorm*Xnorm+Bnorm)*N*eps));
151 if (isnan(Rnorm / ((Anorm * Xnorm + Bnorm) * N * eps)) || (Rnorm / ((Anorm * Xnorm + Bnorm) * N * eps) > 10.0) ) {
152 printf(
"-- The solution is suspicious ! \n");
156 printf(
"-- The solution is CORRECT ! \n");
162 return info_solution;