core_cgbelr.c

Go to the documentation of this file.

#include <lapacke.h>
#include "common.h"
/***************************************************************************/
/***************************************************************************/
#define A(_m, _n)  (PLASMA_Complex32_t *)plasma_geteltaddr(A, ((_m)-1), ((_n)-1), eltsize)
#define V(_m)      &(V[(_m)-1])
#define TAU(_m)    &(TAU[(_m)-1])
int 
CORE_cgbelr(int uplo, int N,
            PLASMA_desc *A,
            PLASMA_Complex32_t *V,
            PLASMA_Complex32_t *TAU,
            int st,
            int ed,
            int eltsize)
{
    int    NB, J1, J2;
    int    len1, len2, t1ed, t2st;
    int    i;
    static PLASMA_Complex32_t zzero = 0.0;
    PLASMA_desc vA=*A;

    /* Check input arguments */
    if (N < 0) {
        coreblas_error(2, "Illegal value of N");
        return -2;
    }
    if (ed <= st) {
        coreblas_error(6, "Illegal value of st and ed (internal)");
        return -6;
    }

    /* Quick return */
    if (N == 0)
        return PLASMA_SUCCESS;

    NB = A->mb;
    if( uplo == PlasmaLower ){
        /* ========================
         *       LOWER CASE
         * ========================*/
        for (i = ed; i >= st+1 ; i--){
            /* generate Householder to annihilate a(i+k-1,i) within the band*/
            *V(i)          = *A(i, (st-1));
            *A(i, (st-1))  = zzero;
            LAPACKE_clarfg_work( 2, A((i-1),(st-1)), V(i), 1, TAU(i));

            /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
            J1    = st;
            J2    = i-2;
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1; 
            len2  = J2-t2st+1;
            if(len1>0)CORE_clarfx2(PlasmaLeft, len1 , *V(i), conjf(*TAU(i)), A(i-1, J1  ), ELTLDD(vA, (i-1)),  A(i,  J1 ), ELTLDD(vA, i) );
            if(len2>0)CORE_clarfx2(PlasmaLeft, len2 , *V(i), conjf(*TAU(i)), A(i-1, t2st), ELTLDD(vA, (i-1)),  A(i, t2st), ELTLDD(vA, i) );
            CORE_clarfx2ce(PlasmaLower, V(i), TAU(i), A(i-1,i-1), A(i,i-1), A(i,i)); 
        }
        /* APPLY RIGHT ON THE REMAINING ELEMENT OF KERNEL 1 */
        for (i = ed; i >= st+1 ; i--){
            J1    = i+1;
            J2    = min(ed,N);
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1; 
            len2  = J2-t2st+1;
            if(len1>0)CORE_clarfx2(PlasmaRight, len1, conjf(*V(i)), conjf(*TAU(i)), A(J1,i-1),   ELTLDD(vA, J1)  , A(J1  , i), ELTLDD(vA, J1)   );
            if(len2>0)CORE_clarfx2(PlasmaRight, len2, conjf(*V(i)), conjf(*TAU(i)), A(t2st,i-1), ELTLDD(vA, t2st), A(t2st, i), ELTLDD(vA, t2st) );
        }
    } else {
        /* ========================
         *       UPPER CASE
         * ========================*/
        for (i = ed; i >= st+1 ; i--){
            /* generate Householder to annihilate a(i+k-1,i) within the band*/
            *V(i)          = *A((st-1),  i);
            *A((st-1),  i) = zzero;
            LAPACKE_clarfg_work( 2, A((st-1), (i-1)), V(i), 1, TAU(i));

            /* apply reflector from the left (horizontal row) and from the right for only the diagonal 2x2.*/
            J1    = st;
            J2    = i-2;
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1; 
            len2  = J2-t2st+1;
            if(len1>0)CORE_clarfx2(PlasmaRight, len1, conjf(*V(i)), conjf(*TAU(i)), A(J1,i-1),   ELTLDD(vA, J1)  , A(J1  , i), ELTLDD(vA, J1)   );
            if(len2>0)CORE_clarfx2(PlasmaRight, len2, conjf(*V(i)), conjf(*TAU(i)), A(t2st,i-1), ELTLDD(vA, t2st), A(t2st, i), ELTLDD(vA, t2st) );
            CORE_clarfx2ce(PlasmaUpper, V(i), TAU(i), A((i-1),(i-1)), A((i-1), i), A(i,i));
        }
        /* APPLY LEFT ON THE REMAINING ELEMENT OF KERNEL 1*/
        for (i = ed; i >= st+1 ; i--){
            J1    = i+1;
            J2    = min(ed,N);
            t1ed  = (J2/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1; 
            len2  = J2-t2st+1;
            if(len1>0)CORE_clarfx2(PlasmaLeft, len1 , *V(i), conjf(*TAU(i)), A(i-1, J1  ), ELTLDD(vA, (i-1)),  A(i,  J1 ), ELTLDD(vA, i) );
            if(len2>0)CORE_clarfx2(PlasmaLeft, len2 , *V(i), conjf(*TAU(i)), A(i-1, t2st), ELTLDD(vA, (i-1)),  A(i, t2st), ELTLDD(vA, i) );
        }
    }  /* end of else for the upper case*/

    return PLASMA_SUCCESS;

}