diff --git a/testing/testing_zgels.c b/testing/testing_zgels.c index 1259502af39c374e83d1e28936c80cf804b9cb0f..46da823c103411ec88f9c4afabf1b4ae9228cae0 100644 --- a/testing/testing_zgels.c +++ b/testing/testing_zgels.c @@ -34,170 +34,6 @@ #include <coreblas.h> #include "testing_zauxiliary.h" -enum blas_order_type { - blas_rowmajor = 101, - blas_colmajor = 102 }; - -enum blas_uplo_type { - blas_upper = 121, - blas_lower = 122 }; - -enum blas_cmach_type { - blas_base = 151, - blas_t = 152, - blas_rnd = 153, - blas_ieee = 154, - blas_emin = 155, - blas_emax = 156, - blas_eps = 157, - blas_prec = 158, - blas_underflow = 159, - blas_overflow = 160, - blas_sfmin = 161}; - -enum blas_norm_type { - blas_one_norm = 171, - blas_real_one_norm = 172, - blas_two_norm = 173, - blas_frobenius_norm = 174, - blas_inf_norm = 175, - blas_real_inf_norm = 176, - blas_max_norm = 177, - blas_real_max_norm = 178 }; - -static void -BLAS_error(char *rname, int err, int val, int x) { - fprintf( stderr, "%s %d %d %d\n", rname, err, val, x ); - abort(); -} - -static -void -BLAS_zsy_norm(enum blas_order_type order, enum blas_norm_type norm, - enum blas_uplo_type uplo, int n, const MORSE_Complex64_t *a, int lda, double *res) { - int i, j; double anorm, v; - char rname[] = "BLAS_zsy_norm"; - - if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 ); - - if (norm == blas_inf_norm) { - anorm = 0.0; - if (blas_upper == uplo) { - for (i = 0; i < n; ++i) { - v = 0.0; - for (j = 0; j < i; ++j) { - v += cabs( a[j + i * lda] ); - } - for (j = i; j < n; ++j) { - v += cabs( a[i + j * lda] ); - } - if (v > anorm) - anorm = v; - } - } else { - BLAS_error( rname, -3, norm, 0 ); - return; - } - } else { - BLAS_error( rname, -2, norm, 0 ); - return; - } - - if (res) *res = anorm; -} - -static -void -BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm, - int m, int n, const MORSE_Complex64_t *a, int lda, double *res) { - int i, j; float anorm, v; - char rname[] = "BLAS_zge_norm"; - - if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 ); - - if (norm == blas_frobenius_norm) { - anorm = 0.0f; - for (j = n; j; --j) { - for (i = m; i; --i) { - v = a[0]; - anorm += v * v; - a++; - } - a += lda - m; - } - anorm = sqrt( anorm ); - } else if (norm == blas_inf_norm) { - anorm = 0.0f; - for (i = 0; i < m; ++i) { - v = 0.0f; - for (j = 0; j < n; ++j) { - v += cabs( a[i + j * lda] ); - } - if (v > anorm) - anorm = v; - } - } else { - BLAS_error( rname, -2, norm, 0 ); - return; - } - - if (res) *res = anorm; -} - -static -double -BLAS_dpow_di(double x, int n) { - double rv = 1.0; - - if (n < 0) { - n = -n; - x = 1.0 / x; - } - - for (; n; n >>= 1, x *= x) { - if (n & 1) - rv *= x; - } - - return rv; -} - -static -double -BLAS_dfpinfo(enum blas_cmach_type cmach) { - double eps = 1.0, r = 1.0, o = 1.0, b = 2.0; - int t = 53, l = 1024, m = -1021; - char rname[] = "BLAS_dfpinfo"; - - if ((sizeof eps) == sizeof(float)) { - t = 24; - l = 128; - m = -125; - } else { - t = 53; - l = 1024; - m = -1021; - } - - /* for (i = 0; i < t; ++i) eps *= half; */ - eps = BLAS_dpow_di( b, -t ); - /* for (i = 0; i >= m; --i) r *= half; */ - r = BLAS_dpow_di( b, m-1 ); - - o -= eps; - /* for (i = 0; i < l; ++i) o *= b; */ - o = (o * BLAS_dpow_di( b, l-1 )) * b; - - switch (cmach) { - case blas_eps: return eps; - case blas_sfmin: return r; - default: - BLAS_error( rname, -1, cmach, 0 ); - break; - } - return 0.0; -} - static int check_orthogonality(int, int, int, MORSE_Complex64_t*, double); static int check_factorization(int, int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, MORSE_Complex64_t*, double); static int check_solution(int, int, int, MORSE_Complex64_t*, int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, double); @@ -268,7 +104,7 @@ int testing_zgels(int argc, char **argv) MORSE_Alloc_Workspace_zgels(M, N, &T, 1, 1); memset(T->mat, 0, (T->llm*T->lln)*sizeof(MORSE_Complex64_t)); - eps = BLAS_dfpinfo( blas_eps ); + eps = LAPACKE_dlamch_work('e'); /*---------------------------------------------------------- * TESTING ZGELS @@ -508,16 +344,18 @@ static int check_orthogonality(int M, int N, int LDQ, MORSE_Complex64_t *Q, doub /* Build the idendity matrix USE DLASET?*/ MORSE_Complex64_t *Id = (MORSE_Complex64_t *) malloc(minMN*minMN*sizeof(MORSE_Complex64_t)); memset((void*)Id, 0, minMN*minMN*sizeof(MORSE_Complex64_t)); - for (i = 0; i < minMN; i++) + for (i = 0; i < minMN; i++) { Id[i*minMN+i] = (MORSE_Complex64_t)1.0; + } /* Perform Id - Q'Q */ - if (M >= N) + if (M >= N) { cblas_zherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N); - else + } + else { cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M); - - BLAS_zsy_norm( blas_colmajor, blas_inf_norm, blas_upper, minMN, Id, minMN, &normQ ); + } + normQ = LAPACKE_zlansy_work( LAPACK_COL_MAJOR, 'i', 'u', minMN, Id, minMN, work ); printf("============\n"); printf("Checking the orthogonality of Q \n"); @@ -583,8 +421,8 @@ static int check_factorization(int M, int N, MORSE_Complex64_t *A1, MORSE_Comple for (j = 0 ; j < N; j++) Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i]; - BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, Residual, M, &Rnorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, A2, LDA, &Anorm ); + Rnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', M, N, Residual, M, work ); + Anorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', M, N, A2, LDA, work ); if (M >= N) { printf("============\n"); @@ -626,9 +464,9 @@ static int check_solution(int M, int N, int NRHS, MORSE_Complex64_t *A, int LDA, alpha = 1.0; beta = -1.0; - BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, A, LDA, &Anorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B, LDB, &Bnorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, NRHS, X, LDB, &Xnorm ); + Anorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', M, N, A, LDA, work ); + Bnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', N, NRHS, B, LDB, work ); + Xnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', M, NRHS, X, LDB, work ); cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, CBLAS_SADDR(alpha), A, LDA, X, LDB, CBLAS_SADDR(beta), B, LDB); @@ -636,14 +474,14 @@ static int check_solution(int M, int N, int NRHS, MORSE_Complex64_t *A, int LDA, MORSE_Complex64_t *Residual = (MORSE_Complex64_t *)malloc(M*NRHS*sizeof(MORSE_Complex64_t)); memset((void*)Residual, 0, M*NRHS*sizeof(MORSE_Complex64_t)); cblas_zgemm(CblasColMajor, CblasConjTrans, CblasNoTrans, N, NRHS, M, CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Residual, M); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, NRHS, Residual, M, &Rnorm ); + Rnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', M, NRHS, Residual, M, work ); free(Residual); } else { MORSE_Complex64_t *Residual = (MORSE_Complex64_t *)malloc(N*NRHS*sizeof(MORSE_Complex64_t)); memset((void*)Residual, 0, N*NRHS*sizeof(MORSE_Complex64_t)); cblas_zgemm(CblasColMajor, CblasConjTrans, CblasNoTrans, N, NRHS, M, CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Residual, N); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, Residual, N, &Rnorm ); + Rnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', N, NRHS, Residual, N, work ); free(Residual); } diff --git a/testing/testing_zgesv_incpiv.c b/testing/testing_zgesv_incpiv.c index dcc2a25d2687bae71baff909b0f192ad9930fe4e..b7702f153efe5c24c8a0f702e0ed8aa22e66e3a1 100644 --- a/testing/testing_zgesv_incpiv.c +++ b/testing/testing_zgesv_incpiv.c @@ -33,131 +33,6 @@ #include <coreblas.h> #include "testing_zauxiliary.h" -enum blas_order_type { - blas_rowmajor = 101, - blas_colmajor = 102 }; - -enum blas_cmach_type { - blas_base = 151, - blas_t = 152, - blas_rnd = 153, - blas_ieee = 154, - blas_emin = 155, - blas_emax = 156, - blas_eps = 157, - blas_prec = 158, - blas_underflow = 159, - blas_overflow = 160, - blas_sfmin = 161}; - -enum blas_norm_type { - blas_one_norm = 171, - blas_real_one_norm = 172, - blas_two_norm = 173, - blas_frobenius_norm = 174, - blas_inf_norm = 175, - blas_real_inf_norm = 176, - blas_max_norm = 177, - blas_real_max_norm = 178 }; - -static void -BLAS_error(char *rname, int err, int val, int x) { - fprintf( stderr, "%s %d %d %d\n", rname, err, val, x ); - abort(); -} - -static -void -BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm, - int m, int n, const MORSE_Complex64_t *a, int lda, double *res) { - int i, j; float anorm, v; - char rname[] = "BLAS_zge_norm"; - - if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 ); - - if (norm == blas_frobenius_norm) { - anorm = 0.0f; - for (j = n; j; --j) { - for (i = m; i; --i) { - v = a[0]; - anorm += v * v; - a++; - } - a += lda - m; - } - anorm = sqrt( anorm ); - } else if (norm == blas_inf_norm) { - anorm = 0.0f; - for (i = 0; i < m; ++i) { - v = 0.0f; - for (j = 0; j < n; ++j) { - v += cabs( a[i + j * lda] ); - } - if (v > anorm) - anorm = v; - } - } else { - BLAS_error( rname, -2, norm, 0 ); - return; - } - - if (res) *res = anorm; -} - -static -double -BLAS_dpow_di(double x, int n) { - double rv = 1.0; - - if (n < 0) { - n = -n; - x = 1.0 / x; - } - - for (; n; n >>= 1, x *= x) { - if (n & 1) - rv *= x; - } - - return rv; -} - -static -double -BLAS_dfpinfo(enum blas_cmach_type cmach) { - double eps = 1.0, r = 1.0, o = 1.0, b = 2.0; - int t = 53, l = 1024, m = -1021; - char rname[] = "BLAS_dfpinfo"; - - if ((sizeof eps) == sizeof(float)) { - t = 24; - l = 128; - m = -125; - } else { - t = 53; - l = 1024; - m = -1021; - } - - /* for (i = 0; i < t; ++i) eps *= half; */ - eps = BLAS_dpow_di( b, -t ); - /* for (i = 0; i >= m; --i) r *= half; */ - r = BLAS_dpow_di( b, m-1 ); - - o -= eps; - /* for (i = 0; i < l; ++i) o *= b; */ - o = (o * BLAS_dpow_di( b, l-1 )) * b; - - switch (cmach) { - case blas_eps: return eps; - case blas_sfmin: return r; - default: - BLAS_error( rname, -1, cmach, 0 ); - break; - } - return 0.0; -} - static int check_solution(int, int , MORSE_Complex64_t *, int, MORSE_Complex64_t *, MORSE_Complex64_t *, int, double); int testing_zgesv_incpiv(int argc, char **argv) @@ -199,7 +74,7 @@ int testing_zgesv_incpiv(int argc, char **argv) return -2; } - eps = BLAS_dfpinfo(blas_eps); + eps = LAPACKE_dlamch_work( 'e' ); /*---------------------------------------------------------- * TESTING ZGESV @@ -352,12 +227,12 @@ static int check_solution(int N, int NRHS, MORSE_Complex64_t *A1, int LDA, MORSE alpha = 1.0; beta = -1.0; - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B2, LDB, &Xnorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Bnorm ); + Xnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B2, LDB ); + Anorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N, A1, LDA ); + Bnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B1, LDB ); cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, CBLAS_SADDR(alpha), A1, LDA, B2, LDB, CBLAS_SADDR(beta), B1, LDB); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Rnorm ); + Rnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B1, LDB ); if (getenv("MORSE_TESTING_VERBOSE")) printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm ); diff --git a/testing/testing_zposv.c b/testing/testing_zposv.c index 82e0d5e26ffa7eade5aacb343e587078c3a66d56..f2e1ead1759acce639fef1d56e95310263d542f0 100644 --- a/testing/testing_zposv.c +++ b/testing/testing_zposv.c @@ -33,131 +33,6 @@ #include <coreblas.h> #include "testing_zauxiliary.h" -enum blas_order_type { - blas_rowmajor = 101, - blas_colmajor = 102 }; - -enum blas_cmach_type { - blas_base = 151, - blas_t = 152, - blas_rnd = 153, - blas_ieee = 154, - blas_emin = 155, - blas_emax = 156, - blas_eps = 157, - blas_prec = 158, - blas_underflow = 159, - blas_overflow = 160, - blas_sfmin = 161}; - -enum blas_norm_type { - blas_one_norm = 171, - blas_real_one_norm = 172, - blas_two_norm = 173, - blas_frobenius_norm = 174, - blas_inf_norm = 175, - blas_real_inf_norm = 176, - blas_max_norm = 177, - blas_real_max_norm = 178 }; - -static void -BLAS_error(char *rname, int err, int val, int x) { - fprintf( stderr, "%s %d %d %d\n", rname, err, val, x ); - abort(); -} - -static -void -BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm, - int m, int n, const MORSE_Complex64_t *a, int lda, double *res) { - int i, j; float anorm, v; - char rname[] = "BLAS_zge_norm"; - - if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 ); - - if (norm == blas_frobenius_norm) { - anorm = 0.0f; - for (j = n; j; --j) { - for (i = m; i; --i) { - v = a[0]; - anorm += v * v; - a++; - } - a += lda - m; - } - anorm = sqrt( anorm ); - } else if (norm == blas_inf_norm) { - anorm = 0.0f; - for (i = 0; i < m; ++i) { - v = 0.0f; - for (j = 0; j < n; ++j) { - v += cabs( a[i + j * lda] ); - } - if (v > anorm) - anorm = v; - } - } else { - BLAS_error( rname, -2, norm, 0 ); - return; - } - - if (res) *res = anorm; -} - -static -double -BLAS_dpow_di(double x, int n) { - double rv = 1.0; - - if (n < 0) { - n = -n; - x = 1.0 / x; - } - - for (; n; n >>= 1, x *= x) { - if (n & 1) - rv *= x; - } - - return rv; -} - -static -double -BLAS_dfpinfo(enum blas_cmach_type cmach) { - double eps = 1.0, r = 1.0, o = 1.0, b = 2.0; - int t = 53, l = 1024, m = -1021; - char rname[] = "BLAS_dfpinfo"; - - if ((sizeof eps) == sizeof(float)) { - t = 24; - l = 128; - m = -125; - } else { - t = 53; - l = 1024; - m = -1021; - } - - /* for (i = 0; i < t; ++i) eps *= half; */ - eps = BLAS_dpow_di( b, -t ); - /* for (i = 0; i >= m; --i) r *= half; */ - r = BLAS_dpow_di( b, m-1 ); - - o -= eps; - /* for (i = 0; i < l; ++i) o *= b; */ - o = (o * BLAS_dpow_di( b, l-1 )) * b; - - switch (cmach) { - case blas_eps: return eps; - case blas_sfmin: return r; - default: - BLAS_error( rname, -1, cmach, 0 ); - break; - } - return 0.0; -} - static int check_factorization(int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, int , double); static int check_solution(int, int, MORSE_Complex64_t*, int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, double); @@ -198,7 +73,7 @@ int testing_zposv(int argc, char **argv) return -2; } - eps = BLAS_dfpinfo( blas_eps ); + eps = LAPACKE_dlamch_work( 'e' ); uplo = MorseUpper; trans1 = uplo == MorseUpper ? MorseConjTrans : MorseNoTrans; @@ -371,8 +246,8 @@ static int check_factorization(int N, MORSE_Complex64_t *A1, MORSE_Complex64_t * for (j = 0; j < N; j++) Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i]; - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, Residual, N, &Rnorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm ); + Rnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N, Residual, N ); + Anorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N, A1, LDA ); printf("============\n"); printf("Checking the Cholesky Factorization \n"); @@ -407,12 +282,12 @@ static int check_solution(int N, int NRHS, MORSE_Complex64_t *A1, int LDA, MORSE alpha = 1.0; beta = -1.0; - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B2, LDB, &Xnorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Bnorm ); + Xnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B2, LDB ); + Anorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N, A1, LDA ); + Bnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B1, LDB ); cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, CBLAS_SADDR(alpha), A1, LDA, B2, LDB, CBLAS_SADDR(beta), B1, LDB); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B1, LDB, &Rnorm ); + Rnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B1, LDB ); if (getenv("MORSE_TESTING_VERBOSE")) printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm ); diff --git a/testing/testing_zpotri.c b/testing/testing_zpotri.c index f16dccb77b569e3f9d6e8db2508dd685d0f3e526..379db6631896321a20e1dc95e47cdee1c8d003f6 100644 --- a/testing/testing_zpotri.c +++ b/testing/testing_zpotri.c @@ -33,141 +33,6 @@ #include <coreblas.h> #include "testing_zauxiliary.h" -enum blas_order_type { - blas_rowmajor = 101, - blas_colmajor = 102 }; - -enum blas_cmach_type { - blas_base = 151, - blas_t = 152, - blas_rnd = 153, - blas_ieee = 154, - blas_emin = 155, - blas_emax = 156, - blas_eps = 157, - blas_prec = 158, - blas_underflow = 159, - blas_overflow = 160, - blas_sfmin = 161}; - -enum blas_norm_type { - blas_one_norm = 171, - blas_real_one_norm = 172, - blas_two_norm = 173, - blas_frobenius_norm = 174, - blas_inf_norm = 175, - blas_real_inf_norm = 176, - blas_max_norm = 177, - blas_real_max_norm = 178 }; - -static void -BLAS_error(char *rname, int err, int val, int x) { - fprintf( stderr, "%s %d %d %d\n", rname, err, val, x ); - abort(); -} - -static -void -BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm, - int m, int n, const MORSE_Complex64_t *a, int lda, double *res) { - int i, j; float anorm, v; - char rname[] = "BLAS_zge_norm"; - - if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 ); - - if (norm == blas_frobenius_norm) { - anorm = 0.0f; - for (j = n; j; --j) { - for (i = m; i; --i) { - v = a[0]; - anorm += v * v; - a++; - } - a += lda - m; - } - anorm = sqrt( anorm ); - } else if (norm == blas_inf_norm) { - anorm = 0.0f; - for (i = 0; i < m; ++i) { - v = 0.0f; - for (j = 0; j < n; ++j) { - v += cabs( a[i + j * lda] ); - } - if (v > anorm) - anorm = v; - } - } else if (norm == blas_one_norm) { - anorm = 0.0f; - for (i = 0; i < m; ++i) { - v = 0.0f; - for (j = 0; j < n; ++j) { - v += cabs( a[i + j * lda] ); - } - if (v > anorm) - anorm = v; - } - } else { - BLAS_error( rname, -2, norm, 0 ); - return; - } - - if (res) *res = anorm; -} - -static -double -BLAS_dpow_di(double x, int n) { - double rv = 1.0; - - if (n < 0) { - n = -n; - x = 1.0 / x; - } - - for (; n; n >>= 1, x *= x) { - if (n & 1) - rv *= x; - } - - return rv; -} - -static -double -BLAS_dfpinfo(enum blas_cmach_type cmach) { - double eps = 1.0, r = 1.0, o = 1.0, b = 2.0; - int t = 53, l = 1024, m = -1021; - char rname[] = "BLAS_dfpinfo"; - - if ((sizeof eps) == sizeof(float)) { - t = 24; - l = 128; - m = -125; - } else { - t = 53; - l = 1024; - m = -1021; - } - - /* for (i = 0; i < t; ++i) eps *= half; */ - eps = BLAS_dpow_di( b, -t ); - /* for (i = 0; i >= m; --i) r *= half; */ - r = BLAS_dpow_di( b, m-1 ); - - o -= eps; - /* for (i = 0; i < l; ++i) o *= b; */ - o = (o * BLAS_dpow_di( b, l-1 )) * b; - - switch (cmach) { - case blas_eps: return eps; - case blas_sfmin: return r; - default: - BLAS_error( rname, -1, cmach, 0 ); - break; - } - return 0.0; -} - static int check_factorization(int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, int, double); static int check_inverse(int, MORSE_Complex64_t *, MORSE_Complex64_t *, int, int, double); @@ -202,13 +67,13 @@ int testing_zpotri(int argc, char **argv) return -2; } - eps = BLAS_dfpinfo( blas_eps ); + eps = LAPACKE_dlamch_work( 'e' ); uplo = MorseUpper; /*------------------------------------------------------------- - * TESTING ZPOTRI - */ + * TESTING ZPOTRI + */ /* Initialize A1 and A2 for Symmetric Positif Matrix */ MORSE_zplghe( (double)N, MorseUpperLower, N, A1, LDA, 51 ); @@ -289,10 +154,10 @@ static int check_factorization(int N, MORSE_Complex64_t *A1, MORSE_Complex64_t * /* Compute the Residual || A -L'L|| */ for (i = 0; i < N; i++) for (j = 0; j < N; j++) - Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i]; + Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i]; - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, Residual, N, &Rnorm ); - BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, N, A1, LDA, &Anorm ); + Rnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', N, N, Residual, N, work ); + Anorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'i', N, N, A1, LDA, work ); printf("============\n"); printf("Checking the Cholesky Factorization \n"); @@ -331,29 +196,34 @@ static int check_inverse(int N, MORSE_Complex64_t *A1, MORSE_Complex64_t *A2, in /* Rebuild the other part of the inverse matrix */ if(uplo == MorseUpper){ - for(j=0; j<N; j++) - for(i=0; i<j; i++) - *(A2+j+i*LDA) = *(A2+i+j*LDA); - cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, N, N, CBLAS_SADDR(alpha), A2, LDA, A1, LDA, CBLAS_SADDR(beta), work, N); + for(j=0; j<N; j++) { + for(i=0; i<j; i++) { + *(A2+j+i*LDA) = *(A2+i+j*LDA); + } + } + cblas_zhemm(CblasColMajor, CblasLeft, CblasUpper, N, N, CBLAS_SADDR(alpha), A2, LDA, A1, LDA, CBLAS_SADDR(beta), work, N); } else { - for(j=0; j<N; j++) - for(i=j; i<N; i++) - *(A2+j+i*LDA) = *(A2+i+j*LDA); - cblas_zhemm(CblasColMajor, CblasLeft, CblasLower, N, N, CBLAS_SADDR(alpha), A2, LDA, A1, LDA, CBLAS_SADDR(beta), work, N); + for(j=0; j<N; j++) { + for(i=j; i<N; i++) { + *(A2+j+i*LDA) = *(A2+i+j*LDA); + } + } + cblas_zhemm(CblasColMajor, CblasLeft, CblasLower, N, N, CBLAS_SADDR(alpha), A2, LDA, A1, LDA, CBLAS_SADDR(beta), work, N); } /* Add the identity matrix to work */ for(i=0; i<N; i++) *(work+i+i*N) = *(work+i+i*N) + zone; - BLAS_zge_norm( blas_colmajor, blas_one_norm, N, N, work, N, &Rnorm ); - BLAS_zge_norm( blas_colmajor, blas_one_norm, N, N, A1, LDA, &Anorm ); - BLAS_zge_norm( blas_colmajor, blas_one_norm, N, N, A2, LDA, &Ainvnorm ); + Rnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'o', N, N, work, N ); + Anorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'o', N, N, A1, LDA ); + Ainvnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'o', N, N, A2, LDA ); - if (getenv("MORSE_TESTING_VERBOSE")) - printf( "||A||_1=%f\n||Ainv||_1=%f\n||Id - A*Ainv||_1=%e\n", Anorm, Ainvnorm, Rnorm ); + if (getenv("MORSE_TESTING_VERBOSE")) { + printf( "||A||_1=%f\n||Ainv||_1=%f\n||Id - A*Ainv||_1=%e\n", Anorm, Ainvnorm, Rnorm ); + } result = Rnorm / ( (Anorm*Ainvnorm)*N*eps ) ; printf("============\n"); @@ -363,7 +233,7 @@ static int check_inverse(int N, MORSE_Complex64_t *A1, MORSE_Complex64_t *A2, in if ( isnan(Ainvnorm) || isinf(Ainvnorm) || isnan(result) || isinf(result) || (result > 60.0) ) { printf("-- The inverse is suspicious ! \n"); info_inverse = 1; - } + } else{ printf("-- The inverse is CORRECT ! \n"); info_inverse = 0;