testing_zgeqrf_qdwh.c 8.09 KB
Newer Older
1
/**
2 3
 *
 * @file testing_zgeqrf_qdwh.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
 *
10
 ***
11
 *
12
 * @brief Chameleon zgeqrf_qdwh testing
13
 *
Mathieu Faverge's avatar
Mathieu Faverge committed
14
 * @version 1.0.0
15
 * @comment This file has been automatically generated
16
 *          from Plasma 2.5.0 for CHAMELEON 1.0.0
17 18 19 20 21 22 23 24
 * @author Bilel Hadri
 * @author Hatem Ltaief
 * @author Mathieu Faverge
 * @author Emmanuel Agullo
 * @author Cedric Castagnede
 * @date 2010-11-15
 * @precisions normal z -> c d s
 *
25
 */
26 27 28 29 30
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

31
#include <chameleon.h>
32 33 34 35
#include <coreblas/cblas.h>
#include <coreblas/lapacke.h>
#include <coreblas.h>
#include <coreblas/coreblas_z.h>
36 37
#include "testing_zauxiliary.h"

38 39
static int check_orthogonality(int, int, const CHAMELEON_Complex64_t*, int, double);
static int check_factorization(int, int, const CHAMELEON_Complex64_t*, int, const CHAMELEON_Complex64_t*, int, CHAMELEON_Complex64_t*, int, double);
40 41 42 43 44

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

Mathieu Faverge's avatar
Mathieu Faverge committed
45 46
    if ( argc != 2 ) {
        USAGE("GEQRF_QDWH", "optid M",
47
              "   - optid: Take into account the fact that A2 is Id or not\n"
Mathieu Faverge's avatar
Mathieu Faverge committed
48
              "   - M    : number of rows of the matrix A1 and A2\n");
49 50 51 52 53 54 55 56
        return -1;
    }

    int optid = atoi(argv[0]) ? 1: 0;
    int M  = atoi(argv[1]);
    int MxM = M * M;
    int LDA = 2*M;
    double eps;
Mathieu Faverge's avatar
Mathieu Faverge committed
57
    int info_ortho, info_factorization;
58 59 60 61 62 63 64 65 66 67

    /**
     * Compute A = QR with
     *
     * A = [ A1 ]  and Q = [ Q1 ]
     *     [ A2 ]        = [ Q2 ]
     *
     * and where A1 is the same size as A2
     *
     */
