testing_zposv.c 11.8 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zposv.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 zposv 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 Bilel Hadri, 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_solution( int N, int NRHS, CHAMELEON_Complex64_t *A1, int LDA,
                CHAMELEON_Complex64_t *B1, CHAMELEON_Complex64_t *B2, int LDB, double eps );
42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

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

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

    int N     = atoi(argv[0]);
    int LDA   = atoi(argv[1]);
    int NRHS  = atoi(argv[2]);
    int LDB   = atoi(argv[3]);
    double eps;
Mathieu Faverge's avatar
Mathieu Faverge committed
62
    cham_uplo_t uplo;
63
    int info_solution, info_factorization;
Mathieu Faverge's avatar
Mathieu Faverge committed
64
    cham_trans_t trans1, trans2;
65

Mathieu Faverge's avatar
Mathieu Faverge committed
66 67 68 69
    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 *B1 = (CHAMELEON_Complex64_t *)malloc(LDB*NRHS*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *B2 = (CHAMELEON_Complex64_t *)malloc(LDB*NRHS*sizeof(CHAMELEON_Complex64_t));
70 71

    /* Check if unable to allocate memory */
72 73 74 75
    if ( (!A1) || (!A2)|| (!B1) || (!B2) )
    {
        free(A1); free(A2);
        free(B1); free(B2);
76 77 78 79
        printf("Out of Memory \n ");
        return -2;
    }

80
    eps = LAPACKE_dlamch_work( 'e' );
81

Mathieu Faverge's avatar
Mathieu Faverge committed
82 83 84
    uplo = ChamUpper;
    trans1 = uplo == ChamUpper ? ChamConjTrans : ChamNoTrans;
    trans2 = uplo == ChamUpper ? ChamNoTrans : ChamConjTrans;
85

86 87 88 89
    /*-------------------------------------------------------------
    *  TESTING ZPOSV
    */

90
    /* Initialize A1 and A2 for Symmetric Positive Matrix */
Mathieu Faverge's avatar
Mathieu Faverge committed
91 92
    CHAMELEON_zplghe( (double)N, ChamUpperLower, N, A1, LDA, 51 );
    CHAMELEON_zlacpy( ChamUpperLower, N, N, A1, LDA, A2, LDA );
93 94

    /* Initialize B1 and B2 */
Mathieu Faverge's avatar
Mathieu Faverge committed
95 96
    CHAMELEON_zplrnt( N, NRHS, B1, LDB, 371 );
    CHAMELEON_zlacpy( ChamUpperLower, N, NRHS, B1, LDB, B2, LDB );
97 98

    printf("\n");
99
    printf("------ TESTS FOR CHAMELEON ZPOSV ROUTINE -------  \n");
100 101 102 103 104 105 106
    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
107 108
    /* CHAMELEON ZPOSV */
    CHAMELEON_zposv(uplo, N, NRHS, A2, LDA, B2, LDB);
109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129

    /* Check the factorization and the solution */
    info_factorization = check_factorization( N, A1, A2, LDA, uplo, eps);
    info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

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

    /*-------------------------------------------------------------
    *  TESTING ZPOTRF + ZPOTRS
    */

    /* Initialize A1 and A2 for Symmetric Positif Matrix */
Mathieu Faverge's avatar
Mathieu Faverge committed
130 131
    CHAMELEON_zplghe( (double)N, ChamUpperLower, N, A1, LDA, 51 );
    CHAMELEON_zlacpy( ChamUpperLower, N, N, A1, LDA, A2, LDA );
132 133

    /* Initialize B1 and B2 */
Mathieu Faverge's avatar
Mathieu Faverge committed
134 135
    CHAMELEON_zplrnt( N, NRHS, B1, LDB, 371 );
    CHAMELEON_zlacpy( ChamUpperLower, N, NRHS, B1, LDB, B2, LDB );
136

Mathieu Faverge's avatar
Mathieu Faverge committed
137 138 139
    /* Cham routines */
    CHAMELEON_zpotrf(uplo, N, A2, LDA);
    CHAMELEON_zpotrs(uplo, N, NRHS, A2, LDA, B2, LDB);
140 141

    printf("\n");
142
    printf("------ TESTS FOR CHAMELEON ZPOTRF + ZPOTRS ROUTINE -------  \n");
143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169
    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 factorization and the solution */
    info_factorization = check_factorization( N, A1, A2, LDA, uplo, eps);
    info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

    if ((info_solution == 0)&(info_factorization == 0)){
        printf("***************************************************\n");
        printf(" ---- TESTING ZPOTRF + ZPOTRS ............ PASSED !\n");
        printf("***************************************************\n");
    }
    else{
        printf("****************************************************\n");
        printf(" - TESTING ZPOTRF + ZPOTRS ... FAILED !\n");
        printf("****************************************************\n");
    }

    /*-------------------------------------------------------------
    *  TESTING ZPOTRF + ZPTRSM + ZTRSM
    */

    /* Initialize A1 and A2 for Symmetric Positif Matrix */
Mathieu Faverge's avatar
Mathieu Faverge committed
170 171
    CHAMELEON_zplghe( (double)N, ChamUpperLower, N, A1, LDA, 51 );
    CHAMELEON_zlacpy( ChamUpperLower, N, N, A1, LDA, A2, LDA );
172 173

    /* Initialize B1 and B2 */
Mathieu Faverge's avatar
Mathieu Faverge committed
174 175
    CHAMELEON_zplrnt( N, NRHS, B1, LDB, 371 );
    CHAMELEON_zlacpy( ChamUpperLower, N, NRHS, B1, LDB, B2, LDB );
176

Mathieu Faverge's avatar
Mathieu Faverge committed
177 178 179
    /* CHAMELEON routines */
    CHAMELEON_zpotrf(uplo, N, A2, LDA);
    CHAMELEON_ztrsm(ChamLeft, uplo, trans1, ChamNonUnit,
180
                 N, NRHS, 1.0, A2, LDA, B2, LDB);
Mathieu Faverge's avatar
Mathieu Faverge committed
181
    CHAMELEON_ztrsm(ChamLeft, uplo, trans2, ChamNonUnit,
182 183 184
                 N, NRHS, 1.0, A2, LDA, B2, LDB);

    printf("\n");
185
    printf("------ TESTS FOR CHAMELEON ZPOTRF + ZTRSM + ZTRSM  ROUTINE -------  \n");
186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
    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 factorization and the solution */
    info_factorization = check_factorization( N, A1, A2, LDA, uplo, eps);
    info_solution = check_solution(N, NRHS, A1, LDA, B1, B2, LDB, eps);

    if ((info_solution == 0)&(info_factorization == 0)){
        printf("***************************************************\n");
        printf(" ---- TESTING ZPOTRF + ZTRSM + ZTRSM ..... PASSED !\n");
        printf("***************************************************\n");
    }
    else{
        printf("***************************************************\n");
        printf(" - TESTING ZPOTRF + ZTRSM + ZTRSM ... FAILED !\n");
        printf("***************************************************\n");
    }

208
    free(A1); free(A2); free(B1); free(B2);
209 210 211 212 213 214 215 216

    return hres;
}


/*------------------------------------------------------------------------
 *  Check the factorization of the matrix A2
 */
Mathieu Faverge's avatar
Mathieu Faverge committed
217 218 219
static int
check_factorization( int N, CHAMELEON_Complex64_t *A1, CHAMELEON_Complex64_t *A2, int LDA,
                     cham_uplo_t uplo, double eps )
220 221
{
    double Anorm, Rnorm;
Mathieu Faverge's avatar
Mathieu Faverge committed
222
    CHAMELEON_Complex64_t alpha;
223 224 225
    int info_factorization;
    int i,j;

Mathieu Faverge's avatar
Mathieu Faverge committed
226 227 228
    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));
229 230
    double *work              = (double *)malloc(N*sizeof(double));

Mathieu Faverge's avatar
Mathieu Faverge committed
231 232
    memset((void*)L1, 0, N*N*sizeof(CHAMELEON_Complex64_t));
    memset((void*)L2, 0, N*N*sizeof(CHAMELEON_Complex64_t));
233 234 235 236 237 238

    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
239
    if (uplo == ChamUpper){
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254
        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++)
           Residual[j*N+i] = L2[j*N+i] - Residual[j*N+i];

255 256
    Rnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N, Residual, N );
    Anorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N, A1, LDA );
257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279

    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 solution of the linear system
 */
Mathieu Faverge's avatar
Mathieu Faverge committed
280 281 282
static int
check_solution( int N, int NRHS, CHAMELEON_Complex64_t *A1, int LDA,
                CHAMELEON_Complex64_t *B1, CHAMELEON_Complex64_t *B2, int LDB, double eps )
283 284 285
{
    int info_solution;
    double Rnorm, Anorm, Xnorm, Bnorm, result;
Mathieu Faverge's avatar
Mathieu Faverge committed
286
    CHAMELEON_Complex64_t alpha, beta;
287 288 289 290 291
    double *work = (double *)malloc(N*sizeof(double));

    alpha = 1.0;
    beta  = -1.0;

292 293 294
    Xnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B2, LDB );
    Anorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, N,    A1, LDA );
    Bnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B1, LDB );
295 296

    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, N, NRHS, N, CBLAS_SADDR(alpha), A1, LDA, B2, LDB, CBLAS_SADDR(beta), B1, LDB);
297
    Rnorm = LAPACKE_zlange( LAPACK_COL_MAJOR, 'i', N, NRHS, B1, LDB );
298

Mathieu Faverge's avatar
Mathieu Faverge committed
299
    if (getenv("CHAMELEON_TESTING_VERBOSE"))
300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319
      printf( "||A||_oo=%f\n||X||_oo=%f\n||B||_oo=%f\n||A X - B||_oo=%e\n", Anorm, Xnorm, Bnorm, Rnorm );

    result = Rnorm / ( (Anorm*Xnorm+Bnorm)*N*eps ) ;
    printf("============\n");
    printf("Checking the Residual of the solution \n");
    printf("-- ||Ax-B||_oo/((||A||_oo||x||_oo+||B||_oo).N.eps) = %e \n", result);

    if (  isnan(Xnorm) || isinf(Xnorm) || isnan(result) || isinf(result) || (result > 60.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;
}