core_shbelr.c

Go to the documentation of this file.

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

/***************************************************************************/
/***************************************************************************/
#define A(_m, _n)  (float *)plasma_geteltaddr(A, ((_m)-1), ((_n)-1), eltsize)
#define V(_m)      &(V[(_m)-1])
#define TAU(_m)    &(TAU[(_m)-1])
int 
CORE_shbelr(int uplo, int N,
            PLASMA_desc *A,
            float *V,
            float *TAU,
            int st,
            int ed,
            int eltsize)
{
    int    NB, J1, J2;
    int    len1, len2, t1ed, t2st;
    int    i;
    static float 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(23, "Illegal value of st and ed (internal)");
        return -23;
    }

    /* 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_slarfg_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; /* can be negative */
          len2  = J2-t2st+1;
          if(len1>0)CORE_slarfx2(PlasmaLeft, len1 , *V(i), (*TAU(i)), A(i-1, J1  ), ELTLDD(vA, i-1),  A(i,  J1 ), ELTLDD(vA, i) );
          if(len2>0)CORE_slarfx2(PlasmaLeft, len2 , *V(i), (*TAU(i)), A(i-1, t2st), ELTLDD(vA, i-1),  A(i, t2st), ELTLDD(vA, i) );
          CORE_slarfx2c(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; /* can be negative */
          len2  = J2-t2st+1;
          if(len1>0)CORE_slarfx2(PlasmaRight, len1, *V(i), *TAU(i), A(J1,  i-1), ELTLDD(vA, J1)  , A(J1  , i), ELTLDD(vA, J1)   );
          if(len2>0)CORE_slarfx2(PlasmaRight, len2, *V(i), *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_slarfg_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; /* can be negative */
          len2  = J2-t2st+1;
          if(len1>0)CORE_slarfx2(PlasmaRight, len1, (*V(i)), (*TAU(i)), A(J1,  i-1), ELTLDD(vA, J1)  , A(J1  , i), ELTLDD(vA, J1)   );
          if(len2>0)CORE_slarfx2(PlasmaRight, len2, (*V(i)), (*TAU(i)), A(t2st,i-1), ELTLDD(vA, t2st), A(t2st, i), ELTLDD(vA, t2st) );
          CORE_slarfx2c(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; /* can be negative */
          len2  = J2-t2st+1;
          if(len1>0)CORE_slarfx2(PlasmaLeft, len1 , (*V(i)), *TAU(i), A(i-1, J1  ), ELTLDD(vA, i-1),  A(i,  J1 ), ELTLDD(vA, i) );
          if(len2>0)CORE_slarfx2(PlasmaLeft, len2 , (*V(i)), *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;

}