zgesvd.c

Go to the documentation of this file.

#include <lapacke.h>
#include "common.h"

/***************************************************************************/
int PLASMA_zgesvd(PLASMA_enum jobu, PLASMA_enum jobvt, int M, int N,
                  PLASMA_Complex64_t *A, int LDA,
                  double *S,
                  PLASMA_Complex64_t *U, int LDU,
                  PLASMA_Complex64_t *VT, int LDVT,
                  PLASMA_desc *descT)
{
    int NB, IB, MT, NT;
    int status;
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    PLASMA_desc descA, descU, descVT;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zgesvd", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    
    /* Tune NB & IB depending on M & N; Set NBNB */
    status = plasma_tune(PLASMA_FUNC_ZGESVD, M, N, 0);
    if (status != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zgesvd", "plasma_tune() failed");
        return status;
    }

    /* Set MT, NT */
    NB = PLASMA_NB;
    IB = PLASMA_IB;
    MT = (M%NB==0) ? (M/NB) : (M/NB+1);
    NT = (N%NB==0) ? (N/NB) : (N/NB+1);

    /* Check input arguments */
    if (jobu != PlasmaNoVec  && jobu !=PlasmaVec) {
        plasma_error("PLASMA_zgesvd", "illegal value of jobu");
        return -1;
    }
    if (jobvt != PlasmaNoVec && jobvt != PlasmaVec) {
        plasma_error("PLASMA_zgesvd", "illegal value of jobvt");
        return -2;
    }
    if (M < 0) {
        plasma_error("PLASMA_zgesvd", "illegal value of M");
        return -3;
    }
    if (N < 0) {
        plasma_error("PLASMA_zgesvd", "illegal value of N");
        return -4;
    }
    if (LDA < max(1, M)) {
        plasma_error("PLASMA_zgesvd", "illegal value of LDA");
        return -6;
    }
    if (LDU < 1) {
        plasma_error("PLASMA_zgesvd", "illegal value of LDU");
        return -9;
    }
    if (LDVT < 1) {
        plasma_error("PLASMA_zgesvd", "illegal value of LDVT");
        return -11;
    }
    if ( (plasma_desc_check(descT) != PLASMA_SUCCESS) || 
         ( descT->m != MT*IB ) || (descT->n != NT*NB) ) {
        plasma_error("PLASMA_zgesvd", "invalid T descriptor");
        return -12;
    }
    /* Quick return */
    if (min(M, N) == 0) {
        return PLASMA_SUCCESS;
    }

    if (jobu == PlasmaVec) {
        plasma_error("PLASMA_zgesvd", "computing the singular vectors is not supported in this version");
        return -1;
    }
    if (jobvt == PlasmaVec) {
        plasma_error("PLASMA_zgesvd", "computing the singular vectors is not supported in this version");
        return -2;
    }

    plasma_sequence_create(plasma, &sequence);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_zooplap2tile( descA,   A, NB, NB,  LDA, N, 0, 0, M, N, plasma_desc_mat_free(&(descA)) );
        if (jobu == PlasmaVec){
            plasma_zooplap2tile( descU,   U, NB, NB,  LDU, M, 0, 0, M, M, plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descU)));
        }
        if (jobvt == PlasmaVec){
            plasma_zooplap2tile( descVT, VT, NB, NB, LDVT, N, 0, 0, N, N, plasma_desc_mat_free(&(descA)); plasma_desc_mat_free(&(descU)); plasma_desc_mat_free(&(descVT)));
        }
    } else {
        plasma_ziplap2tile( descA,   A, NB, NB,  LDA, N, 0, 0, M, N);
        if (jobu == PlasmaVec){
            plasma_ziplap2tile( descU,   U, NB, NB,  LDU, M, 0, 0, M, M);
        }
        if (jobvt == PlasmaVec){
            plasma_ziplap2tile( descVT, VT, NB, NB, LDVT, N, 0, 0, N, N);
        }
    }

    /* Call the tile interface */
    PLASMA_zgesvd_Tile_Async(jobu, jobvt, &descA, S, &descU, &descVT, descT, sequence, &request);

    if ( PLASMA_TRANSLATION == PLASMA_OUTOFPLACE ) {
        plasma_zooptile2lap( descA,   A, NB, NB,  LDA, N );
        if (jobu == PlasmaVec){
            plasma_zooptile2lap( descU,   U, NB, NB,  LDU, M );
        }
        if (jobvt == PlasmaVec){
            plasma_zooptile2lap( descVT, VT, NB, NB, LDVT, N );
        }
        plasma_dynamic_sync();
        plasma_desc_mat_free(&descA);
        if (jobu == PlasmaVec){
            plasma_desc_mat_free(&descU);
        }
        if (jobvt == PlasmaVec){
            plasma_desc_mat_free(&descVT);
        }
    } else {
        plasma_ziptile2lap( descA,   A, NB, NB,  LDA, N );
        if (jobu == PlasmaVec){
            plasma_ziptile2lap( descU,   U, NB, NB,  LDU, M );
        }
        if (jobvt == PlasmaVec){
            plasma_ziptile2lap( descVT, VT, NB, NB, LDVT, N );
        }
        plasma_dynamic_sync();
    }

    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************/
