testing_zpemv.c 8.16 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zpemv.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
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
10
 ***
11
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
12
 * @brief Chameleon zpemv testing
13
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
14
 * @version 1.0.0
15 16 17 18 19 20 21
 * @author Dulceneia Becker
 * @author Mathieu Faverge
 * @author Emmanuel Agullo
 * @author Cedric Castagnede
 * @date 2011-10-06
 * @precisions normal z -> c d s
 *
22
 */
23 24 25 26 27
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

Mathieu Faverge's avatar
Mathieu Faverge committed
28
#include <chameleon.h>
29 30 31
#include <coreblas/cblas.h>
#include <coreblas/lapacke.h>
#include <coreblas.h>
32 33 34 35 36
#include "testing_zauxiliary.h"

/*--------------------------------------------------------------
 * Check the pemv
 */
Mathieu Faverge's avatar
Mathieu Faverge committed
37
static int check_solution(cham_trans_t trans, cham_store_t storev,
38
                          int M, int N, int L,
Mathieu Faverge's avatar
Mathieu Faverge committed
39 40 41 42 43
                          CHAMELEON_Complex64_t alpha, CHAMELEON_Complex64_t *A, int LDA,
                          CHAMELEON_Complex64_t *X, int INCX,
                          CHAMELEON_Complex64_t beta,  CHAMELEON_Complex64_t *Y0, int INCY0,
                          CHAMELEON_Complex64_t *Y,  int INCY,
                          CHAMELEON_Complex64_t *W, double *Rnorm)
44 45 46 47
{
    int k;
    double eps = LAPACKE_dlamch_work('e');
    double *work;
Mathieu Faverge's avatar
Mathieu Faverge committed
48
    CHAMELEON_Complex64_t mzone = -1.0;
49 50

    /* Copy x to w */
Mathieu Faverge's avatar
Mathieu Faverge committed
51
    if ( trans == ChamNoTrans ) {
52 53 54 55
        k = N;
    } else {
        k = M;
    }
56

57 58
    work = (double *)malloc(k * sizeof(double));
    cblas_zcopy(k, Y0, INCY0, W, 1);
59

60
    /* w = a A x + b w */
61
    cblas_zgemv(CblasColMajor, (CBLAS_TRANSPOSE)trans,
62
                M, N,
63 64
                CBLAS_SADDR(alpha), A,  LDA,
                X,  INCX,
65 66 67 68
                CBLAS_SADDR(beta),  W,  1);

    /* y - w */
    cblas_zaxpy(k, CBLAS_SADDR(mzone), Y, INCY, W, 1);
69

70 71
    /* Max Norm */
    *Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, 'm', 1, k, W, 1, work);
72

73 74 75 76 77
    if ( (*Rnorm / (M*N)) > eps) {
        return 1;
    } else {
        return 0;
    }
78 79

    (void)L; (void)storev;
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
}

/*--------------------------------------------------------------
 * Testing ZPEMV
 */
