28 #define max(a, b) ((a) > (b) ? (a) : (b))
31 #define min(a, b) ((a) < (b) ? (a) : (b))
49 int info_ortho, info_factorization;
53 float *A1 = (
float *)malloc(LDA*N*
sizeof(
float));
54 float *A2 = (
float *)malloc(LDA*N*
sizeof(
float));
55 float *
Q = (
float *)malloc(LDA*N*
sizeof(
float));
59 if ((!A1)||(!A2)||(!Q)){
60 printf(
"Out of Memory \n ");
66 printf(
"-- PLASMA is initialized to run on %d cores. \n",cores);
72 LAPACKE_slarnv_work(
IONE,
ISEED, LDAxN, A1);
73 for (i = 0; i < M; i++)
74 for (j = 0; j < N; j++)
75 A2[LDA*j+i] = A1[LDA*j+i] ;
81 memset((
void*)Q, 0, LDA*N*
sizeof(
float));
82 for (i = 0; i < K; i++)
91 if ((info_ortho != 0)|(info_factorization != 0)|(info != 0))
92 printf(
"-- Error in SGELQF example ! \n");
94 printf(
"-- Run of SGELQF example successful ! \n");
96 free(A1); free(A2); free(Q); free(T);
113 int minMN =
min(M, N);
115 float *work = (
float *)malloc(minMN*
sizeof(
float));
117 eps = LAPACKE_slamch_work(
'e');
122 float *Id = (
float *) malloc(minMN*minMN*
sizeof(
float));
123 memset((
void*)Id, 0, minMN*minMN*
sizeof(
float));
124 for (i = 0; i < minMN; i++)
125 Id[i*minMN+i] = (
float)1.0;
129 cblas_ssyrk(
CblasColMajor,
CblasUpper,
CblasTrans, N, M, alpha, Q, LDQ, beta, Id, N);
131 cblas_ssyrk(
CblasColMajor,
CblasUpper,
CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);
135 printf(
"============\n");
136 printf(
"Checking the orthogonality of Q \n");
137 printf(
"||Id-Q'*Q||_oo / (N*eps) = %e \n",normQ/(minMN*eps));
139 if ( isnan(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 10.0) ) {
140 printf(
"-- Orthogonality is suspicious ! \n");
144 printf(
"-- Orthogonality is CORRECT ! \n");
148 free(work); free(Id);
161 int info_factorization;
165 eps = LAPACKE_slamch_work(
'e');
167 float *Ql = (
float *)malloc(M*N*
sizeof(
float));
168 float *Residual = (
float *)malloc(M*N*
sizeof(
float));
169 float *work = (
float *)malloc(
max(M,N)*
sizeof(float));
176 float *R = (
float *)malloc(N*N*
sizeof(
float));
177 memset((
void*)R, 0, N*N*
sizeof(
float));
178 LAPACKE_slacpy_work(LAPACK_COL_MAJOR,
'u', M, N, A2, LDA, R, N);
181 memset((
void*)Ql, 0, M*N*
sizeof(
float));
182 cblas_sgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, N, N, (alpha), Q, LDA, R, N, (beta), Ql, M);
187 float *
L = (
float *)malloc(M*M*
sizeof(
float));
188 memset((
void*)L, 0, M*M*
sizeof(
float));
189 LAPACKE_slacpy_work(LAPACK_COL_MAJOR,
'l', M, N, A2, LDA, L, M);
192 memset((
void*)Ql, 0, M*N*
sizeof(
float));
193 cblas_sgemm(
CblasColMajor,
CblasNoTrans,
CblasNoTrans, M, N, M, (alpha), L, M, Q, LDA, (beta), Ql, M);
198 for (i = 0; i < M; i++)
199 for (j = 0 ; j < N; j++)
200 Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];
206 printf(
"============\n");
207 printf(
"Checking the QR Factorization \n");
208 printf(
"-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
211 printf(
"============\n");
212 printf(
"Checking the LQ Factorization \n");
213 printf(
"-- ||A-LQ||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
216 if (isnan(Rnorm / (Anorm * N *eps)) || (Rnorm / (Anorm * N * eps) > 10.0) ) {
217 printf(
"-- Factorization is suspicious ! \n");
218 info_factorization = 1;
221 printf(
"-- Factorization is CORRECT ! \n");
222 info_factorization = 0;
225 free(work); free(Ql); free(Residual);
227 return info_factorization;