PLASMA  2.4.5
PLASMA - Parallel Linear Algebra for Scalable Multi-core Architectures
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Groups
testing_ssygv.c
Go to the documentation of this file.
1 
15 #include <stdlib.h>
16 #include <stdio.h>
17 #include <string.h>
18 #include <math.h>
19 
20 #include <plasma.h>
21 #include <cblas.h>
22 #include <lapacke.h>
23 #include <core_blas.h>
24 #include "testing_smain.h"
25 
26 #undef COMPLEX
27 #define REAL
28 
29 static int check_orthogonality(int, int, float*, int, float);
30 static int check_reduction(int, int, int, int, float*, float*, int, float*, int, float*, float);
31 #if 0
32 static int check_solution(int, int, int, float*, int, float*, float*, int, float);
33 #endif
34 static int check_solution(int, float*, float*, float);
35 
36 int testing_ssygv(int argc, char **argv)
37 {
38  /* Check for number of arguments*/
39  if (argc != 3) {
40  USAGE("HEGV", "N LDA LDB",
41  " - N : size of the matrices A and B\n"
42  " - LDA : leading dimension of the matrix A\n"
43  " - LDB : leading dimension of the matrix B\n");
44  return -1;
45  }
46 
47  float eps = LAPACKE_slamch_work('e');
48  PLASMA_enum vec = PlasmaNoVec;
49  int N = atoi(argv[0]);
50  int LDA = atoi(argv[1]);
51  int LDB = atoi(argv[2]);
52  int LDQ = LDA;
53  int LDAxN = LDA*N;
54  int LDBxN = LDB*N;
55  int LDQxN = LDQ*N;
56 
57  int info_ortho = 0;
58  int info_solution = 0;
59  int info_reduction = 0;
60  int i, u;
61 
62  float *A1 = (float *)malloc(LDAxN*sizeof(float));
63  float *A2 = (float *)malloc(LDAxN*sizeof(float));
64  float *B1 = (float *)malloc(LDBxN*sizeof(float));
65  float *B2 = (float *)malloc(LDBxN*sizeof(float));
66  float *Q = (float *)malloc(LDQxN*sizeof(float));
67  float *Ainit = (float *)malloc(LDAxN*sizeof(float));
68  float *Binit = (float *)malloc(LDBxN*sizeof(float));
69  float *W1 = (float *)malloc(N*sizeof(float));
70  float *W2 = (float *)malloc(N*sizeof(float));
71  float *work = (float *)malloc(3*N* sizeof(float));
72  PLASMA_desc *T;
73 
74  /* Check if unable to allocate memory */
75  if ((!A1)||(!A2)||(!B1)||(!B2)||(!Q)||(!Ainit)||(!Binit)){
76  printf("Out of Memory \n ");
77  return -2;
78  }
79 
80  /*
81  PLASMA_Disable(PLASMA_AUTOTUNING);
82  PLASMA_Set(PLASMA_TILE_SIZE, 120);
83  PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, 20);
84  */
85 
86  PLASMA_Enable(PLASMA_WARNINGS);
87  PLASMA_Enable(PLASMA_ERRORS);
88  PLASMA_Alloc_Workspace_ssygv(N, N, &T);
89 
90  /*----------------------------------------------------------
91  * TESTING SSYGV
92  */
93 
94  /* Initialize A1 and Ainit */
95  PLASMA_splgsy(0., N, A1, LDA, 5198);
96  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, A1, LDA, Ainit, LDA);
97 
98  /* Initialize B1 and Binit */
99  PLASMA_splgsy((float)N, N, B1, LDB, 4321 );
100  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, B1, LDB, Binit, LDB);
101 
102  printf("\n");
103  printf("------ TESTS FOR PLASMA SSYGV ROUTINE ------- \n");
104  printf(" Size of the Matrix %d by %d\n", N, N);
105  printf("\n");
106  printf(" The matrix A is randomly generated for each test.\n");
107  printf("============\n");
108  printf(" The relative machine precision (eps) is to be %e \n",eps);
109  printf(" Computational tests pass if scaled residuals are less than 60.\n");
110 
111  /*----------------------------------------------------------
112  * TESTING SSYGV
113  */
114 
115  for (i=0; i<3; i++) {
116  for (u=0; u<2; u++) {
117  LAPACKE_slaset_work(LAPACK_COL_MAJOR, 'A', LDA, N, 0., 1., Q, LDA);
118 
119  memcpy(A2, Ainit, LDAxN*sizeof(float));
120  memcpy(B2, Binit, LDBxN*sizeof(float));
121 
122  PLASMA_ssygv(itype[i], vec, uplo[u], N, A2, LDA, B2, LDB, W2, T, Q, LDQ);
123 
124  /* Check the orthogonality, reduction and the eigen solutions */
125  if (vec == PlasmaVec)
126  info_ortho = check_orthogonality(N, N, Q, LDA, eps);
127  /*
128  * WARNING: For now, Q is associated to Band tridiagonal reduction and
129  * not to the final tridiagonal reduction, so we can not call the check
130  */
131  if (0)
132  info_reduction = check_reduction(itype[i], uplo[u], N, 1, A1, A2, LDA, B2, LDB, Q, eps);
133 
134  memcpy(A1, Ainit, LDAxN*sizeof(float));
135  memcpy(B1, Binit, LDBxN*sizeof(float));
136 
137  LAPACKE_ssygv( LAPACK_COL_MAJOR,
138  itype[i], lapack_const(vec), lapack_const(uplo[u]),
139  N, A1, LDA, B1, LDB, W1);
140 
141  /*info_solution = check_solution(N, N, N, A1, LDA, B1, B2, LDB, eps);*/
142  info_solution = check_solution(N, W1, W2, eps);
143 
144  if ( (info_ortho == 0) & (info_reduction == 0) & (info_solution == 0)) {
145  printf("***************************************************\n");
146  printf(" ---- TESTING SSYGV (%s, %s) ...................... PASSED !\n", itypestr[i], uplostr[u]);
147  printf("***************************************************\n");
148  }
149  else {
150  printf("************************************************\n");
151  printf(" - TESTING SSYGV (%s, %s) ... FAILED !\n", itypestr[i], uplostr[u]);
152  printf("************************************************\n");
153  }
154  }
155  }
156 
157  PLASMA_Dealloc_Handle_Tile(&T);
158  free(A1);
159  free(A2);
160  free(B1);
161  free(B2);
162  free(Q);
163  free(Ainit);
164  free(Binit);
165  free(W1);
166  free(W2);
167  free(work);
168 
169  return 0;
170 }
171 
172 /*-------------------------------------------------------------------
173  * Check the orthogonality of Q
174  */
175 static int check_orthogonality(int M, int N, float *Q, int LDQ, float eps)
176 {
177  float alpha = 1.0;
178  float beta = -1.0;
179  float normQ, result;
180  int info_ortho;
181  int minMN = min(M, N);
182  float *work = (float *)malloc(minMN*sizeof(float));
183 
184  /* Build the idendity matrix */
185  float *Id = (float *) malloc(minMN*minMN*sizeof(float));
186  LAPACKE_slaset_work(LAPACK_COL_MAJOR, 'A', minMN, minMN, 0., 1., Id, minMN);
187 
188  /* Perform Id - Q'Q */
189  if (M >= N)
190  cblas_ssyrk(CblasColMajor, CblasUpper, CblasTrans, N, M, alpha, Q, LDQ, beta, Id, N);
191  else
192  cblas_ssyrk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);
193 
194  normQ = LAPACKE_slansy_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), 'U', minMN, Id, minMN, work);
195 
196  result = normQ / (minMN * eps);
197  printf("============\n");
198  printf("Checking the orthogonality of Q \n");
199  printf("||Id-Q'*Q||_oo / (minMN*eps) = %e \n", result);
200 
201  if ( isnan(result) || isinf(result) || (result > 60.0) ) {
202  printf("-- Orthogonality is suspicious ! \n");
203  info_ortho=1;
204  }
205  else {
206  printf("-- Orthogonality is CORRECT ! \n");
207  info_ortho=0;
208  }
209 
210  free(work); free(Id);
211 
212  return info_ortho;
213 }
214 
215 /*------------------------------------------------------------
216  * Check the tridiagonal reduction
217  */
218 
219 static int check_reduction(int itype, int uplo, int N, int bw,
220  float *A1, float *A2, int LDA,
221  float *B2, int LDB, float *Q, float eps )
222 {
223  float alpha = 1.0;
224  float beta = 0.0;
225  float Anorm, Rnorm, result;
226  int info_reduction;
227  int i, j;
228  char *str;
229 
230  float *Aorig = (float *)malloc(N*N*sizeof(float));
231  float *Residual = (float *)malloc(N*N*sizeof(float));
232  float *T = (float *)malloc(N*N*sizeof(float));
233  float *work = (float *)malloc(N*sizeof(float));
234 
235  memset((void*)T, 0, N*N*sizeof(float));
236  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, ' ', N, N, A1, LDA, Residual, N);
237 
238  /* Rebuild the T */
239  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, lapack_const(uplo), N, N, A2, LDA, T, N);
240 
241  if (uplo == PlasmaLower) {
242  /* Set the reflectors to 0 */
243  for (i = bw+1; i < N; i++)
244  for (j = 0 ; (j < N) && (j < i-bw); j++)
245  T[j*N+i] = 0.;
246 
247  /* Copy the lower part to the upper part to rebuild the symmetry */
248  for (i = 0; i < N; i++)
249  for (j = 0 ; j < i; j++)
250  T[i*N+j] = (T[j*N+i]);
251  } else {
252  /* Set the reflectors to 0 */
253  for (j = bw+1; j < N; j++)
254  for (i = 0 ; (i < N) && (i < j-bw); i++)
255  T[j*N+i] = 0.;
256  /* Copy the upper part to the lower part to rebuild the symmetry */
257  for (i = 0; i < N; i++)
258  for (j = i+1 ; j < N; j++)
259  T[i*N+j] = (T[j*N+i]);
260  }
261 
262  memset((void*)Aorig, 0, N*N*sizeof(float));
263 
264  if (itype == 1) {
265  if (uplo == PlasmaLower) {
266  str = "L*Q*T*Q'*L'";
267 
268  /* Compute Aorig=Q*T*Q' */
269  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, N, N, (alpha), Q, LDA, T, N, (beta), Aorig, N);
270  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasTrans, N, N, N, (alpha), Aorig, N, Q, LDA, (beta), T, N);
271 
272  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, PlasmaUpperLower, N, N, T, N, Aorig, N);
273 
274  cblas_strmm(CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
275  cblas_strmm(CblasColMajor, CblasRight, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
276 
277  }
278  else {
279  str = "U'*Q*T*Q'*U";
280 
281  /* Compute Aorig=Q'*T*Q */
282  cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, N, (alpha), Q, LDA, T, N, (beta), Aorig, N);
283  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, N, N, (alpha), Aorig, N, Q, LDA, (beta), T, N);
284 
285  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, T, N, Aorig, N);
286 
287  cblas_strmm(CblasColMajor, CblasLeft, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
288  cblas_strmm(CblasColMajor, CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
289  }
290  }
291  else {
292  if (uplo == PlasmaLower) {
293  str = "inv(L')*Q*A2*Q'*inv(L)";
294 
295  /* Compute Aorig=Q*T*Q' */
296  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, N, N, (alpha), Q, LDA, T, N, (beta), Aorig, N);
297  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasTrans, N, N, N, (alpha), Aorig, N, Q, LDA, (beta), T, N );
298 
299  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, T, N, Aorig, N);
300 
301  cblas_strsm(CblasColMajor, CblasLeft, CblasLower, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
302  cblas_strsm(CblasColMajor, CblasRight, CblasLower, CblasNoTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
303  }
304  else{
305  str = "inv(U)*Q*A2*Q'*inv(U')";
306 
307  /* Compute Aorig=Q'*T*Q */
308  cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, N, N, (alpha), Q, LDA, T, N, (beta), Aorig, N);
309  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, N, N, (alpha), Aorig, N, Q, LDA, (beta), T, N);
310 
311  LAPACKE_slacpy_work(LAPACK_COL_MAJOR, 'A', N, N, T, N, Aorig, N);
312 
313  cblas_strsm(CblasColMajor, CblasLeft, CblasUpper, CblasNoTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
314  cblas_strsm(CblasColMajor, CblasRight, CblasUpper, CblasTrans, CblasNonUnit, N, N, (alpha), B2, LDB, Aorig, N);
315  }
316  }
317 
318  /* Compute the Residual */
319  for (i = 0; i < N; i++)
320  for (j = 0 ; j < N; j++)
321  Residual[j*N+i] = A1[j*LDA+i]-Aorig[j*N+i];
322 
323  Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, Residual, N, work);
324  Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, N, A2, LDA, work);
325 
326  result = Rnorm / (Anorm * N *eps);
327  printf("============\n");
328  printf("Checking the tridiagonal reduction \n");
329  printf("-- ||A-%s||_oo/(||A||_oo.N.eps) = %e \n", str, result );
330 
331  if (isnan(result) || isinf(result) || (result > 60.0) ) {
332  printf("-- Reduction is suspicious ! \n");
333  info_reduction = 1;
334  }
335  else {
336  printf("-- Reduction is CORRECT ! \n");
337  info_reduction = 0;
338  }
339 
340  free(Aorig); free(Residual); free(T); free(work);
341 
342  return info_reduction;
343 }
344 
345 /*--------------------------------------------------------------
346  * Check the solution
347  */
348 static int check_solution(int N, float *E1, float *E2, float eps)
349 {
350  int info_solution, i;
351  float resid;
352  float maxtmp;
353  float maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
354  float maxeig = max( fabs(E1[0]), fabs(E2[0]) );
355  for (i = 1; i < N; i++){
356  resid = fabs(fabs(E1[i])-fabs(E2[i]));
357  maxtmp = max(fabs(E1[i]), fabs(E2[i]));
358 
359  /* Update */
360  maxeig = max(maxtmp, maxeig);
361  maxel = max(resid, maxel );
362  }
363 
364  maxel = maxel / (maxeig * N * eps);
365  printf(" ======================================================\n");
366  printf(" | D - eigcomputed | / (|D| * N * eps) : %15.3E \n", maxel );
367  printf(" ======================================================\n");
368 
369  printf("============\n");
370  printf("Checking the eigenvalues of A\n");
371  if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
372  printf("-- The eigenvalues are suspicious ! \n");
373  info_solution = 1;
374  }
375  else{
376  printf("-- The eigenvalues are CORRECT ! \n");
377  info_solution = 0;
378  }
379 
380  return info_solution;
381 }
382 
383 #if 0
384 static int check_solution(int M, int N, int NRHS, float *A1, int LDA,
385  float *B1, float *B2, int LDB, float eps)
386 {
387  float alpha = 1.0;
388  float beta = 0.0;
389  float *Residual;
390  int info_solution;
391  int maxMN = max(M, N);
392  float Rnorm, Anorm, Xnorm, Bnorm, result;
393  float *work = (float *)malloc( max(maxMN, NRHS)* sizeof(float));
394 
395  Anorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, N, A1, LDA, work);
396  Xnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), M, NRHS, B2, LDB, work);
397  Bnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), N, NRHS, B1, LDB, work);
398 
399  cblas_sgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N,
400  (alpha), A1, LDA,
401  B2, LDB,
402  (beta), B1, LDB);
403 
404  Residual = (float *)malloc(maxMN*NRHS*sizeof(float));
405  memset((void*)Residual, 0, maxMN*NRHS*sizeof(float));
406  cblas_sgemm(CblasColMajor, CblasTrans, CblasNoTrans, N, NRHS, M,
407  (alpha), A1, LDA,
408  B1, LDB,
409  (beta), Residual, maxMN);
410 
411  Rnorm = LAPACKE_slange_work(LAPACK_COL_MAJOR, lapack_const(PlasmaInfNorm), maxMN, NRHS, Residual, maxMN, work);
412  free(Residual);
413  free(work);
414 
415  if (getenv("PLASMA_TESTING_VERBOSE"))
416  printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );
417 
418  result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
419  printf("============\n");
420  printf("Checking the Residual of the solution \n");
421  printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);
422 
423  if ( isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
424  printf("-- The solution is suspicious ! \n");
425  info_solution = 1;
426  }
427  else{
428  printf("-- The solution is CORRECT ! \n");
429  info_solution = 0;
430  }
431  return info_solution;
432 }
433 #endif