36 USAGE(
"HEGST",
"N LDA LDB",
37 " - N : size of the matrices A and B\n"
38 " - LDA : leading dimension of the matrix A\n"
39 " - LDB : leading dimension of the matrix B\n");
43 double eps = LAPACKE_dlamch_work(
'e');
44 int N = atoi(argv[0]);
45 int LDA = atoi(argv[1]);
46 int LDB = atoi(argv[2]);
47 int info_transformation, info_factorization;
60 if ((!A1)||(!A2)||(!B1)||(!B2)||(!Ainit)||(!Binit)){
61 printf(
"Out of Memory \n ");
71 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'A', N, N, A1, LDA, Ainit, LDA);
75 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'A', N, N, B1, LDB, Binit, LDB);
78 printf(
"------ TESTS FOR PLASMA ZHEGST ROUTINE ------- \n");
79 printf(
" Size of the Matrix %d by %d\n", N, N);
81 printf(
" The matrices A and B are randomly generated for each test.\n");
82 printf(
"============\n");
83 printf(
" The relative machine precision (eps) is to be %e \n",eps);
84 printf(
" Computational tests pass if scaled residuals are less than 60.\n");
100 info_transformation = check_transformation(
itype[i],
uplo[u], N, A1, A2, LDA, B2, LDB, eps);
102 if ( (info_transformation == 0) && (info_factorization == 0) ) {
103 printf(
"***************************************************\n");
104 printf(
" ---- TESTING ZHEGST (%s, %s) ....... PASSED !\n",
itypestr[i],
uplostr[u]);
105 printf(
"***************************************************\n");
108 printf(
"************************************************\n");
109 printf(
" - TESTING ZHEGST (%s, %s) ... FAILED !\n",
itypestr[i],
uplostr[u]);
110 printf(
"************************************************\n");
130 double Anorm, Rnorm, result;
131 int info_factorization;
137 double *work = (
double *)malloc(N*
sizeof(
double));
142 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
' ', N, N, A1, LDA, Residual, N);
145 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
lapack_const(uplo), N, N, A2, LDA, L1, N);
146 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
lapack_const(uplo), N, N, A2, LDA, L2, N);
148 cblas_ztrmm(
CblasColMajor,
CblasLeft, (
CBLAS_UPLO)uplo,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), L1, N, L2, N);
150 cblas_ztrmm(
CblasColMajor,
CblasRight, (
CBLAS_UPLO)uplo,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), L1, N, L2, N);
153 for (i = 0; i < N; i++)
154 for (j = 0; j < N; j++)
155 Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];
160 result = Rnorm / ( Anorm * N * eps );
161 printf(
"============\n");
162 printf(
"Checking the Cholesky Factorization \n");
163 printf(
"-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n", result);
165 if ( isnan(result) || isinf(result) || (result > 60.0) ){
166 printf(
"-- Factorization is suspicious ! \n");
167 info_factorization = 1;
170 printf(
"-- Factorization is CORRECT ! \n");
171 info_factorization = 0;
174 free(Residual); free(L1); free(L2); free(work);
176 return info_factorization;
186 double Anorm, Rnorm, result;
187 int info_transformation;
193 double *work = (
double *)malloc(N*
sizeof(
double));
195 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'a', N, N, A1, LDA, Residual, N);
196 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
lapack_const(uplo), N, N, A2, LDA, Aorig, N);
200 for (j = 0; j < N; j++)
201 for (i = j+1; i < N; i++)
202 Aorig[j+i*N] =
conj(Aorig[i+j*N]);
204 for (i = 0; i < N; i++)
205 for (j = i+1; j < N; j++)
206 Aorig[j+i*N] =
conj(Aorig[i+j*N]);
212 cblas_ztrmm(
CblasColMajor,
CblasLeft,
CblasLower,
CblasNoTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
213 cblas_ztrmm(
CblasColMajor,
CblasRight,
CblasLower,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
217 cblas_ztrmm(
CblasColMajor,
CblasLeft,
CblasUpper,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
218 cblas_ztrmm(
CblasColMajor,
CblasRight,
CblasUpper,
CblasNoTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
223 str =
"inv(L')*A2*inv(L)";
224 cblas_ztrsm(
CblasColMajor,
CblasLeft,
CblasLower,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
225 cblas_ztrsm(
CblasColMajor,
CblasRight,
CblasLower,
CblasNoTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
228 str =
"inv(U)*A2*inv(U')";
229 cblas_ztrsm(
CblasColMajor,
CblasLeft,
CblasUpper,
CblasNoTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
230 cblas_ztrsm(
CblasColMajor,
CblasRight,
CblasUpper,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), B2, LDB, Aorig, N);
236 for (i = 0; i < N; i++)
237 for (j = 0; j < N; j++)
238 Residual[j*N+i] = Aorig[j*N+i] - Residual[j*N+i];
243 result = Rnorm / (Anorm * N *eps);
244 printf(
"============\n");
245 printf(
"Checking the global transformation \n");
246 printf(
"-- ||A1-%s||_oo/(||A1||_oo.N.eps) = %e \n", str, result );
248 if (isnan(result) || isinf(result) || (result > 60.0) ) {
249 printf(
"-- Transformation is suspicious ! \n");
250 info_transformation = 1;
253 printf(
"-- Transformation is CORRECT ! \n");
254 info_transformation = 0;
257 free(Residual); free(Aorig); free(work);
259 return info_transformation;