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,
float*,
int,
float*,
float*,
int);
53 int LDBxNRHS = LDB*NRHS;
55 float *A1 = (
float *)malloc(LDA*N*
sizeof(
float));
56 float *A2 = (
float *)malloc(LDA*N*
sizeof(
float));
57 float *B1 = (
float *)malloc(LDB*NRHS*
sizeof(
float));
58 float *B2 = (
float *)malloc(LDB*NRHS*
sizeof(
float));
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_slarnv_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_slarnv_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] ;
95 if ((info_solution != 0)|(info != 0))
96 printf(
"-- Error in SGELQS example ! \n");
98 printf(
"-- Run of SGELQS example successful ! \n");
100 free(A1); free(A2); free(B1); free(B2); free(T);
111 int check_solution(
int M,
int N,
int NRHS,
float *A1,
int LDA,
float *B1,
float *B2,
int LDB)
114 float Rnorm, Anorm, Xnorm, Bnorm;
116 float *work = (
float *)malloc(
max(M, N)*
sizeof(float));
119 eps = LAPACKE_slamch_work(
'e');
128 cblas_sgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, NRHS, N, (alpha), A1, LDA, B2, LDB, (beta), B1, LDB);
131 float *Residual = (
float *)malloc(M*NRHS*
sizeof(
float));
132 memset((
void*)Residual, 0, M*NRHS*
sizeof(
float));
133 cblas_sgemm(
CblasColMajor,
CblasTrans,
CblasNoTrans, N, NRHS, M, (alpha), A1, LDA, B1, LDB, (beta), Residual, M);
138 float *Residual = (
float *)malloc(N*NRHS*
sizeof(
float));
139 memset((
void*)Residual, 0, N*NRHS*
sizeof(
float));
140 cblas_sgemm(
CblasColMajor,
CblasTrans,
CblasNoTrans, N, NRHS, M, (alpha), A1, LDA, B1, LDB, (beta), Residual, N);
145 printf(
"============\n");
146 printf(
"Checking the Residual of the solution \n");
147 printf(
"-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||)_oo.N.eps) = %e \n",Rnorm/((Anorm*Xnorm+Bnorm)*N*eps));
149 if (isnan(Rnorm / ((Anorm * Xnorm + Bnorm) * N * eps)) || (Rnorm / ((Anorm * Xnorm + Bnorm) * N * eps) > 10.0) ) {
150 printf(
"-- The solution is suspicious ! \n");
154 printf(
"-- The solution is CORRECT ! \n");
160 return info_solution;