testing_zher2k.c 7.34 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zher2k.c
4
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
5 6
 * @copyright 2009-2014 The University of Tennessee and The University of
 *                      Tennessee Research Foundation. All rights reserved.
Mathieu Faverge's avatar
Mathieu Faverge committed
7
 * @copyright 2012-2018 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8
 *                      Univ. Bordeaux. All rights reserved.
9
 *
10
 ***
11
 *
12
 * @brief Chameleon zher2k testing
13
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
14
 * @version 1.0.0
15
 * @comment This file has been automatically generated
16
 *          from Plasma 2.5.0 for CHAMELEON 1.0.0
17 18 19 20 21 22
 * @author Mathieu Faverge
 * @author Emmanuel Agullo
 * @author Cedric Castagnede
 * @date 2010-11-15
 * @precisions normal z -> c
 *
23
 */
24 25 26 27 28
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

29
#include <chameleon.h>
30 31 32
#include <coreblas/cblas.h>
#include <coreblas/lapacke.h>
#include <coreblas.h>
33 34
#include "testing_zauxiliary.h"

35 36 37
static int check_solution(cham_uplo_t uplo, cham_trans_t trans, int N, int K,
                          CHAMELEON_Complex64_t alpha, CHAMELEON_Complex64_t *A, int LDA,
                          CHAMELEON_Complex64_t *B, int LDB,
Mathieu Faverge's avatar
Mathieu Faverge committed
38
                          double beta,  CHAMELEON_Complex64_t *Cref, CHAMELEON_Complex64_t *Ccham, int LDC);
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56


int testing_zher2k(int argc, char **argv)
{
    int hres = 0;
    /* Check for number of arguments*/
    if ( argc != 7 ){
        USAGE("HER2K", "alpha beta M N LDA LDB LDC",
              "   - alpha : alpha coefficient\n"
              "   - beta : beta coefficient\n"
              "   - N : number of columns and rows of matrix C and number of row of matrix A and B\n"
              "   - K : number of columns of matrix A and B\n"
              "   - LDA : leading dimension of matrix A\n"
              "   - LDB : leading dimension of matrix B\n"
              "   - LDC : leading dimension of matrix C\n");
        return -1;
    }

57
    CHAMELEON_Complex64_t alpha = (CHAMELEON_Complex64_t) atol(argv[0]);
58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
    double beta  = (double) atol(argv[1]);
    int N     = atoi(argv[2]);
    int K     = atoi(argv[3]);
    int LDA   = atoi(argv[4]);
    int LDB   = atoi(argv[5]);
    int LDC   = atoi(argv[6]);
    int NKmax = max(N, K);

    double eps;
    int info_solution;
    int u, t;
    size_t LDAxK = LDA*NKmax;
    size_t LDBxK = LDB*NKmax;
    size_t LDCxN = LDC*N;

73 74 75 76 77
    CHAMELEON_Complex64_t *A      = (CHAMELEON_Complex64_t *)malloc(LDAxK*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *B      = (CHAMELEON_Complex64_t *)malloc(LDBxK*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *C      = (CHAMELEON_Complex64_t *)malloc(LDCxN*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *Cinit  = (CHAMELEON_Complex64_t *)malloc(LDCxN*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *Cfinal = (CHAMELEON_Complex64_t *)malloc(LDCxN*sizeof(CHAMELEON_Complex64_t));
78 79

    /* Check if unable to allocate memory */
80 81 82
    if ( (!A) || (!B) || (!C) || (!Cinit) || (!Cfinal) ){
        free(A); free(B); free(C);
        free(Cinit); free(Cfinal);
83 84 85 86 87 88 89
        printf("Out of Memory \n ");
        return -2;
    }

    eps = LAPACKE_dlamch_work('e');

    printf("\n");
90
    printf("------ TESTS FOR CHAMELEON ZHER2K ROUTINE -------  \n");
91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106
    printf("            Size of the Matrix C %d by %d\n", N, K);
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf("============\n");
    printf(" The relative machine precision (eps) is to be %e \n",eps);
    printf(" Computational tests pass if scaled residuals are less than 10.\n");

    /*----------------------------------------------------------
    *  TESTING ZHER2K
    */

    /* Initialize A,B */
    LAPACKE_zlarnv_work(IONE, ISEED, LDAxK, A);
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxK, B);

    /* Initialize C */
107
    CHAMELEON_zplghe( (double)0., ChamUpperLower, N, C, LDC, 51 );
108 109 110

    for (u=0; u<2; u++) {
        for (t=0; t<3; t++) {
111
            if (trans[t] == ChamTrans) continue;
112

113 114
            memcpy(Cinit,  C, LDCxN*sizeof(CHAMELEON_Complex64_t));
            memcpy(Cfinal, C, LDCxN*sizeof(CHAMELEON_Complex64_t));
115

116 117
            /* CHAMELEON ZHER2K */
            CHAMELEON_zher2k(uplo[u], trans[t], N, K, alpha, A, LDA, B, LDB, beta, Cfinal, LDC);
118 119

            /* Check the solution */
120
            info_solution = check_solution(uplo[u], trans[t], N, K,
121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
                                           alpha, A, LDA, B, LDB, beta, Cinit, Cfinal, LDC);

            if (info_solution == 0) {
                printf("***************************************************\n");
                printf(" ---- TESTING ZHER2K (%5s, %s) ........... PASSED !\n", uplostr[u], transstr[t]);
                printf("***************************************************\n");
            }
            else {
                printf("************************************************\n");
                printf(" - TESTING ZHER2K (%5s, %s) ... FAILED !\n", uplostr[u], transstr[t]);    hres++;
                printf("************************************************\n");
            }
        }
    }

    free(A); free(B); free(C);
    free(Cinit); free(Cfinal);

    return hres;
}

/*--------------------------------------------------------------
 * Check the solution
 */

146 147 148
static int check_solution(cham_uplo_t uplo, cham_trans_t trans, int N, int K,
                          CHAMELEON_Complex64_t alpha, CHAMELEON_Complex64_t *A, int LDA,
                          CHAMELEON_Complex64_t *B, int LDB,
Mathieu Faverge's avatar
Mathieu Faverge committed
149
                          double beta, CHAMELEON_Complex64_t *Cref, CHAMELEON_Complex64_t *Ccham, int LDC)
150 151
{
    int info_solution;
Mathieu Faverge's avatar
Mathieu Faverge committed
152
    double Anorm, Bnorm, Cinitnorm, Cchamnorm, Clapacknorm, Rnorm, result;
153
    double eps;
154
    CHAMELEON_Complex64_t beta_const;
155 156 157 158

    double *work = (double *)malloc(max(N, K)* sizeof(double));

    beta_const  = -1.0;
159
    Anorm       = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I',
160 161
                                      (trans == ChamNoTrans) ? N : K,
                                      (trans == ChamNoTrans) ? K : N, A, LDA, work);
162
    Bnorm       = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I',
163 164
                                      (trans == ChamNoTrans) ? N : K,
                                      (trans == ChamNoTrans) ? K : N, B, LDB, work);
165
    Cinitnorm   = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', N, N, Cref,    LDC, work);
Mathieu Faverge's avatar
Mathieu Faverge committed
166
    Cchamnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', N, N, Ccham, LDC, work);
167

168
    cblas_zher2k(CblasColMajor, (CBLAS_UPLO)uplo, (CBLAS_TRANSPOSE)trans,
169 170
                 N, K, CBLAS_SADDR(alpha), A, LDA, B, LDB, (beta), Cref, LDC);

171
    Clapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', N, N, Cref, LDC, work);
172

Mathieu Faverge's avatar
Mathieu Faverge committed
173
    cblas_zaxpy(LDC*N, CBLAS_SADDR(beta_const), Ccham, 1, Cref, 1);
174

175
    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', N, N, Cref, LDC, work);
176 177

    eps = LAPACKE_dlamch_work('e');
178

Mathieu Faverge's avatar
Mathieu Faverge committed
179 180
    printf("Rnorm %e, Anorm %e, Cinitnorm %e, Cchamnorm %e, Clapacknorm %e\n",
           Rnorm, Anorm, Cinitnorm, Cchamnorm, Clapacknorm);
181 182 183 184

    result = Rnorm / ((Anorm + Bnorm + Cinitnorm) * N * eps);
    printf("============\n");
    printf("Checking the norm of the difference against reference ZHER2K \n");
Mathieu Faverge's avatar
Mathieu Faverge committed
185
    printf("-- ||Ccham - Clapack||_oo/((||A||_oo+||C||_oo).N.eps) = %e \n", result);
186 187 188 189 190 191 192 193 194 195 196 197 198 199

    if (  isnan(Rnorm) || isinf(Rnorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
         printf("-- The solution is suspicious ! \n");
         info_solution = 1;
    }
    else {
         printf("-- The solution is CORRECT ! \n");
         info_solution= 0 ;
    }

    free(work);

    return info_solution;
}