/**
*
* @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_factorization;
/**
* 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_zgeqrf( 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 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;
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 */
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;
}