39 USAGE(
"HEEV",
"N LDA",
40 " - N : size of the matrix A\n"
41 " - LDA : leading dimension of the matrix A\n");
45 int N = atoi(argv[0]);
46 int LDA = atoi(argv[1]);
49 float eps = LAPACKE_slamch_work(
'e');
56 int info_solution = 0;
57 int info_reduction = 0;
65 float *W1 = (
float *)malloc(N*
sizeof(
float));
66 float *W2 = (
float *)malloc(N*
sizeof(
float));
71 if ( (!A2) || (!W1) || (!W2) ){
72 printf(
"Out of Memory \n ");
91 LAPACKE_clatms_work( LAPACK_COL_MAJOR, N, N,
102 LAPACKE_slasrt_work(
'I', N, W1 );
104 if (getenv(
"PLASMA_TESTING_VERBOSE"))
106 printf(
"Eigenvalues original\n");
107 for (i = 0; i <
min(N,25); i++){
108 printf(
"%f \n", W1[i]);
118 LAPACKE_clacpy_work(LAPACK_COL_MAJOR,
'A', N, N, A2, LDA, A1, LDA);
124 if (getenv(
"PLASMA_TESTING_VERBOSE"))
126 printf(
"Eigenvalues computed\n");
127 for (i = 0; i <
min(N,25); i++){
128 printf(
"%f \n", W2[i]);
134 printf(
"------ TESTS FOR PLASMA CHEEV ROUTINE ------- \n");
135 printf(
" Size of the Matrix %d by %d\n", N, N);
137 printf(
" The matrix A is randomly generated for each test.\n");
138 printf(
"============\n");
139 printf(
" The relative machine precision (eps) is to be %e \n",eps);
140 printf(
" Computational tests pass if scaled residuals are less than 60.\n");
153 if ( (info_solution == 0) & (info_ortho == 0) & (info_reduction == 0) ) {
154 printf(
"***************************************************\n");
155 printf(
" ---- TESTING CHEEV ...................... PASSED !\n");
156 printf(
"***************************************************\n");
159 printf(
"************************************************\n");
160 printf(
" - TESTING CHEEV ... FAILED !\n");
161 printf(
"************************************************\n");
169 if (Q != NULL) free(Q);
170 if (A1 != NULL) free(A1);
185 int minMN =
min(M, N);
186 float *work = (
float *)malloc(minMN*
sizeof(
float));
190 LAPACKE_claset_work(LAPACK_COL_MAJOR,
'A', minMN, minMN, 0., 1., Id, minMN);
194 cblas_cherk(
CblasColMajor,
CblasUpper,
CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N);
196 cblas_cherk(
CblasColMajor,
CblasUpper,
CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);
200 result = normQ / (minMN * eps);
201 printf(
"============\n");
202 printf(
"Checking the orthogonality of Q \n");
203 printf(
"||Id-Q'*Q||_oo / (minMN*eps) = %e \n", result);
205 if ( isnan(result) || isinf(result) || (result > 60.0) ) {
206 printf(
"-- Orthogonality is suspicious ! \n");
210 printf(
"-- Orthogonality is CORRECT ! \n");
214 free(work); free(Id);
227 float Anorm, Rnorm, result;
234 float *work = (
float *)malloc(N*
sizeof(
float));
239 LAPACKE_clacpy_work(LAPACK_COL_MAJOR,
lapack_const(uplo), N, N, A2, LDA, T, N);
241 printf(
"============\n");
242 printf(
"Checking the tridiagonal reduction \n");
247 for (i = 0; i < N; i++)
248 for (j =
max(0, i-bw); j < i; j++)
249 T[i*N+j] = conjf(T[j*N+i]);
253 cblas_cgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, N, N, N,
CBLAS_SADDR(alpha), Q, LDA, T, N,
CBLAS_SADDR(beta), Aorig, N);
254 cblas_cgemm(
CblasColMajor,
CblasNoTrans,
CblasConjTrans, N, N, N,
CBLAS_SADDR(alpha), Aorig, N, Q, LDA,
CBLAS_SADDR(beta), T, N);
259 for (i = 0; i < N; i++)
260 for (j = i+1 ; j <
min(i+bw, N); j++)
261 T[i*N+j] = conjf(T[j*N+i]);
265 cblas_cgemm(
CblasColMajor,
CblasConjTrans,
CblasNoTrans, N, N, N,
CBLAS_SADDR(alpha), Q, LDA, T, N,
CBLAS_SADDR(beta), Aorig, N);
266 cblas_cgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, N, N, N,
CBLAS_SADDR(alpha), Aorig, N, Q, LDA,
CBLAS_SADDR(beta), T, N);
270 for (i = 0; i < N; i++)
271 for (j = 0 ; j < N; j++)
272 Residual[j*N+i] = A1[j*LDA+i] - T[j*N+i];
277 result = Rnorm / ( Anorm * N * eps);
279 printf(
"-- ||A-Q*T*Q'||_oo/(||A||_oo.N.eps) = %e \n", result);
281 printf(
"-- ||A-Q'*T*Q||_oo/(||A||_oo.N.eps) = %e \n", result);
283 if ( isnan(result) || isinf(result) || (result > 60.0) ) {
284 printf(
"-- Reduction is suspicious ! \n");
288 printf(
"-- Reduction is CORRECT ! \n");
292 free(Aorig); free(Residual); free(T);
294 return info_reduction;
300 int info_solution, i;
303 float maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
304 float maxeig =
max( fabs(E1[0]), fabs(E2[0]) );
305 for (i = 1; i < N; i++){
306 resid = fabs(fabs(E1[i])-fabs(E2[i]));
307 maxtmp =
max(fabs(E1[i]), fabs(E2[i]));
310 maxeig =
max(maxtmp, maxeig);
311 maxel =
max(resid, maxel );
314 maxel = maxel / (maxeig * N * eps);
315 printf(
" ======================================================\n");
316 printf(
" | D - eigcomputed | / (|D| * N * eps) : %15.3E \n", maxel );
317 printf(
" ======================================================\n");
319 printf(
"============\n");
320 printf(
"Checking the eigenvalues of A\n");
321 if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
322 printf(
"-- The eigenvalues are suspicious ! \n");
326 printf(
"-- The eigenvalues are CORRECT ! \n");
330 return info_solution;