testing_zpotri.c 8.25 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zpotri.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.
7 8
 * @copyright 2012-2014 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria,
 *                      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 zpotri testing
13
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
14
 * @version 1.0.0
15
 * @comment This file has been automatically generated
Mathieu Faverge's avatar
Mathieu Faverge committed
16
 *          from Plasma 2.5.0 for CHAMELEON 1.0.0
17 18 19 20 21 22 23
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Emmanuel Agullo
 * @author Cedric Castagnede
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
24
 */
25 26 27 28 29
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

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

Mathieu Faverge's avatar
Mathieu Faverge committed
36 37 38 39 40 41
static int
check_factorization( int N, CHAMELEON_Complex64_t *A1, CHAMELEON_Complex64_t *A2, int LDA,
                     cham_uplo_t uplo, double eps );
static int
check_inverse( int N, CHAMELEON_Complex64_t *A1, CHAMELEON_Complex64_t *A2, int LDA,
               cham_uplo_t uplo, double eps );
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

int testing_zpotri(int argc, char **argv)
{
    int hres = 0;

    /* Check for number of arguments*/
    if (argc != 2){
        USAGE("POTRI", "N LDA",
              "   - N    : the size of the matrix\n"
              "   - LDA  : leading dimension of the matrix A\n");
        return -1;
    }

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    double eps;
Mathieu Faverge's avatar
Mathieu Faverge committed
58
    cham_uplo_t uplo;
59 60
    int info_inverse, info_factorization;

Mathieu Faverge's avatar
Mathieu Faverge committed
61 62 63
    CHAMELEON_Complex64_t *A1   = (CHAMELEON_Complex64_t *)malloc(LDA*N*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *A2   = (CHAMELEON_Complex64_t *)malloc(LDA*N*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *WORK = (CHAMELEON_Complex64_t *)malloc(2*LDA*sizeof(CHAMELEON_Complex64_t));
64 65

    /* Check if unable to allocate memory */
Mathieu Faverge's avatar
Mathieu Faverge committed
66
    if ( (!A1) || (!A2) || (!WORK) )
67 68
    {
        free(A1); free(A2);
Mathieu Faverge's avatar
Mathieu Faverge committed
69
        free(WORK);
70 71 72 73
        printf("Out of Memory \n ");
        return -2;
    }

74
    eps = LAPACKE_dlamch_work( 'e' );
75

Mathieu Faverge's avatar
Mathieu Faverge committed
76
    uplo = ChamUpper;
77

78
    /*-------------------------------------------------------------
79 80
     *  TESTING ZPOTRI
     */
81 82

    /* Initialize A1 and A2 for Symmetric Positif Matrix */
Mathieu Faverge's avatar
Mathieu Faverge committed
83 84
    CHAMELEON_zplghe( (double)N, ChamUpperLower, N, A1, LDA, 51 );
    CHAMELEON_zlacpy( ChamUpperLower, N, N, A1, LDA, A2, LDA );
85 86

    printf("\n");
87
    printf("------ TESTS FOR CHAMELEON ZPOTRI ROUTINE -------  \n");
88 89 90 91 92 93 94
    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");

Mathieu Faverge's avatar
Mathieu Faverge committed
95 96
    /* CHAMELEON ZPOTRF */
    CHAMELEON_zpotrf(uplo, N, A2, LDA);
97 98 99 100

    /* Check the factorization */
    info_factorization = check_factorization( N, A1, A2, LDA, uplo, eps);

Mathieu Faverge's avatar
Mathieu Faverge committed
101 102
    /* CHAMELEON ZPOTRI */
    CHAMELEON_zpotri(uplo, N, A2, LDA);
103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126

    /* Check the inverse */
    info_inverse = check_inverse(N, A1, A2, LDA, uplo, eps);

    if ( (info_inverse == 0) && (info_factorization == 0) ) {
        printf("***************************************************\n");
        printf(" ---- TESTING ZPOTRI ..................... PASSED !\n");
        printf("***************************************************\n");
    }
    else {
        printf("***************************************************\n");
        printf(" - TESTING ZPOTRI ... FAILED !\n");    hres++;
        printf("***************************************************\n");
    }

    free(A1); free(A2); free(WORK);

    return hres;
}


/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
Mathieu Faverge's avatar
Mathieu Faverge committed
127
static int check_factorization(int N, CHAMELEON_Complex64_t *A1, CHAMELEON_Complex64_t *A2, int LDA, cham_uplo_t uplo, double eps)
128 129
{
    double Anorm, Rnorm;
Mathieu Faverge's avatar
Mathieu Faverge committed
130
    CHAMELEON_Complex64_t alpha;
131 132 133
    int info_factorization;
    int i,j;

Mathieu Faverge's avatar
Mathieu Faverge committed
134 135 136
    CHAMELEON_Complex64_t *Residual = (CHAMELEON_Complex64_t *)malloc(N*N*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *L1       = (CHAMELEON_Complex64_t *)malloc(N*N*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *L2       = (CHAMELEON_Complex64_t *)malloc(N*N*sizeof(CHAMELEON_Complex64_t));
137 138
    double *work              = (double *)malloc(N*sizeof(double));

Mathieu Faverge's avatar
Mathieu Faverge committed
139 140
    memset((void*)L1, 0, N*N*sizeof(CHAMELEON_Complex64_t));
    memset((void*)L2, 0, N*N*sizeof(CHAMELEON_Complex64_t));
141 142 143 144 145 146

    alpha= 1.0;

    LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,' ', N, N, A1, LDA, Residual, N);

    /* Dealing with L'L or U'U  */
Mathieu Faverge's avatar
Mathieu Faverge committed
147
    if (uplo == ChamUpper){
148 149 150 151 152 153 154 155 156 157 158 159 160
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasLeft, CblasUpper, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }
    else{
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L1, N);
        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', N, N, A2, LDA, L2, N);
        cblas_ztrmm(CblasColMajor, CblasRight, CblasLower, CblasConjTrans, CblasNonUnit, N, N, CBLAS_SADDR(alpha), L1, N, L2, N);
    }

    /* Compute the Residual || A -L'L|| */
    for (i = 0; i < N; i++)
        for (j = 0; j < N; j++)
161
            Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];
162

163 164
    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 );
165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188

    printf("============\n");
    printf("Checking the Cholesky Factorization \n");
    printf("-- ||L'L-A||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));

    if ( isnan(Rnorm/(Anorm*N*eps)) || isinf(Rnorm/(Anorm*N*eps)) || (Rnorm/(Anorm*N*eps) > 60.0) ){
        printf("-- Factorization is suspicious ! \n");
        info_factorization = 1;
    }
    else{
        printf("-- Factorization is CORRECT ! \n");
        info_factorization = 0;
    }

    free(Residual); free(L1); free(L2); free(work);

    return info_factorization;
}


/*------------------------------------------------------------------------
 *  Check the accuracy of the computed inverse
 */

Mathieu Faverge's avatar
Mathieu Faverge committed
189
static int check_inverse(int N, CHAMELEON_Complex64_t *A1, CHAMELEON_Complex64_t *A2, int LDA, cham_uplo_t uplo, double eps )
190 191 192 193
{
    int info_inverse;
    int i, j;
    double Rnorm, Anorm, Ainvnorm, result;
Mathieu Faverge's avatar
Mathieu Faverge committed
194 195
    CHAMELEON_Complex64_t alpha, beta, zone;
    CHAMELEON_Complex64_t *work = (CHAMELEON_Complex64_t *)malloc(N*N*sizeof(CHAMELEON_Complex64_t));
196 197 198 199 200 201

    alpha = -1.0;
    beta  = 0.0;
    zone = 1.0;

    /* Rebuild the other part of the inverse matrix */
Mathieu Faverge's avatar
Mathieu Faverge committed
202
    if(uplo == ChamUpper){
203 204 205 206 207 208
        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);
209 210 211

    }
    else {
212 213 214 215 216 217
        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);
218
    }
219

220 221 222 223
    /* Add the identity matrix to work */
    for(i=0; i<N; i++)
        *(work+i+i*N) = *(work+i+i*N) + zone;

224 225 226
    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 );
227

Mathieu Faverge's avatar
Mathieu Faverge committed
228
    if (getenv("CHAMELEON_TESTING_VERBOSE")) {
229 230
        printf( "||A||_1=%f\n||Ainv||_1=%f\n||Id - A*Ainv||_1=%e\n", Anorm, Ainvnorm, Rnorm );
    }
231 232 233 234 235 236 237 238 239

    result = Rnorm / ( (Anorm*Ainvnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the inverse \n");
    printf("-- ||Id - A*Ainv||_1/((||A||_1||Ainv||_1).N.eps) = %e \n", result);

    if (  isnan(Ainvnorm) || isinf(Ainvnorm) || isnan(result) || isinf(result) || (result > 60.0) ) {
        printf("-- The inverse is suspicious ! \n");
        info_inverse = 1;
240
    }
241 242 243 244 245 246 247 248 249
    else{
        printf("-- The inverse is CORRECT ! \n");
        info_inverse = 0;
    }

    free(work);

    return info_inverse;
}