core_cherfb.c

Go to the documentation of this file.

#include <lapacke.h>
#include "common.h"
#undef REAL
#define COMPLEX

/***************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cherfb = PCORE_cherfb
#define CORE_cherfb PCORE_cherfb
#define CORE_cunmlq PCORE_cunmlq
#define CORE_cunmqr PCORE_cunmqr
int  CORE_cunmlq(int side, int trans,
                 int M, int N, int IB, int K,
                 PLASMA_Complex32_t *V, int LDV,
                 PLASMA_Complex32_t *T, int LDT,
                 PLASMA_Complex32_t *C, int LDC,
                 PLASMA_Complex32_t *WORK, int LDWORK);
int  CORE_cunmqr(int side, int trans,
                 int M, int N, int K, int IB,
                 PLASMA_Complex32_t *V, int LDV,
                 PLASMA_Complex32_t *T, int LDT,
                 PLASMA_Complex32_t *C, int LDC,
                 PLASMA_Complex32_t *WORK, int LDWORK);
#endif
int CORE_cherfb( PLASMA_enum uplo, int n,
                 int k, int ib, int nb,
                 PLASMA_Complex32_t *A, int lda,
                 PLASMA_Complex32_t *T, int ldt,
                 PLASMA_Complex32_t *C, int ldc,
                 PLASMA_Complex32_t *WORK, int ldwork )
{
    int i, j;

    if (uplo == PlasmaLower) {
        /* Rebuild the symmetric block: WORK <- C */
        for (j = 0; j < n; j++)
            for (i = j; i < n; i++){
                *(WORK + i + j * ldwork) = *(C + i + j*ldc);
                if (i > j){
                    *(WORK + j + i * ldwork) =  *(WORK + i + j * ldwork);
#ifdef COMPLEX
                    LAPACKE_clacgv_work(1, WORK + j + i * ldwork, ldwork);
#endif
                }
            }
        
        /* Left */
        CORE_cunmqr(PlasmaLeft, PlasmaConjTrans, n, n, k, ib, 
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);
        /* Right */
        CORE_cunmqr(PlasmaRight, PlasmaNoTrans, n, n, k, ib, 
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);
        
        /* 
         * Copy back the final result to the lower part of C 
         */
        /* C = WORK */
        for (j = 0; j < n; j++)
            for (i = j; i < n; i++)
                *(C + i + j*ldc) = *(WORK + i + j * ldwork);
    }
    else {
        /* Rebuild the symmetric block: WORK <- C */
        for (i = 0; i < n; i++)
            for (j = i; j < n; j++){
                *(WORK + i + j * ldwork) = *(C + i + j*ldc);
                if (j > i){
                    *(WORK + j + i * ldwork) =  *(WORK + i + j * ldwork);
#ifdef COMPLEX
                    LAPACKE_clacgv_work(1, WORK + j + i * ldwork, ldwork);
#endif
                }
            }
        
        /* Right */
        CORE_cunmlq(PlasmaRight, PlasmaConjTrans, n, n, k, ib, 
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);
        /* Left */
        CORE_cunmlq(PlasmaLeft, PlasmaNoTrans, n, n, k, ib, 
                    A, lda, T, ldt, WORK, ldwork, WORK+nb*ldwork, ldwork);
        
        /* 
         * Copy back the final result to the upper part of C 
         */
        /* C = WORK */
        for (i = 0; i < n; i++)
            for (j = i; j < n; j++)
                *(C + i + j*ldc) = *(WORK + i + j * ldwork);
    }
    return 0;
}


/***************************************************************************/
void QUARK_CORE_cherfb(Quark *quark, Quark_Task_Flags *task_flags,
                       PLASMA_enum uplo,
                       int n, int k, int ib, int nb,
                       PLASMA_Complex32_t *A, int lda,
                       PLASMA_Complex32_t *T, int ldt,
                       PLASMA_Complex32_t *C, int ldc)
{
    QUARK_Insert_Task(
        quark, CORE_cherfb_quark, task_flags,
        sizeof(PLASMA_enum),                     &uplo,  VALUE,
        sizeof(int),                             &n,     VALUE,
        sizeof(int),                             &k,     VALUE,
        sizeof(int),                             &ib,    VALUE,
        sizeof(int),                             &nb,    VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,        A,          uplo == PlasmaUpper ? INOUT|QUARK_REGION_U : INOUT|QUARK_REGION_L,
        sizeof(int),                             &lda,   VALUE,
        sizeof(PLASMA_Complex32_t)*ib*nb,        T,          INPUT,
        sizeof(int),                             &ldt,   VALUE,
        sizeof(PLASMA_Complex32_t)*nb*nb,        C,          uplo == PlasmaUpper ? INOUT|QUARK_REGION_D|QUARK_REGION_U : INOUT|QUARK_REGION_D|QUARK_REGION_L,
        sizeof(int),                             &ldc,   VALUE,
        sizeof(PLASMA_Complex32_t)*2*nb*nb,    NULL,         SCRATCH,
        sizeof(int),                             &nb,    VALUE,
        0);
}

/***************************************************************************/
#if defined(PLASMA_HAVE_WEAK)
#pragma weak CORE_cherfb_quark = PCORE_cherfb_quark
#define CORE_cherfb_quark PCORE_cherfb_quark
#endif
void CORE_cherfb_quark(Quark *quark)
{
    PLASMA_enum uplo;
    int n;
    int k;
    int ib;
    int nb;
    PLASMA_Complex32_t *A;
    int lda;
    PLASMA_Complex32_t *T;
    int ldt;
    PLASMA_Complex32_t *C;
    int ldc;
    PLASMA_Complex32_t *WORK;
    int ldwork;

    quark_unpack_args_13(quark, uplo, n, k, ib, nb, A, lda, T, ldt, C, ldc, WORK, ldwork);
    CORE_cherfb(uplo, n, k, ib, nb, A, lda, T, ldt, C, ldc, WORK, ldwork);
}