int PLASMA_zgesvd_Tile(PLASMA_enum jobu, PLASMA_enum jobvt, PLASMA_desc *A,
                      double *S, PLASMA_desc *U, PLASMA_desc *VT, PLASMA_desc *T)
{
    plasma_context_t *plasma;
    PLASMA_sequence *sequence = NULL;
    PLASMA_request request = PLASMA_REQUEST_INITIALIZER;
    int status;

    plasma = plasma_context_self();
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zgesvd_Tile", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    plasma_sequence_create(plasma, &sequence);
    PLASMA_zgesvd_Tile_Async(jobu, jobvt, A, S, U, VT, T, sequence, &request);
    plasma_dynamic_sync();
    status = sequence->status;
    plasma_sequence_destroy(plasma, sequence);
    return status;
}

/***************************************************************************/
int PLASMA_zgesvd_Tile_Async(PLASMA_enum jobu, PLASMA_enum jobvt, PLASMA_desc *A,
                             double *S, PLASMA_desc *U, PLASMA_desc *VT, PLASMA_desc *T,
                             PLASMA_sequence *sequence, PLASMA_request *request)
{
    PLASMA_desc descA  = *A;
    PLASMA_desc descT  = *T;
    double *E;
    int minMN = min(descA.m, descA.n);
    int NCVT = 0;
    int NRU = 0;
    int NCC = 0;

    plasma_context_t *plasma;

    plasma = plasma_context_self();

    if (jobu != PlasmaNoVec  && jobu !=PlasmaVec) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "illegal value of jobu");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (jobvt != PlasmaNoVec && jobvt != PlasmaVec) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "illegal value of jobvt");
        return PLASMA_ERR_NOT_SUPPORTED;
    }
    if (plasma == NULL) {
        plasma_fatal_error("PLASMA_zgesvd_Tile_Async", "PLASMA not initialized");
        return PLASMA_ERR_NOT_INITIALIZED;
    }
    if (sequence == NULL) {
        plasma_fatal_error("PLASMA_zgesvd_Tile_Async", "NULL sequence");
        return PLASMA_ERR_UNALLOCATED;
    }
    if (request == NULL) {
        plasma_fatal_error("PLASMA_zgesvd_Tile_Async", "NULL request");
        return PLASMA_ERR_UNALLOCATED;
    }
    /* Check sequence status */
    if (sequence->status == PLASMA_SUCCESS)
        request->status = PLASMA_SUCCESS;
    else
        return plasma_request_fail(sequence, request, PLASMA_ERR_SEQUENCE_FLUSHED);

    /* Check descriptors for correctness */
    if (plasma_desc_check(&descA) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "invalid first descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if ((jobu != PlasmaNoVec) && (plasma_desc_check(U) != PLASMA_SUCCESS)) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "invalid second descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if ((jobvt != PlasmaNoVec) && (plasma_desc_check(VT) != PLASMA_SUCCESS) ) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "invalid third descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (plasma_desc_check(&descT) != PLASMA_SUCCESS) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "invalid fourth descriptor");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    /* Check input arguments */
    if (descA.nb != descA.mb) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if (( (jobu != PlasmaNoVec) && (U->nb != U->mb) )  || ( (jobvt != PlasmaNoVec) && (VT->nb != VT->mb) )) {
        plasma_error("PLASMA_zgesvd_Tile_Async", "only square tiles supported");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    if ((jobu == PlasmaVec) || (jobvt == PlasmaVec) ){
        plasma_error("PLASMA_zgesvd_Tile_Async", "computing the singular vectors is not supported in this version");
        return plasma_request_fail(sequence, request, PLASMA_ERR_ILLEGAL_VALUE);
    }
    
    E = (double *)plasma_shared_alloc(plasma, minMN-1, PlasmaRealDouble);
 
    /* 
     * Reduction to bidiagonal form with a two-stage approach.
     */

    /* 
     * 1: Reduction to BAND bidiagonal form
     * May be further optimized using the algo described in Trefethen
     */
    /* if (plasma->householder == PLASMA_FLAT_HOUSEHOLDER) { */
        plasma_dynamic_call_4(plasma_pzgerbb,
            PLASMA_desc, descA,
            PLASMA_desc, descT,
            PLASMA_sequence*, sequence,
            PLASMA_request*, request);
    /* } */
    /* else { */
    /*     plasma_dynamic_call_4(plasma_pzgerbbrh, */
    /*         PLASMA_desc, descA, */
    /*         PLASMA_desc, descT, */
    /*         PLASMA_sequence*, sequence, */
    /*         PLASMA_request*, request); */
    /* } */

    /* Build the U of the first stage */
    /* if (jobu == PlasmaVec){ */
    /*    /\* Initialize U to Identity *\/ */
    /*    plasma_dynamic_call_6(plasma_pzlaset, */
    /*        PLASMA_enum, PlasmaUpperLower, */
    /*        PLASMA_Complex64_t, 0.0, */
    /*        PLASMA_Complex64_t, 1.0, */
    /*        PLASMA_desc, descU, */
    /*        PLASMA_sequence*, sequence, */
    /*        PLASMA_request*, request); */
    /*    /\* Accumulate the transformations from the first stage *\/ */
    /*    if (plasma->householder == PLASMA_FLAT_HOUSEHOLDER) { */
    /*        plasma_dynamic_call_6(plasma_pzungbr, */
    /*            PLASMA_enum, PlasmaLeft, */
    /*            PLASMA_desc, descA, */
    /*            PLASMA_desc, descU, */
    /*            PLASMA_desc, descT, */
    /*            PLASMA_sequence*, sequence, */
    /*            PLASMA_request*, request); */
    /*    } */
    /*    else { */
    /*        plasma_dynamic_call_6(plasma_pzungbrrh, */
    /*            PLASMA_enum, PlasmaLeft, */
    /*            PLASMA_desc, descA, */
    /*            PLASMA_desc, descU, */
    /*            PLASMA_desc, descT, */
    /*            PLASMA_sequence*, sequence, */
    /*            PLASMA_request*, request); */
    /*    } */
    /* } */
    
    /* Build the VT of the first stage */
    /* if (jobvt == PlasmaVec){ */
    /*     /\* Initialize VT to Identity *\/ */
    /*    plasma_dynamic_call_6(plasma_pzlaset, */
    /*        PLASMA_enum, PlasmaUpperLower, */
    /*        PLASMA_Complex64_t, 0.0, */
    /*        PLASMA_Complex64_t, 1.0, */
    /*        PLASMA_desc, descVT, */
    /*        PLASMA_sequence*, sequence, */
    /*        PLASMA_request*, request); */

    /*    /\* Accumulate the transformations from the first stage *\/ */
    /*    if (plasma->householder == PLASMA_FLAT_HOUSEHOLDER) { */
    /*        plasma_dynamic_call_6(plasma_pzungbr, */
    /*            PLASMA_enum, PlasmaRight, */
    /*            PLASMA_desc, descA, */
    /*            PLASMA_desc, descVT, */
    /*            PLASMA_desc, descT, */
    /*            PLASMA_sequence*, sequence, */
    /*            PLASMA_request*, request); */
    /*    } */
    /*    else { */
    /*        plasma_dynamic_call_6(plasma_pzungbrrh, */
    /*            PLASMA_enum, PlasmaRight, */
    /*            PLASMA_desc, descA, */
    /*            PLASMA_desc, descVT, */
    /*            PLASMA_desc, descT, */
    /*            PLASMA_sequence*, sequence, */
    /*            PLASMA_request*, request); */
    /*    } */
    /* } */

    /*
     * Set the V's to zero before the 2nd stage i.e., bulge chasing
     */
    plasma_dynamic_call_5(plasma_pzlaset2,
        PLASMA_enum, PlasmaLower,
        PLASMA_Complex64_t, 0.0,
        PLASMA_desc, descA.m >= descA.n ? descA : plasma_desc_submatrix(descA, descA.mb, 0, descA.m-descA.mb, descA.n),
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    plasma_dynamic_call_5(plasma_pzlaset2,
        PLASMA_enum, PlasmaUpper,
        PLASMA_Complex64_t, 0.0,
        PLASMA_desc, descA.m >= descA.n ? plasma_desc_submatrix(descA, 0, descA.nb, descA.m, descA.n-descA.nb) : descA,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /* 
     * 2: Reduction from BAND bidiagonal to the final condensed form
     */
    plasma_dynamic_call_7(plasma_pzgbrdb,
        PLASMA_enum, descA.m >= descA.n ? PlasmaUpper : PlasmaLower,
        PLASMA_desc, descA,
        double*, S,
        double*, E,
        PLASMA_desc, descT,
        PLASMA_sequence*, sequence,
        PLASMA_request*, request);

    /*
     * Compute the singular values ONLY for now
     */
    plasma_dynamic_sync();
    if (descA.m >= descA.n)
       LAPACKE_zbdsqr(
              LAPACK_COL_MAJOR, lapack_const(PlasmaUpper),
              minMN, NCVT, NRU, NCC,
              S, E,
              NULL, 1, NULL, 1, NULL, 1 );
    else {
       LAPACKE_zbdsqr(
              LAPACK_COL_MAJOR, lapack_const(PlasmaLower),
              minMN, NCVT, NRU, NCC,
              S, E,
              NULL, 1, NULL, 1, NULL, 1 );
    }

    plasma_shared_free(plasma, E);

    return PLASMA_SUCCESS;
}