testing_zheevd.c 9.92 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zheevd.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 zheevd testing
13
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
14
 * @version 1.0.0
15 16 17 18 19
 * @author Hatem Ltaief
 * @author Azzam Haidar
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
20
 */
21 22 23 24 25
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

Mathieu Faverge's avatar
Mathieu Faverge committed
26
#include <chameleon.h>
27 28 29
#include <coreblas/cblas.h>
#include <coreblas/lapacke.h>
#include <coreblas.h>
30 31
#include "testing_zauxiliary.h"

Mathieu Faverge's avatar
Mathieu Faverge committed
32 33
static int check_orthogonality(int, int, CHAMELEON_Complex64_t*, int, double);
static int check_reduction(cham_uplo_t, int, int, CHAMELEON_Complex64_t*, double*, int, CHAMELEON_Complex64_t*, double);
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
static int check_solution(int, double*, double*, double);

int testing_zheevd(int argc, char **argv)
{
    /* Check for number of arguments*/
    if (argc != 3) {
        USAGE("HEEVD", "MODE N LDA",
              "   - MODE : mode used to generate the matrix\n"
              "   - N    : size of the matrix A\n"
              "   - LDA  : leading dimension of the matrix A\n");
        return -1;
    }

    int    mode  = atoi(argv[0]);
    int    N     = atoi(argv[1]);
    int    LDA   = atoi(argv[2]);
    double eps   = LAPACKE_dlamch_work('e');
    double dmax  = 1.0;
    double rcond = 1.0e6;
    int INFO     = -1;

Mathieu Faverge's avatar
Mathieu Faverge committed
55 56
    cham_uplo_t uplo   = ChamUpper;
    cham_job_t vec     = ChamVec;
57 58 59 60 61
    int info_ortho     = 0;
    int info_solution  = 0;
    int info_reduction = 0;
    int LDAxN          = LDA * N;

Mathieu Faverge's avatar
Mathieu Faverge committed
62 63
    CHAMELEON_Complex64_t *A1   = NULL;
    CHAMELEON_Complex64_t *A2   = (CHAMELEON_Complex64_t *)malloc(LDAxN * sizeof(CHAMELEON_Complex64_t));
64 65
    double            *W1   = (double *)malloc(N * sizeof(double));
    double            *W2   = (double *)malloc(N * sizeof(double));
Mathieu Faverge's avatar
Mathieu Faverge committed
66 67
    CHAMELEON_Complex64_t *work = (CHAMELEON_Complex64_t *)malloc(3* N * sizeof(CHAMELEON_Complex64_t));
    CHAM_desc_t      *T;
68 69

    /* Check if unable to allocate memory */
70 71 72 73 74
    if ( (!A2) || (!W1) || (!W2) || !(work) )
    {
        free(A2);
        free(W1); free(W2);
        free(work);
75 76 77 78
        printf("Out of Memory \n ");
        return -2;
    }

Mathieu Faverge's avatar
Mathieu Faverge committed
79
    CHAMELEON_Alloc_Workspace_zheevd(N, N, &T, 1, 1);
80 81 82 83 84 85 86 87 88 89 90 91

    /*----------------------------------------------------------
    *  TESTING ZHEEVD
    */
    /* Initialize A1 */
    if (mode == 0){
        int i;
        for (i=0; i<N; i++){
            W1[i] = (double )i+1;
        }
    }
    LAPACKE_zlatms_work( LAPACK_COL_MAJOR, N, N,
Mathieu Faverge's avatar
Mathieu Faverge committed
92 93
                         chameleon_lapack_const(ChamDistSymmetric), ISEED,
                         chameleon_lapack_const(ChamHermGeev), W1, mode, rcond,
94
                         dmax, N, N,
Mathieu Faverge's avatar
Mathieu Faverge committed
95
                         chameleon_lapack_const(ChamNoPacking), A2, LDA, work );
96 97 98 99 100 101 102 103

    /*
     * Sort the eigenvalue because when computing the tridiag
     * and then the eigenvalue of the DSTQR are sorted.
     * So to avoid testing fail when having good results W1 should be sorted
     */
    LAPACKE_dlasrt_work( 'I', N, W1 );

Mathieu Faverge's avatar
Mathieu Faverge committed
104 105
    if ( vec == ChamVec ) {
        A1 = (CHAMELEON_Complex64_t *)malloc(LDAxN*sizeof(CHAMELEON_Complex64_t));
106 107 108 109 110 111

        /* Copy A2 into A1 */
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, 'A', N, N, A2, LDA, A1, LDA);
    }

    /*
Mathieu Faverge's avatar
Mathieu Faverge committed
112
     * CHAMELEON ZHEEVD
113
     */