int testing_zpemv(int argc, char **argv)
{
    int hres = 0;
    /* Check for number of arguments*/
    if ( argc != 1) {
        USAGE("PEMV", "N",
              "   - N      : number of columns\n");
        return -1;
    }

    /* Args */
    int arg_n = atoi(argv[0]);

    /* Local variables */
Mathieu Faverge's avatar
Mathieu Faverge committed
99 100
    CHAMELEON_Complex64_t *A, *X, *Y, *A0, *Y0, *work;
    CHAMELEON_Complex64_t alpha, beta, alpha0, beta0;
101 102 103 104 105 106 107
    int n    = arg_n;
    int lda  = arg_n;

    int info_solution = 0;
    int i, j, k, t;
    int nbtests = 0;
    int nfails = 0;
Mathieu Faverge's avatar
Mathieu Faverge committed
108
    cham_store_t storev;
109 110 111 112 113 114 115 116 117
    int l = 0;
    int m = n;
    int incx = 1;
    int incy = 1;
    char *cstorev;
    double rnorm;
    double eps = LAPACKE_dlamch_work('e');

    /* Allocate Data */
Mathieu Faverge's avatar
Mathieu Faverge committed
118 119 120 121 122 123
    A    = (CHAMELEON_Complex64_t *)malloc(lda*n*sizeof(CHAMELEON_Complex64_t));
    A0   = (CHAMELEON_Complex64_t *)malloc(lda*n*sizeof(CHAMELEON_Complex64_t));
    X    = (CHAMELEON_Complex64_t *)malloc(lda*n*sizeof(CHAMELEON_Complex64_t));
    Y    = (CHAMELEON_Complex64_t *)malloc(lda*n*sizeof(CHAMELEON_Complex64_t));
    Y0   = (CHAMELEON_Complex64_t *)malloc(    n*sizeof(CHAMELEON_Complex64_t));
    work = (CHAMELEON_Complex64_t *)malloc(  2*n*sizeof(CHAMELEON_Complex64_t));
124 125 126 127 128

    LAPACKE_zlarnv_work(1, ISEED, 1, &alpha0);
    LAPACKE_zlarnv_work(1, ISEED, 1, &beta0 );

    /* Check if unable to allocate memory */
129 130 131 132
    if ( (!A) || (!A0) || (!X) || (!Y) || (!Y0) || (!work) ) {
        free(A); free(A0);
        free(X); free(Y); free(Y0);
        free(work);
133
        printf("Out of Memory \n ");
134
        return -2;
135
    }
136 137

    /* Initialize Data */
Mathieu Faverge's avatar
Mathieu Faverge committed
138 139 140
    CHAMELEON_zplrnt(n, n, A,  lda, 479 );
    CHAMELEON_zplrnt(n, n, X,  lda, 320 );
    CHAMELEON_zplrnt(n, 1, Y0, n,   573 );
141 142

    printf("\n");
143
    printf("------ TESTS FOR CHAMELEON ZPEMV ROUTINE -------  \n");
144 145 146 147 148 149 150 151
    printf("\n");
    printf(" The matrix A is randomly generated for each test.\n");
    printf(" The relative machine precision (eps) is %e \n",eps);
    printf(" Computational tests pass if scaled residual is less than eps.\n");
    printf("\n");

    nfails = 0;
    for (i=0; i<6; i++) {
152

153 154 155 156 157 158 159 160 161 162 163
        /* m and n cannot be greater than lda (arg_n) */
        switch (i) {
        case 0: l = 0;       m = arg_n;   n = m;        break;
        case 1: l = 0;       m = arg_n;   n = arg_n/2;  break; /**/
        case 2: l = arg_n;   m = l;       n = l;        break;
        case 3: l = arg_n/2; m = l;       n = arg_n;    break;
        case 4: l = arg_n/2; m = arg_n-l; n = l;        break;
        case 5: l = arg_n/3; m = arg_n-l; n = arg_n/2;  break; /**/
        }

        /* Colwise ConjTrans & Rowwise NoTrans */
164
#if defined(PRECISION_z) || defined(PRECISION_c)
165
        for (t=0; t<3; t++)
166
#else
167
        for (t=0; t<2; t++)
168
#endif
169
        {
170 171 172 173 174
            /* Swap m and n for transpose cases */
            if ( t == 1 ) {
                k = m; m = n; n = k;
            }

175
            LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A', m, n,
176 177
                                 A, lda, A0, lda);

Mathieu Faverge's avatar
Mathieu Faverge committed
178 179
            if ( trans[t] == ChamNoTrans ) {
                storev = ChamRowwise;
180 181 182 183 184 185 186 187 188 189 190
                cstorev = storevstr[0];

                /* zeroed the upper right triangle */
                int64_t i, j;
                for (j=(n-l); j<n; j++) {
                    for (i=0; i<(j-(n-l)); i++) {
                        A0[i+j*lda] = 0.0;
                    }
                }
            }
            else {
Mathieu Faverge's avatar
Mathieu Faverge committed
191
                storev = ChamColumnwise;
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216
                cstorev = storevstr[1];

                /* zeroed the lower left triangle */
                int64_t i, j;
                for (j=0; j<(l-1); j++) {
                    for (i=(m-l+1+j); i<m; i++) {
                        A0[i+j*lda] = 0.0;
                    }
                }
            }

            for (j=0; j<3; j++) {

                /* Choose alpha and beta */
                alpha = ( j==1 ) ? 0.0 : alpha0;
                beta  = ( j==2 ) ? 0.0 : beta0;

                /* incx and incy: 1 or lda */
                for (k=0; k<4; k++) {
                    switch (k) {
                    case 0:  incx = 1;    incy = 1;    break;
                    case 1:  incx = 1;    incy = lda;  break;
                    case 2:  incx = lda;  incy = 1;    break;
                    case 3:  incx = lda;  incy = lda;  break;
                    }
217

218 219
                    /* initialize Y with incy */
                    cblas_zcopy(n, Y0, 1, Y, incy);
220

221
                    /* ZPEMV */
222 223 224 225
                    CORE_zpemv( trans[t], storev, m, n, l,
                                alpha, A, lda,
                                X, incx,
                                beta,  Y, incy,
226
                                work);
227

228
                    /* Check the solution */
229 230
                    info_solution = check_solution(trans[t], storev,
                                                   m, n, l,
231
                                                   alpha, A0,  lda,
232 233 234
                                                   X,   incx,
                                                   beta,  Y0,  1,
                                                   Y,   incy,
235
                                                   work, &rnorm);
236

237 238
                    if ( info_solution != 0 ) {
                        nfails++;
239
                        printf("Failed: t=%s, s=%s, M=%3d, N=%3d, L=%3d, alpha=%e, incx=%3d, beta=%e, incy=%3d, rnorm=%e\n",
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267
                               transstr[t], cstorev, m, n, l, creal(alpha), incx, creal(beta), incy, rnorm );
                    }
                    nbtests++;
                }
            }
        }
    }

    if ( nfails )
        printf("%d / %d tests failed\n", nfails, nbtests);

    printf("***************************************************\n");
    if (nfails == 0) {
        printf(" ---- TESTING ZPEMV ...... PASSED !\n");
    }
    else {
        printf(" ---- TESTING ZPEMV ... FAILED !\n");    hres++;
    }
    printf("***************************************************\n");

    free( A0 );
    free( A );
    free( X );
    free( Y0 );
    free( Y );

    return hres;
}