26 int minMN =
min(M, N);
28 double *work = (
double *)malloc(minMN*
sizeof(
double));
30 eps = LAPACKE_dlamch_work(
'e');
37 for (i = 0; i < minMN; i++)
42 cblas_zherk(
CblasColMajor,
CblasUpper,
CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N);
44 cblas_zherk(
CblasColMajor,
CblasUpper,
CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);
46 normQ = LAPACKE_zlansy_work(LAPACK_COL_MAJOR,
'i',
'u', minMN, Id, minMN, work);
48 printf(
"============\n");
49 printf(
"Checking the orthogonality of Q \n");
50 printf(
"||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));
52 if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
53 printf(
"-- Orthogonality is suspicious ! \n");
57 printf(
"-- Orthogonality is CORRECT ! \n");
74 int info_factorization;
78 eps = LAPACKE_dlamch_work(
'e');
82 double *work = (
double *)malloc(
max(M,N)*
sizeof(double));
91 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'u', M, N, A2, LDA, R, N);
95 cblas_zgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, N, N,
CBLAS_SADDR(alpha), Q, LDA, R, N,
CBLAS_SADDR(beta), Ql, M);
102 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'l', M, N, A2, LDA, L, M);
106 cblas_zgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, N, M,
CBLAS_SADDR(alpha), L, M, Q, LDA,
CBLAS_SADDR(beta), Ql, M);
111 for (i = 0; i < M; i++)
112 for (j = 0 ; j < N; j++)
113 Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];
115 Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, Residual, M, work);
116 Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, A2, LDA, work);
119 printf(
"============\n");
120 printf(
"Checking the QR Factorization \n");
121 printf(
"-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
124 printf(
"============\n");
125 printf(
"Checking the LQ Factorization \n");
126 printf(
"-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
129 if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
130 printf(
"-- Factorization is suspicious ! \n");
131 info_factorization = 1;
134 printf(
"-- Factorization is CORRECT ! \n");
135 info_factorization = 0;
138 free(work); free(Ql); free(Residual);
140 return info_factorization;
151 int info_factorization;
155 eps = LAPACKE_dlamch_work(
'e');
160 double *work = (
double *)malloc(N*
sizeof(
double));
167 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
' ', N, N, A1, LDA, Residual, N);
171 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'u', N, N, A2, LDA, L1, N);
172 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'u', N, N, A2, LDA, L2, N);
173 cblas_ztrmm(
CblasColMajor,
CblasLeft,
CblasUpper,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), L1, N, L2, N);
176 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'l', N, N, A2, LDA, L1, N);
177 LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,
'l', N, N, A2, LDA, L2, N);
178 cblas_ztrmm(
CblasColMajor,
CblasRight,
CblasLower,
CblasConjTrans,
CblasNonUnit, N, N,
CBLAS_SADDR(alpha), L1, N, L2, N);
182 for (i = 0; i < N; i++)
183 for (j = 0; j < N; j++)
184 Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];
186 Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', N, N, Residual, N, work);
187 Anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', N, N, A1, LDA, work);
189 printf(
"============\n");
190 printf(
"Checking the Cholesky Factorization \n");
191 printf(
"-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
193 if ( isnan(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 10.0) ){
194 printf(
"-- Factorization is suspicious ! \n");
195 info_factorization = 1;
198 printf(
"-- Factorization is CORRECT ! \n");
199 info_factorization = 0;
202 free(Residual); free(L1); free(L2); free(work);
204 return info_factorization;
215 double *Cinitnorm,
double *Cplasmanorm,
double *Clapacknorm )
219 double *work = (
double *)malloc(
max(K,
max(M, N))*
sizeof(double));
221 *Cinitnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, Cref, LDC, work);
222 *Cplasmanorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, Cplasma, LDC, work);
227 *Clapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, Cref, LDC, work);
231 Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, Cref, LDC, work);
245 double *Binitnorm,
double *Bplasmanorm,
double *Blapacknorm )
249 double *work = (
double *)malloc(
max(M, NRHS)*
sizeof(double));
252 *Binitnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, NRHS, Bref, LDB, work);
253 *Bplasmanorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'm', M, NRHS, Bplasma, LDB, work);
259 *Blapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'm', M, NRHS, Bref, LDB, work);
263 Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'm', M, NRHS, Bref, LDB, work);
264 Rnorm = Rnorm / *Blapacknorm;
278 double *anorm,
double *bnorm,
double *xnorm )
281 double Rnorm = -1.00;
284 double *work = (
double *)malloc(
max(M, N)*
sizeof(double));
286 *anorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, N, A, LDA, work);
287 *xnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', M, NRHS, X, LDB, work);
288 *bnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', N, NRHS, B, LDB, work);
290 cblas_zgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, NRHS, N,
CBLAS_SADDR(zone), A, LDA, X, LDB,
CBLAS_SADDR(mzone), B, LDB);
292 Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR,
'i', N, NRHS, B, M, work);