Mathieu Faverge's avatar
Mathieu Faverge committed
114
    INFO = CHAMELEON_zheevd(vec, uplo, N, A2, LDA, W2, T);
115 116 117 118 119 120 121

    if (INFO != 0) {
        printf(" ERROR OCCURED INFO %d\n", INFO);
        goto fin;
    }

    printf("\n");
Mathieu Faverge's avatar
Mathieu Faverge committed
122
    printf("------ TESTS FOR CHAMELEON ZHEEVD ROUTINE -------  \n");
123 124 125 126 127 128 129 130
    printf("        Size of the Matrix %d by %d\n", N, 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 60.\n");

    /* Check the orthogonality, reduction and the eigen solutions */
Mathieu Faverge's avatar
Mathieu Faverge committed
131
    if (vec == ChamVec) {
132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148
        info_ortho = check_orthogonality(N, N, A2, LDA, eps);
        info_reduction = check_reduction(uplo, N, 1, A1, W2, LDA, A2, eps);
    }
    info_solution = check_solution(N, W1, W2, eps);

    if ( (info_solution == 0) & (info_ortho == 0) & (info_reduction == 0) ) {
        printf("***************************************************\n");
        printf(" ---- TESTING ZHEEVD ...................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("************************************************\n");
        printf(" - TESTING ZHEEVD ... FAILED !\n");
        printf("************************************************\n");
    }

 fin:
Mathieu Faverge's avatar
Mathieu Faverge committed
149
    CHAMELEON_Dealloc_Workspace(&T);
150 151 152 153 154 155 156 157 158 159 160 161
    free(A2);
    free(W1);
    free(W2);
    free(work);
    if (A1 != NULL) free(A1);

    return 0;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */
Mathieu Faverge's avatar
Mathieu Faverge committed
162
static int check_orthogonality(int M, int N, CHAMELEON_Complex64_t *Q, int LDQ, double eps)
163 164 165 166 167 168 169 170 171
{
    double  done  =  1.0;
    double  mdone = -1.0;
    double  normQ, result;
    int     info_ortho;
    int     minMN = min(M, N);
    double *work = (double *)malloc(minMN*sizeof(double));

    /* Build the idendity matrix */
Mathieu Faverge's avatar
Mathieu Faverge committed
172
    CHAMELEON_Complex64_t *Id = (CHAMELEON_Complex64_t *) malloc(minMN*minMN*sizeof(CHAMELEON_Complex64_t));
173 174 175 176 177 178 179 180
    LAPACKE_zlaset_work(LAPACK_COL_MAJOR, 'A', minMN, minMN, 0., 1., Id, minMN);

    /* Perform Id - Q'Q */
    if (M >= N)
        cblas_zherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, done, Q, LDQ, mdone, Id, minMN);
    else
        cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans,   M, N, done, Q, LDQ, mdone, Id, minMN);

Mathieu Faverge's avatar
Mathieu Faverge committed
181
    normQ = LAPACKE_zlanhe_work(LAPACK_COL_MAJOR, chameleon_lapack_const(ChamInfNorm), 'U', minMN, Id, minMN, work);
182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202

    result = normQ / (minMN * eps);
    printf(" ======================================================\n");
    printf(" ||Id-Q'*Q||_oo / (minMN*eps)          : %15.3E \n",  result );
    printf(" ======================================================\n");

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Orthogonality is suspicious ! \n");
        info_ortho=1;
    }
    else {
        printf("-- Orthogonality is CORRECT ! \n");
        info_ortho=0;
    }
    free(work); free(Id);
    return info_ortho;
}

/*------------------------------------------------------------
 *  Check the reduction
 */
