Add a testing for the QR factorization of two square matrices on top of each...

Add a testing for the QR factorization of two square matrices on top of each other, and the associated generation of the Q that is required in QDWH

testing/CMakeLists.txt

@@ -89,6 +89,7 @@ set(ZSRC
     #testing_zhegst.c
     #testing_zhegv.c
     #testing_zhegvd.c
+    testing_zgeqrf_qdwh.c
     )
 
 # Add include and link directories

testing/testing_zauxiliary.c

@@ -304,6 +304,9 @@ int main (int argc, char **argv)
     /* else if ( strcmp(func, "GETMI") == 0 ) { */
     /*     info += testing_zgetmi( argc, argv ); */
     /* } */
+    else if ( strcmp(func, "GEQRF_QDWH") == 0 ) {
+        info += testing_zgeqrf_qdwh( argc, argv );
+    }
     else {
         fprintf(stderr, "Function unknown\n");
     }

testing/testing_zauxiliary.h

@@ -108,6 +108,8 @@ int testing_zhegst(int argc, char **argv);
 int testing_zgecfi(int argc, char **argv);
 int testing_zgetmi(int argc, char **argv);
 
+int testing_zgeqrf_qdwh(int argc, char **argv);
+
 #ifdef DOUBLE
 int testing_zcposv(int argc, char **argv);
 int testing_zcgesv(int argc, char **argv);

testing/testing_zgeqrf_qdwh.c

+/**
+ *
+ * @copyright (c) 2009-2014 The University of Tennessee and The University
+ *                          of Tennessee Research Foundation.
+ *                          All rights reserved.
+ * @copyright (c) 2012-2014 Inria. All rights reserved.
+ * @copyright (c) 2012-2014 Bordeaux INP, CNRS (LaBRI UMR 5800), Inria, Univ. Bordeaux. All rights reserved.
+ *
+ **/
+
+/**
+ *
+ * @file testing_zgels.c
+ *
+ *  MORSE testing routines
+ *  MORSE is a software package provided by Univ. of Tennessee,
+ *  Univ. of California Berkeley and Univ. of Colorado Denver
+ *
+ * @version 2.5.0
+ * @comment This file has been automatically generated
+ *          from Plasma 2.5.0 for MORSE 1.0.0
+ * @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
+ *
+ **/
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <math.h>
+
+#include <morse.h>
+#include <coreblas/include/cblas.h>
+#include <coreblas/include/lapacke.h>
+#include <coreblas/include/coreblas.h>
+#include <coreblas/include/coreblas_z.h>
+#include "testing_zauxiliary.h"
+
+#undef REAL
+#define COMPLEX
+
+static int check_orthogonality(int, int, const MORSE_Complex64_t*, int, double);
+static int check_factorization(int, int, const MORSE_Complex64_t*, int, const MORSE_Complex64_t*, int, MORSE_Complex64_t*, int, double);
+
+int testing_zgeqrf_qdwh(int argc, char **argv)
+{
+    int hres = 0;
+
+    if ( argc != 4 ) {
+        USAGE("GEQRF_QDWH", "optid M NB LDA",
+              "   - optid: Take into account the fact that A2 is Id or not\n"
+              "   - M    : number of rows of the matrix A1 and A2\n"
+              "   - NB   : tile size\n"
+              "   - IB   : inner tile size\n");
+        return -1;
+    }
+
+    int optid = atoi(argv[0]) ? 1: 0;
+    int M  = atoi(argv[1]);
+    int NB = atoi(argv[2]);
+    int IB = atoi(argv[3]);
+    int MxM = M * M;
+    int LDA = 2*M;
+    double eps;
+    int info_ortho, info_solution, info_factorization;
+    int i, j, rc;
+
+    /**
+     * Compute A = QR with
+     *
+     * A = [ A1 ]  and Q = [ Q1 ]
+     *     [ A2 ]        = [ Q2 ]
+     *
+     * and where A1 is the same size as A2
+     *
+     */
+    MORSE_Complex64_t *A1 = (MORSE_Complex64_t *)malloc(M*M*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *A2 = (MORSE_Complex64_t *)malloc(M*M*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *Q1 = (MORSE_Complex64_t *)malloc(M*M*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *Q2 = (MORSE_Complex64_t *)malloc(M*M*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *A  = (MORSE_Complex64_t *)malloc(2*M*M*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *Q;
+    MORSE_desc_t *T1, *T2;
+
+    /* Check if unable to allocate memory */
+    if ( (!A) || (!A1) || (!A2) || (!Q1) || (!Q2) ){
+        printf("Out of Memory \n ");
+        return -2;
+    }
+
+    MORSE_Disable(MORSE_AUTOTUNING);
+    MORSE_Set(MORSE_TILE_SIZE, NB);
+    MORSE_Set(MORSE_INNER_BLOCK_SIZE, IB);
+
+    MORSE_Alloc_Workspace_zgels(M, M, &T1, 1, 1);
+    MORSE_Alloc_Workspace_zgels(M, M, &T2, 1, 1);
+
+    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 */
+    rc = MORSE_zgeqrf( M, M, A1, M, T1 );
+    assert( rc == 0 );
+    rc = MORSE_ztpqrt( M, M, optid ? M : 0,
+                       A1, M,
+                       A2, M, T2 );
+    assert( rc == 0 );
+
+    /* Generate the Q */
+    rc = MORSE_ztpgqrt( M, M, M, (optid) ? M : 0,
+                        A1, M, T1, A2, M, T2, Q1, M, Q2, M );
+    assert( rc == 0 );
+
+    /* Copy Q in a single matrix */
+    Q = (MORSE_Complex64_t *)malloc(2*M*M*sizeof(MORSE_Complex64_t));
+    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);
+    MORSE_Dealloc_Workspace( &T1 );
+    MORSE_Dealloc_Workspace( &T2 );
+
+    return hres;
+}
+
+/*-------------------------------------------------------------------
+ * Check the orthogonality of Q
+ */
+
+static int
+check_orthogonality( int M, int N,
+                     const MORSE_Complex64_t *Q, int LDQ,
+                     double eps )
+{
+    MORSE_Complex64_t *Id;
+    double alpha, beta;
+    double normQ;
+    int info_ortho;
+    int i;
+    int minMN = min(M, N);
+
+    double *work = (double *)malloc(minMN*sizeof(double));
+
+    alpha = 1.0;
+    beta  = -1.0;
+
+    /* Build the idendity matrix */
+    Id = (MORSE_Complex64_t *) malloc(minMN*minMN*sizeof(MORSE_Complex64_t));
+    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,
+                    const MORSE_Complex64_t *A, int LDA,
+                    const MORSE_Complex64_t *R, int LDR,
+                          MORSE_Complex64_t *Q, int LDQ,
+                    double eps )
+{
+    double Anorm, Rnorm;
+    MORSE_Complex64_t alpha, beta;
+    int info_factorization;
+    int i,j;
+    double *work = (double *)malloc(max(M,N)*sizeof(double));
+
+    alpha = 1.0;
+    beta  = 0.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 */
+    CORE_zgeadd( MorseNoTrans, M, N, -1., A, LDA, 1., Q, LDQ );
+
+    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;
+}