testing_zhemm.c 7.28 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zhemm.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.
PRUVOST Florent's avatar
PRUVOST Florent committed
7
 * @copyright 2012-2019 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
8
 *                      Univ. Bordeaux. All rights reserved.
9
 *
10
 ***
11
 *
12
 * @brief Chameleon zhemm testing
13
 *
PRUVOST Florent's avatar
PRUVOST Florent committed
14
 * @version 0.9.2
15
 * @comment This file has been automatically generated
PRUVOST Florent's avatar
PRUVOST Florent committed
16
 *          from Plasma 2.5.0 for CHAMELEON 0.9.2
17 18 19
 * @author Mathieu Faverge
 * @author Emmanuel Agullo
 * @author Cedric Castagnede
PRUVOST Florent's avatar
PRUVOST Florent committed
20
 * @date 2014-11-16
21 22
 * @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_side_t side, cham_uplo_t uplo, int M, int N,
                          CHAMELEON_Complex64_t alpha, CHAMELEON_Complex64_t *A, int LDA,
                          CHAMELEON_Complex64_t *B, int LDB,
Mathieu Faverge's avatar
Mathieu Faverge committed
38
                          CHAMELEON_Complex64_t 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

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

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

    double eps;
    int info_solution;
    int i, j, s, u;
    int LDAxM = LDA*MNmax;
    int LDBxN = LDB*N;
    int LDCxN = LDC*N;

72 73 74 75 76
    CHAMELEON_Complex64_t *A      = (CHAMELEON_Complex64_t *)malloc(LDAxM*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *B      = (CHAMELEON_Complex64_t *)malloc(LDBxN*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));
77 78

    /* Check if unable to allocate memory */
79 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 ZHEMM ROUTINE -------  \n");
91 92 93 94 95 96 97 98 99 100 101 102
    printf("            Size of the Matrix %d by %d\n", M, N);
    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 ZHEMM
    */

    /* Initialize A */
103
    CHAMELEON_zplghe( (double)0., ChamUpperLower, MNmax, A, LDA, 51 );
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121

    /* Initialize B */
    LAPACKE_zlarnv_work(IONE, ISEED, LDBxN, B);

    /* Initialize C */
    LAPACKE_zlarnv_work(IONE, ISEED, LDCxN, C);

    for (s=0; s<2; s++) {
        for (u=0; u<2; u++) {

            /* Initialize  Cinit / Cfinal */
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cinit[LDC*j+i] = C[LDC*j+i];
            for ( i = 0; i < M; i++)
                for (  j = 0; j < N; j++)
                    Cfinal[LDC*j+i] = C[LDC*j+i];

122 123
            /* CHAMELEON ZHEMM */
            CHAMELEON_zhemm(side[s], uplo[u], M, N, alpha, A, LDA, B, LDB, beta, Cfinal, LDC);
124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150

            /* Check the solution */
            info_solution = check_solution(side[s], uplo[u], M, N, alpha, A, LDA, B, LDB, beta, Cinit, Cfinal, LDC);

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

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

    return hres;
}

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

151 152 153
static int check_solution(cham_side_t side, cham_uplo_t uplo, int M, int N,
                   CHAMELEON_Complex64_t alpha, CHAMELEON_Complex64_t *A, int LDA,
                   CHAMELEON_Complex64_t *B, int LDB,
Mathieu Faverge's avatar
Mathieu Faverge committed
154
                   CHAMELEON_Complex64_t beta, CHAMELEON_Complex64_t *Cref, CHAMELEON_Complex64_t *Ccham, int LDC)
155 156
{
    int info_solution, NrowA;
Mathieu Faverge's avatar
Mathieu Faverge committed
157
    double Anorm, Bnorm, Cinitnorm, Cchamnorm, Clapacknorm, Rnorm;
158
    double eps;
159
    CHAMELEON_Complex64_t beta_const;
160 161 162
    double result;
    double *work = (double *)malloc(max(M, N)* sizeof(double));

163
    beta_const  = (CHAMELEON_Complex64_t)-1.0;
164

165
    NrowA = (side == ChamLeft) ? M : N;
166 167 168
    Anorm       = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', NrowA, NrowA, A,       LDA, work);
    Bnorm       = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M,     N,     B,       LDB, work);
    Cinitnorm   = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M,     N,     Cref,    LDC, work);
Mathieu Faverge's avatar
Mathieu Faverge committed
169
    Cchamnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M,     N,     Ccham, LDC, work);
170 171 172

    cblas_zhemm(CblasColMajor, (CBLAS_SIDE)side, (CBLAS_UPLO)uplo, M, N, CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Cref, LDC);

173
    Clapacknorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M, N, Cref, LDC, work);
174

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

177
    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'I', M, N, Cref, LDC, work);
178 179 180

    eps = LAPACKE_dlamch_work('e');

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

    result = Rnorm / ((Anorm + Bnorm + Cinitnorm) * N * eps);

    printf("============\n");
    printf("Checking the norm of the difference against reference ZHEMM \n");
Mathieu Faverge's avatar
Mathieu Faverge committed
187
    printf("-- ||Ccham - Clapack||_oo/((||A||_oo+||B||_oo+||C||_oo).N.eps) = %e \n", result );
188

Mathieu Faverge's avatar
Mathieu Faverge committed
189
    if ( isinf(Clapacknorm) || isinf(Cchamnorm) || isnan(result) || isinf(result) || (result > 10.0) ) {
190 191 192 193 194 195 196 197 198 199
        printf("-- The solution is suspicious ! \n");
        info_solution = 1;
    }
    else {
        printf("-- The solution is CORRECT ! \n");
        info_solution= 0 ;
    }
    free(work);
    return info_solution;
}