Mathieu Faverge's avatar
Mathieu Faverge committed
203
static int check_reduction(cham_uplo_t uplo, int N, int bw, CHAMELEON_Complex64_t *A, double *D, int LDA, CHAMELEON_Complex64_t *Q, double eps )
204 205
{
    (void) bw;
Mathieu Faverge's avatar
Mathieu Faverge committed
206 207 208 209
    CHAMELEON_Complex64_t zone  =  1.0;
    CHAMELEON_Complex64_t mzone = -1.0;
    CHAMELEON_Complex64_t *TEMP     = (CHAMELEON_Complex64_t *)malloc(N*N*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *Residual = (CHAMELEON_Complex64_t *)malloc(N*N*sizeof(CHAMELEON_Complex64_t));
210 211 212 213 214 215
    double *work = (double *)malloc(N*sizeof(double));
    double Anorm, Rnorm, result;
    int info_reduction;
    int i;

    /* Compute TEMP =  Q * LAMBDA */
Mathieu Faverge's avatar
Mathieu Faverge committed
216
    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, chameleon_lapack_const(ChamUpperLower), N, N, Q, LDA, TEMP, N);
217 218 219 220 221 222 223 224 225 226

    for (i = 0; i < N; i++){
        cblas_zdscal(N, D[i], &(TEMP[i*N]), 1);
    }
    /* Compute Residual = A - Q * LAMBDA * Q^H */
    /* A is Hermetian but both upper and lower
     * are assumed valable here for checking
     * otherwise it need to be symetrized before
     * checking.
     */
Mathieu Faverge's avatar
Mathieu Faverge committed
227
    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR, chameleon_lapack_const(ChamUpperLower), N, N, A, LDA, Residual, N);
228 229
    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasConjTrans, N, N, N, CBLAS_SADDR(mzone), TEMP, N,  Q, LDA, CBLAS_SADDR(zone), Residual,     N);

Mathieu Faverge's avatar
Mathieu Faverge committed
230 231
    Rnorm = LAPACKE_zlange_work(LAPACK_COL_MAJOR, chameleon_lapack_const(ChamOneNorm), N, N, Residual, N,   work);
    Anorm = LAPACKE_zlanhe_work(LAPACK_COL_MAJOR, chameleon_lapack_const(ChamOneNorm), chameleon_lapack_const(uplo), N, A, LDA, work);
232 233

    result = Rnorm / ( Anorm * N * eps);
Mathieu Faverge's avatar
Mathieu Faverge committed
234
    if ( uplo == ChamLower ){
235 236 237 238 239 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 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
        printf(" ======================================================\n");
        printf(" ||A-Q*LAMBDA*Q'||_oo/(||A||_oo.N.eps) : %15.3E \n",  result );
        printf(" ======================================================\n");
    }else{
        printf(" ======================================================\n");
        printf(" ||A-Q'*LAMBDA*Q||_oo/(||A||_oo.N.eps) : %15.3E \n",  result );
        printf(" ======================================================\n");
    }

    if ( isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- Reduction is suspicious ! \n");
        info_reduction = 1;
    }
    else {
        printf("-- Reduction is CORRECT ! \n");
        info_reduction = 0;
    }

    free(TEMP); free(Residual);
    free(work);

    return info_reduction;
}
/*------------------------------------------------------------
 *  Check the eigenvalues
 */
static int check_solution(int N, double *E1, double *E2, double eps)
{
    int info_solution, i;
    double resid;
    double maxtmp;
    double maxel = fabs( fabs(E1[0]) - fabs(E2[0]) );
    double maxeig = max( fabs(E1[0]), fabs(E2[0]) );

    for (i = 1; i < N; i++){
        resid   = fabs(fabs(E1[i])-fabs(E2[i]));
        maxtmp  = max(fabs(E1[i]), fabs(E2[i]));

        /* Update */
        maxeig = max(maxtmp, maxeig);
        maxel  = max(resid,  maxel );
    }

    maxel = maxel / (maxeig * N * eps);
    printf(" ======================================================\n");
    printf(" | D - eigcomputed | / (|D| * N * eps) : %15.3E \n",  maxel );
    printf(" ======================================================\n");

    if ( isnan(maxel) || isinf(maxel) || (maxel > 100) ) {
        printf("-- The eigenvalues are suspicious ! \n");
        info_solution = 1;
    }
    else{
        printf("-- The eigenvalues are CORRECT ! \n");
        info_solution = 0;
    }
    return info_solution;
}