testing_dgesvd.c File Reference

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <plasma.h>
#include <cblas.h>
#include <lapacke.h>
#include <core_blas.h>
#include "testing_dmain.h"

Include dependency graph for testing_dgesvd.c:

Go to the source code of this file.

Macros
#define	REAL

Functions
int	testing_dgesvd (int argc, char **argv)

---

Detailed Description

PLASMA testing routines PLASMA is a software package provided by Univ. of Tennessee, Univ. of California Berkeley and Univ. of Colorado Denver

Version:

2.4.5

Author:

Hatem Ltaief

Date:

2010-11-15 d Tue Nov 22 14:35:48 2011

Definition in file testing_dgesvd.c.

---

Macro Definition Documentation

#define REAL

Definition at line 27 of file testing_dgesvd.c.

---

Function Documentation

int testing_dgesvd	(	int	argc,
		char **	argv
	)

Definition at line 33 of file testing_dgesvd.c.

References check_orthogonality(), check_solution(), ISEED, lapack_const, max, min, PLASMA_Alloc_Workspace_dgesvd(), PLASMA_Dealloc_Handle_Tile(), PLASMA_dgesvd(), PLASMA_Enable(), PLASMA_ERRORS, PLASMA_HOUSEHOLDER_MODE, PLASMA_HOUSEHOLDER_SIZE, PLASMA_Set(), PLASMA_TREE_HOUSEHOLDER, PLASMA_WARNINGS, PlasmaDistUniform, PlasmaLeft, PlasmaNonsymPosv, PlasmaNoPacking, PlasmaNoVec, PlasmaRight, PlasmaVec, T, and USAGE.