68 69 70 71 72 73 74
    CHAMELEON_Complex64_t *A1 = (CHAMELEON_Complex64_t *)malloc(M*M*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *A2 = (CHAMELEON_Complex64_t *)malloc(M*M*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *Q1 = (CHAMELEON_Complex64_t *)malloc(M*M*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *Q2 = (CHAMELEON_Complex64_t *)malloc(M*M*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *A  = (CHAMELEON_Complex64_t *)malloc(2*M*M*sizeof(CHAMELEON_Complex64_t));
    CHAMELEON_Complex64_t *Q;
    CHAM_desc_t *T1, *T2;
75 76 77

    /* Check if unable to allocate memory */
    if ( (!A) || (!A1) || (!A2) || (!Q1) || (!Q2) ){
78 79
        free(A); free(A1); free(A2);
        free(Q1); free(Q2);
80 81 82 83
        printf("Out of Memory \n ");
        return -2;
    }

84 85
    CHAMELEON_Alloc_Workspace_zgels(M, M, &T1, 1, 1);
    CHAMELEON_Alloc_Workspace_zgels(M, M, &T2, 1, 1);
86 87 88 89 90 91 92 93 94 95 96

    eps = LAPACKE_dlamch('e');

    /* Initialize A1, A2, and A */
    LAPACKE_zlarnv_work(IONE, ISEED, MxM, A1);
    LAPACKE_zlaset_work( LAPACK_COL_MAJOR, 'A', M, M, 0., 1., A2, M );

    LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A', M, M, A1, M, A,     LDA );
    LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A', M, M, A2, M, A + M, LDA );

    /* Factorize A */
97 98
    CHAMELEON_zgeqrf( M, M, A1, M, T1 );
    CHAMELEON_ztpqrt( M, M, optid ? M : 0,
Mathieu Faverge's avatar
Mathieu Faverge committed
99 100
                  A1, M,
                  A2, M, T2 );
101 102

    /* Generate the Q */
103
    CHAMELEON_ztpgqrt( M, M, M, (optid) ? M : 0,
Mathieu Faverge's avatar
Mathieu Faverge committed
104
                   A1, M, T1, A2, M, T2, Q1, M, Q2, M );
105 106

    /* Copy Q in a single matrix */
107
    Q = (CHAMELEON_Complex64_t *)malloc(2*M*M*sizeof(CHAMELEON_Complex64_t));
108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137
    LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A', M, M, Q1, M, Q,     LDA );
    free(Q1);
    LAPACKE_zlacpy_work( LAPACK_COL_MAJOR, 'A', M, M, Q2, M, Q + M, LDA );
    free(Q2);

    printf("\n");
    printf("------ TESTS FOR CHAMELEON ZGELS ROUTINE -------  \n");
    printf("            Size of the Matrix %d by %d\n", M, M);
    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, factorization and the solution */
    info_ortho = check_orthogonality( 2*M, M, Q, LDA, eps );
    info_factorization = check_factorization( 2*M, M, A, LDA, A1, M, Q, LDA, eps );

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

    free(A); free(A1); free(A2); free(Q);
138 139
    CHAMELEON_Dealloc_Workspace( &T1 );
    CHAMELEON_Dealloc_Workspace( &T2 );
140 141 142 143 144 145 146 147 148 149

    return hres;
}

/*-------------------------------------------------------------------
 * Check the orthogonality of Q
 */

static int
check_orthogonality( int M, int N,
150
                     const CHAMELEON_Complex64_t *Q, int LDQ,
151 152
                     double eps )
{
153
    CHAMELEON_Complex64_t *Id;
154 155 156 157 158 159 160 161 162 163 164
    double alpha, beta;
    double normQ;
    int info_ortho;
    int minMN = min(M, N);

    double *work = (double *)malloc(minMN*sizeof(double));

    alpha = 1.0;
    beta  = -1.0;

    /* Build the idendity matrix */
165
    Id = (CHAMELEON_Complex64_t *) malloc(minMN*minMN*sizeof(CHAMELEON_Complex64_t));
166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199
    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, beta, Q, LDQ, alpha, Id, N);
    else
        cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, beta, Q, LDQ, alpha, Id, M);

    normQ = LAPACKE_zlansy_work( LAPACK_COL_MAJOR, 'I', 'U', minMN, Id, minMN, work );

    printf("============\n");
    printf("Checking the orthogonality of Q \n");
    printf("||Id-Q'*Q||_oo / (N*eps) = %e \n", normQ/(minMN*eps));

    if ( isnan(normQ / (minMN * eps)) || isinf(normQ / (minMN * eps)) || (normQ / (minMN * eps) > 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 factorization QR
 */

static int
check_factorization(int M, int N,
200 201 202
                    const CHAMELEON_Complex64_t *A, int LDA,
                    const CHAMELEON_Complex64_t *R, int LDR,
                          CHAMELEON_Complex64_t *Q, int LDQ,
203 204 205
                    double eps )
{
    double Anorm, Rnorm;
206
    CHAMELEON_Complex64_t alpha;
207 208 209 210 211 212 213 214 215 216 217 218 219 220 221
    int info_factorization;
    double *work = (double *)malloc(max(M,N)*sizeof(double));

    alpha = 1.0;

    if (M >= N) {
        /* Perform Q = Q * R */
        cblas_ztrmm( CblasColMajor, CblasRight, CblasUpper, CblasNoTrans, CblasNonUnit, M, N, CBLAS_SADDR(alpha), R, LDR, Q, LDQ);
    }
    else {
        /* Perform Q = L * Q */
        cblas_ztrmm( CblasColMajor, CblasLeft, CblasLower, CblasNoTrans, CblasNonUnit, M, N, CBLAS_SADDR(alpha), R, LDR, Q, LDQ);
    }

    /* Compute the Residual */
222
    CORE_zgeadd( ChamNoTrans, M, N, -1., A, LDA, 1., Q, LDQ );
223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250

    Rnorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'I', M, N, Q, LDQ, work );
    Anorm = LAPACKE_zlange_work( LAPACK_COL_MAJOR, 'I', M, N, A, LDA, work );

    if (M >= N) {
        printf("============\n");
        printf("Checking the QR Factorization \n");
        printf("-- ||A-QR||_oo/(||A||_oo.N.eps) = %e \n",Rnorm/(Anorm*N*eps));
    }
    else {
        printf("============\n");
        printf("Checking the LQ Factorization \n");
        printf("-- ||A-LQ||_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(work);

    return info_factorization;
}