need to add the reference to libhqr for the testfile

compute/zgeqrf_param.c

@@ -39,13 +39,16 @@
  *
  *******************************************************************************
  *
+ * @param[in] qrtree
+ *          The tree used for the factorization
+ *
  * @param[in] M
  *          The number of rows of the matrix A. M >= 0.
  *
  * @param[in] N
  *          The number of columns of the matrix A.  N >= 0.
  *
- * @param[in] A
+ * @param[in, out] A
  *          On entry, the M-by-N matrix A.
  *          On exit, the elements on and above the diagonal of the array contain the min(M,N)-by-N
  *          upper trapezoidal matrix R (R is upper triangular if M >= N); the elements below the
@@ -55,7 +58,11 @@
  * @param[in] LDA
  *          The leading dimension of the array A. LDA >= max(1,M).
  *
- * @param[in] descTS
+ * @param[out] descTS
+ *          On exit, auxiliary factorization data, required by MORSE_zgeqrs to solve the system
+ *          of equations.
+ *
+ * @param[out] descTT
  *          On exit, auxiliary factorization data, required by MORSE_zgeqrs to solve the system
  *          of equations.
  *
@@ -167,7 +174,11 @@ int MORSE_zgeqrf_param(const libhqr_tree_t *qrtree, int M, int N,
  *          diagonal represent the unitary matrix Q as a product of elementary reflectors stored
  *          by tiles.
  *
- * @param[out] T
+ * @param[out] TS
+ *          On exit, auxiliary factorization data, required by MORSE_zgeqrs to solve the system
+ *          of equations.
+ *
+ * @param[out] TT
  *          On exit, auxiliary factorization data, required by MORSE_zgeqrs to solve the system
  *          of equations.
  *

testing/CMakeLists.txt

@@ -64,6 +64,7 @@ set(ZSRC
     # LAPACK
     ##################
     testing_zgels.c
+    testing_zgeqrf_param.c
     #testing_zgesv.c
     testing_zgesv_incpiv.c
     #testing_zgetri.c
@@ -90,7 +91,6 @@ set(ZSRC
     #testing_zhegv.c
     #testing_zhegvd.c
     testing_zgeqrf_qdwh.c
-    testing_zgeqrfhqr_qdwh.c
     )
 
 # Add include and link directories

testing/testing_zgeqrf_param.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_zgeqrf_param.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
+ * @author Boucherie Raphael
+ * @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 "testing_zauxiliary.h"
+
+#undef REAL
+#define COMPLEX
+
+enum blas_order_type {
+            blas_rowmajor = 101,
+            blas_colmajor = 102 };
+
+enum blas_uplo_type  {
+            blas_upper = 121,
+            blas_lower = 122 };
+
+enum blas_cmach_type {
+            blas_base      = 151,
+            blas_t         = 152,
+            blas_rnd       = 153,
+            blas_ieee      = 154,
+            blas_emin      = 155,
+            blas_emax      = 156,
+            blas_eps       = 157,
+            blas_prec      = 158,
+            blas_underflow = 159,
+            blas_overflow  = 160,
+            blas_sfmin     = 161};
+
+enum blas_norm_type {
+            blas_one_norm       = 171,
+            blas_real_one_norm  = 172,
+            blas_two_norm       = 173,
+            blas_frobenius_norm = 174,
+            blas_inf_norm       = 175,
+            blas_real_inf_norm  = 176,
+            blas_max_norm       = 177,
+            blas_real_max_norm  = 178 };
+
+static void
+BLAS_error(char *rname, int err, int val, int x) {
+  fprintf( stderr, "%s %d %d %d\n", rname, err, val, x );
+  abort();
+}
+
+static
+void
+BLAS_zsy_norm(enum blas_order_type order, enum blas_norm_type norm,
+  enum blas_uplo_type uplo, int n, const MORSE_Complex64_t *a, int lda, double *res) {
+  int i, j; double anorm, v;
+  char rname[] = "BLAS_zsy_norm";
+
+  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );
+
+  if (norm == blas_inf_norm) {
+    anorm = 0.0;
+    if (blas_upper == uplo) {
+      for (i = 0; i < n; ++i) {
+        v = 0.0;
+        for (j = 0; j < i; ++j) {
+          v += cabs( a[j + i * lda] );
+        }
+        for (j = i; j < n; ++j) {
+          v += cabs( a[i + j * lda] );
+        }
+        if (v > anorm)
+          anorm = v;
+      }
+    } else {
+      BLAS_error( rname, -3, norm, 0 );
+      return;
+    }
+  } else {
+    BLAS_error( rname, -2, norm, 0 );
+    return;
+  }
+
+  if (res) *res = anorm;
+}
+
+static
+void
+BLAS_zge_norm(enum blas_order_type order, enum blas_norm_type norm,
+  int m, int n, const MORSE_Complex64_t *a, int lda, double *res) {
+  int i, j; float anorm, v;
+  char rname[] = "BLAS_zge_norm";
+
+  if (order != blas_colmajor) BLAS_error( rname, -1, order, 0 );
+
+  if (norm == blas_frobenius_norm) {
+    anorm = 0.0f;
+    for (j = n; j; --j) {
+      for (i = m; i; --i) {
+        v = a[0];
+        anorm += v * v;
+        a++;
+      }
+      a += lda - m;
+    }
+    anorm = sqrt( anorm );
+  } else if (norm == blas_inf_norm) {
+    anorm = 0.0f;
+    for (i = 0; i < m; ++i) {
+      v = 0.0f;
+      for (j = 0; j < n; ++j) {
+        v += cabs( a[i + j * lda] );
+      }
+      if (v > anorm)
+        anorm = v;
+    }
+  } else {
+    BLAS_error( rname, -2, norm, 0 );
+    return;
+  }
+
+  if (res) *res = anorm;
+}
+
+static
+double
+BLAS_dpow_di(double x, int n) {
+  double rv = 1.0;
+
+  if (n < 0) {
+    n = -n;
+    x = 1.0 / x;
+  }
+
+  for (; n; n >>= 1, x *= x) {
+    if (n & 1)
+      rv *= x;
+  }
+
+  return rv;
+}
+
+static
+double
+BLAS_dfpinfo(enum blas_cmach_type cmach) {
+  double eps = 1.0, r = 1.0, o = 1.0, b = 2.0;
+  int t = 53, l = 1024, m = -1021;
+  char rname[] = "BLAS_dfpinfo";
+
+  if ((sizeof eps) == sizeof(float)) {
+    t = 24;
+    l = 128;
+    m = -125;
+  } else {
+    t = 53;
+    l = 1024;
+    m = -1021;
+  }
+
+  /* for (i = 0; i < t; ++i) eps *= half; */
+  eps = BLAS_dpow_di( b, -t );
+  /* for (i = 0; i >= m; --i) r *= half; */
+  r = BLAS_dpow_di( b, m-1 );
+
+  o -= eps;
+  /* for (i = 0; i < l; ++i) o *= b; */
+  o = (o * BLAS_dpow_di( b, l-1 )) * b;
+
+  switch (cmach) {
+    case blas_eps: return eps;
+    case blas_sfmin: return r;
+    default:
+      BLAS_error( rname, -1, cmach, 0 );
+      break;
+  }
+  return 0.0;
+}
+
+static int check_orthogonality(int, int, int, MORSE_Complex64_t*, double);
+static int check_factorization(int, int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, MORSE_Complex64_t*, double);
+static int check_solution(int, int, int, MORSE_Complex64_t*, int, MORSE_Complex64_t*, MORSE_Complex64_t*, int, double);
+
+int testing_zgeqrf_param(int argc, char **argv)
+{
+    int hres = 0;
+    int mode = 0;
+
+    if ( argc < 1 ){
+        goto usage;
+    } else {
+        mode = atoi(argv[0]);
+    }
+
+    /* Check for number of arguments*/
+    if ( ((mode == 0) && (argc != 6)) ||
+         ((mode != 0) && (argc != 7)) ){
+      usage:
+        USAGE("GEQRF_PARAM", "MODE M N LDA NRHS LDB [RH]",
+              "   - MODE : 0: flat, 1: tree (RH needed)\n"
+              "   - M    : number of rows of the matrix A\n"
+              "   - N    : number of columns of the matrix A\n"
+              "   - LDA  : leading dimension of the matrix A\n"
+              "   - NRHS : number of RHS\n"
+              "   - LDB  : leading dimension of the matrix B\n"
+              "   - RH   : Size of each subdomains\n");
+        return -1;
+    }
+
+    int M     = atoi(argv[1]);
+    int N     = atoi(argv[2]);
+    int LDA   = max( atoi(argv[3]), M );
+    int NRHS  = atoi(argv[4]);
+    int LDB   = max( max( atoi(argv[5]), M ), N );
+    int rh;
+    libhqr_tree_t qrtree;
+    libhqr_tiledesc_t matrix;
+
+    int K = min(M, N);
+    double eps;
+    int info_ortho, info_solution, info_factorization;
+    int i,j;
+    int LDAxN = LDA*N;
+    int LDBxNRHS = LDB*NRHS;
+    int domino, tsrr, llvl, hlvl, qr_a, qr_p;
+
+    MORSE_Complex64_t *A1 = (MORSE_Complex64_t *)malloc(LDA*N*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *A2 = (MORSE_Complex64_t *)malloc(LDA*N*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *B1 = (MORSE_Complex64_t *)malloc(LDB*NRHS*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *B2 = (MORSE_Complex64_t *)malloc(LDB*NRHS*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *Q  = (MORSE_Complex64_t *)malloc(LDA*N*sizeof(MORSE_Complex64_t));
+    MORSE_desc_t *TS;
+    MORSE_desc_t *TT;
+
+    /* Check if unable to allocate memory */
+    if ((!A1)||(!A2)||(!B1)||(!B2)||(!Q)){
+        printf("Out of Memory \n ");
+        return -2;
+    }
+
+    if ( mode ) {
+        rh = atoi(argv[6]);
+
+        MORSE_Set(MORSE_HOUSEHOLDER_MODE, MORSE_TREE_HOUSEHOLDER);
+        MORSE_Set(MORSE_HOUSEHOLDER_SIZE, rh);
+    }
+
+    MORSE_Alloc_Workspace_zgels(M, N, &TS, 1, 1);
+    MORSE_Alloc_Workspace_zgels(M, N, &TT, 1, 1);
+    memset(TS->mat, 0, (TS->llm*TS->lln)*sizeof(MORSE_Complex64_t));
+    memset(TT->mat, 0, (TT->llm*TT->lln)*sizeof(MORSE_Complex64_t));
+    eps = BLAS_dfpinfo( blas_eps );
+
+    /*----------------------------------------------------------
+    *  TESTING ZGEQRF_PARAM
+    */
+
+    /* Initialize matrix */
+    matrix.mt = TS->mt;
+    matrix.nt = TS->nt;
+    matrix.nodes = 1;
+    matrix.p = 1;
+    /* Initialize qrtree */
+    libhqr_hqr_init( &qrtree, 0, &matrix, llvl, hlvl, qr_a, qr_p, domino, tsrr);
+
+    /* Initialize A1 and A2 */
+    LAPACKE_zlarnv_work(IONE, ISEED, LDAxN, A1);
+    for (i = 0; i < M; i++)
+        for (j = 0; j < N; j++)
+            A2[LDA*j+i] = A1[LDA*j+i] ;
+
+    /* Initialize B1 and B2 */
+    memset(B2, 0, LDB*NRHS*sizeof(MORSE_Complex64_t));
+    LAPACKE_zlarnv_work(IONE, ISEED, LDBxNRHS, B1);
+    for (i = 0; i < M; i++)
+        for (j = 0; j < NRHS; j++)
+            B2[LDB*j+i] = B1[LDB*j+i] ;
+
+    /* MORSE ZGEQRF_PARAM */
+    MORSE_zgeqrf_param(&qrtree, M, N, A2, LDA, TS, TT);
+
+    /* MORSE ZGELS */
+    if (M >= N)
+       /* Building the economy-size Q */
+       MORSE_zungqr(M, N, K, A2, LDA, TS, Q, LDA);
+    else
+       /* Building the economy-size Q */
+       MORSE_zunglq(M, N, K, A2, LDA, TS, Q, LDA);
+
+    printf("\n");
+    printf("------ TESTS FOR CHAMELEON ZGEQRF_PARAM ROUTINE -------  \n");
+    printf("            Size of the Matrix %d by %d\n", M, 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, factorization and the solution */
+    info_ortho = check_orthogonality(M, N, LDA, Q, eps);
+    info_factorization = check_factorization(M, N, A1, A2, LDA, Q, eps);
+    info_solution = check_solution(M, N, NRHS, A1, LDA, B1, B2, LDB, eps);
+
+    if ((info_solution == 0)&(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");
+    }
+
+    libhqr_matrix_finalize( &qrtree );
+    free(A1); free(A2); free(B1); free(B2); free(Q);
+    MORSE_Dealloc_Workspace( &TS );
+    MORSE_Dealloc_Workspace( &TT );
+
+    return hres;
+}
+
+/*-------------------------------------------------------------------
+ * Check the orthogonality of Q
+ */
+
+static int check_orthogonality(int M, int N, int LDQ, MORSE_Complex64_t *Q, double eps)
+{
+    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 USE DLASET?*/
+    MORSE_Complex64_t *Id = (MORSE_Complex64_t *) malloc(minMN*minMN*sizeof(MORSE_Complex64_t));
+    memset((void*)Id, 0, minMN*minMN*sizeof(MORSE_Complex64_t));
+    for (i = 0; i < minMN; i++)
+        Id[i*minMN+i] = (MORSE_Complex64_t)1.0;
+
+    /* Perform Id - Q'Q */
+    if (M >= N)
+        cblas_zherk(CblasColMajor, CblasUpper, CblasConjTrans, N, M, alpha, Q, LDQ, beta, Id, N);
+    else
+        cblas_zherk(CblasColMajor, CblasUpper, CblasNoTrans, M, N, alpha, Q, LDQ, beta, Id, M);
+
+    BLAS_zsy_norm( blas_colmajor, blas_inf_norm, blas_upper, minMN, Id, minMN, &normQ );
+
+    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, MORSE_Complex64_t *A1, MORSE_Complex64_t *A2, int LDA, MORSE_Complex64_t *Q, double eps )
+{
+    double Anorm, Rnorm;
+    MORSE_Complex64_t alpha, beta;
+    int info_factorization;
+    int i,j;
+
+    MORSE_Complex64_t *Ql       = (MORSE_Complex64_t *)malloc(M*N*sizeof(MORSE_Complex64_t));
+    MORSE_Complex64_t *Residual = (MORSE_Complex64_t *)malloc(M*N*sizeof(MORSE_Complex64_t));
+    double *work              = (double *)malloc(max(M,N)*sizeof(double));
+
+    alpha=1.0;
+    beta=0.0;
+
+    if (M >= N) {
+        /* Extract the R */
+        MORSE_Complex64_t *R = (MORSE_Complex64_t *)malloc(N*N*sizeof(MORSE_Complex64_t));
+        memset((void*)R, 0, N*N*sizeof(MORSE_Complex64_t));
+        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'u', M, N, A2, LDA, R, N);
+
+        /* Perform Ql=Q*R */
+        memset((void*)Ql, 0, M*N*sizeof(MORSE_Complex64_t));
+        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, N, CBLAS_SADDR(alpha), Q, LDA, R, N, CBLAS_SADDR(beta), Ql, M);
+        free(R);
+    }
+    else {
+        /* Extract the L */
+        MORSE_Complex64_t *L = (MORSE_Complex64_t *)malloc(M*M*sizeof(MORSE_Complex64_t));
+        memset((void*)L, 0, M*M*sizeof(MORSE_Complex64_t));
+        LAPACKE_zlacpy_work(LAPACK_COL_MAJOR,'l', M, N, A2, LDA, L, M);
+
+    /* Perform Ql=LQ */
+        memset((void*)Ql, 0, M*N*sizeof(MORSE_Complex64_t));
+        cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, N, M, CBLAS_SADDR(alpha), L, M, Q, LDA, CBLAS_SADDR(beta), Ql, M);
+        free(L);
+    }
+
+    /* Compute the Residual */
+    for (i = 0; i < M; i++)
+        for (j = 0 ; j < N; j++)
+            Residual[j*M+i] = A1[j*LDA+i]-Ql[j*M+i];
+
+    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, Residual, M, &Rnorm );
+    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N, A2, LDA, &Anorm );
+
+    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); free(Ql); free(Residual);
+
+    return info_factorization;
+}
+
+/*--------------------------------------------------------------
+ * Check the solution
+ */
+
+static int check_solution(int M, int N, int NRHS, MORSE_Complex64_t *A, int LDA, MORSE_Complex64_t *B, MORSE_Complex64_t *X, int LDB, double eps)
+{
+    int info_solution;
+    double Rnorm, Anorm, Xnorm, Bnorm;
+    MORSE_Complex64_t alpha, beta;
+    double result;
+    double *work = (double *)malloc(max(M, N)* sizeof(double));
+
+    alpha = 1.0;
+    beta  = -1.0;
+
+    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, N,    A, LDA, &Anorm );
+    BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, B, LDB, &Bnorm );
+    BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, NRHS, X, LDB, &Xnorm );
+
+    cblas_zgemm(CblasColMajor, CblasNoTrans, CblasNoTrans, M, NRHS, N, CBLAS_SADDR(alpha), A, LDA, X, LDB, CBLAS_SADDR(beta), B, LDB);
+
+    if (M >= N) {
+       MORSE_Complex64_t *Residual = (MORSE_Complex64_t *)malloc(M*NRHS*sizeof(MORSE_Complex64_t));
+       memset((void*)Residual, 0, M*NRHS*sizeof(MORSE_Complex64_t));
+       cblas_zgemm(CblasColMajor, CblasConjTrans, CblasNoTrans, N, NRHS, M, CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Residual, M);
+       BLAS_zge_norm( blas_colmajor, blas_inf_norm, M, NRHS, Residual, M, &Rnorm );
+       free(Residual);
+    }
+    else {
+       MORSE_Complex64_t *Residual = (MORSE_Complex64_t *)malloc(N*NRHS*sizeof(MORSE_Complex64_t));
+       memset((void*)Residual, 0, N*NRHS*sizeof(MORSE_Complex64_t));
+       cblas_zgemm(CblasColMajor, CblasConjTrans, CblasNoTrans, N, NRHS, M, CBLAS_SADDR(alpha), A, LDA, B, LDB, CBLAS_SADDR(beta), Residual, N);
+       BLAS_zge_norm( blas_colmajor, blas_inf_norm, N, NRHS, Residual, N, &Rnorm );
+       free(Residual);
+    }
+
+    if (getenv("MORSE_TESTING_VERBOSE"))
+       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;
+}

testing/testing_zgeqrfhqr_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_zgeqrfhqr_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;
-
-    /**
-     * 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 */
-    MORSE_zgeqrfhqr( M, M, A1, M, T1 );
-    MORSE_ztpqrt( M, M, optid ? M : 0,
-                  A1, M,
-                  A2, M, T2 );
-
-    /* Generate the Q */
-    MORSE_ztpgqrt( M, M, M, (optid) ? M : 0,
-                   A1, M, T1, A2, M, T2, Q1, M, Q2, M );
-
-    /* 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;
-}