{
    int tree = 0;

    if ( argc < 1 ){
        goto usage;
    } else {
        tree = atoi(argv[0]);
    }

    /* Check for number of arguments*/
    if ( ((tree == 0) && (argc != 4)) ||
         ((tree != 0) && (argc != 5)) ){
      usage:
        USAGE("GESVD", "MODE M N LDA [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"
              "   - RH   : Size of each subdomains\n");
        return -1;
    }

    int M   = atoi(argv[1]);
    int N   = atoi(argv[2]);
    int LDA = atoi(argv[3]);
    int rh;
    if ( tree ) {
        rh = atoi(argv[4]);
        PLASMA_Set(PLASMA_HOUSEHOLDER_MODE, PLASMA_TREE_HOUSEHOLDER);
        PLASMA_Set(PLASMA_HOUSEHOLDER_SIZE, rh);
    }

    if (LDA < M){
        printf("LDA should be >= M !\n");
        return -1;
    }

    double eps  = LAPACKE_dlamch_work('e');
    double dmax = 1.0;
    PLASMA_enum vecu  = PlasmaNoVec;
    PLASMA_enum vecvt = PlasmaNoVec;
    int info_orthou    = 0;
    int info_orthovt   = 0;
    int info_solution  = 0;
    int info_reduction = 0;
    int minMN = min(M, N);
    int mode  = 4;
    double rcond = (double) minMN;

    double *A1   = (double *)malloc(LDA*N*sizeof(double));
    double *S1               = (double *)            malloc(minMN*sizeof(double));
    double *S2               = (double *)            malloc(minMN*sizeof(double));
    double *work = (double *)malloc(3*max(M, N)* sizeof(double));
    double *A2 = NULL;
    double *U  = NULL;
    double *VT = NULL;
    PLASMA_desc *T;

    /* Check if unable to allocate memory */
    if ( (!A1) || (!S1) || (!S2) || (!work) ) {
        printf("Out of Memory \n ");
        return -2;
    }

    /*
    PLASMA_Disable(PLASMA_AUTOTUNING);
    PLASMA_Set(PLASMA_TILE_SIZE, 120);
    PLASMA_Set(PLASMA_INNER_BLOCK_SIZE, 40);
    */

    PLASMA_Enable(PLASMA_WARNINGS);
    PLASMA_Enable(PLASMA_ERRORS);
    PLASMA_Alloc_Workspace_dgesvd(M, N, &T);

    /*----------------------------------------------------------
    *  TESTING DGESVD
    */
    /* Initialize A1 */
    LAPACKE_dlatms_work( LAPACK_COL_MAJOR, M, N,
                         lapack_const(PlasmaDistUniform), ISEED,
                         lapack_const(PlasmaNonsymPosv), S1, mode, rcond,
                         dmax, M, N,
                         lapack_const(PlasmaNoPacking), A1, LDA, work );
    free(work);

    /* Copy A1 for check */
    if ( (vecu == PlasmaVec) && (vecvt == PlasmaVec) ) {
        A2 = (double *)malloc(LDA*N*sizeof(double));
        LAPACKE_dlacpy_work(LAPACK_COL_MAJOR, 'A', M, N, A1, LDA, A2, LDA);
    }
    if ( vecu == PlasmaVec ) {
        U = (double *)malloc(M*M*sizeof(double));
        LAPACKE_dlaset_work(LAPACK_COL_MAJOR, 'A', M, M, 0., 1., U, M);
    }
    if ( vecvt == PlasmaVec ) {
        VT = (double *)malloc(N*N*sizeof(double));
        LAPACKE_dlaset_work(LAPACK_COL_MAJOR, 'A', N, N, 0., 1., VT, N);
    }
 
    /* PLASMA DGESVD */
    PLASMA_dgesvd(vecu, vecvt, M, N, A1, LDA, S2, U, M, VT, N, T);

    if (getenv("PLASMA_TESTING_VERBOSE"))
    {
        int i;
        printf("Eigenvalues original\n");
        for (i = 0; i < min(N,25); i++){
            printf("%f ", S1[i]);
        }
        printf("\n");

        printf("Eigenvalues computed\n");
        for (i = 0; i < min(N,25); i++){
            printf("%f ", S2[i]);
        }
        printf("\n");
    }

    printf("\n");
    printf("------ TESTS FOR PLASMA DGESVD 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, reduction and the singular values */
    if ( vecu == PlasmaVec )
        info_orthou = check_orthogonality(PlasmaLeft, M, M, U, M, eps);

    if ( vecvt == PlasmaVec )
        info_orthovt = check_orthogonality(PlasmaRight, N, N, VT, N, eps);

    /* 
     * WARNING: For now, Q is associated to Band tridiagonal reduction and 
     * not to the final tridiagonal reduction, so we can not call the check
     * Need to accumulate the orthogonal transformations
     * during the bulge chasing to be able to perform the next test!
     */
    if ( (vecu == PlasmaVec) && (vecvt == PlasmaVec) && 0 )
        info_reduction = check_reduction(M, N, A2, A1, LDA, U, M, VT, N, eps);

    info_solution = check_solution(minMN, S1, S2, eps);

    if ( (info_solution == 0) & (info_orthou == 0) & 
         (info_orthovt == 0) & (info_reduction == 0) ) {
        if (M >= N) {
           printf("***************************************************\n");
           printf(" ---- TESTING DGESVD .. M >= N ........... PASSED !\n");
           printf("***************************************************\n");
        }
        else {
           printf("***************************************************\n");
           printf(" ---- TESTING DGESVD .. M < N ............ PASSED !\n");
           printf("***************************************************\n");
        }
    }
    else {
        if (M >= N) {
           printf("************************************************\n");
           printf(" - TESTING DGESVD .. M >= N .. FAILED !\n");
           printf("************************************************\n");
        }
        else {
           printf("************************************************\n");
           printf(" - TESTING DGESVD .. M < N .. FAILED !\n");
           printf("************************************************\n");
        }
    }

    if ( A2 != NULL ) free(A2); 
    if ( U  != NULL ) free(U); 
    if ( VT != NULL ) free(VT); 
    free(A1); free(S1); free(S2);
    PLASMA_Dealloc_Handle_Tile(&T);

    return 0;
}

Here is the call graph for this function:

Here is the caller graph for this function: