C++ TAU函数代码示例

void   Flow::CorrectFlow(FP T, FP p, FP ref_val, FixedValue fv) {
    FP res_p = 1.0, res_t = 1.0;
    int iter=0;
    if(fv == FV_MACH) {
       do {
           MACH(ref_val);
           t0 = T/TAU();
           p0 = p/PF();
           res_p = fabs((p0-p/PF())/p0);
           res_t = fabs((t0-T/TAU())/t0);
           Wg(ref_val*Asound());
           iter++;
         } while ((res_p > 0.0001 || res_t > 0.0001) && iter < 100);
       /*
       MACH(ref_val);
       t0 = T/TAU();
       p0 = p/PF();
       */
    } else if(fv == FV_VELOCITY) {
       do {
           MACH(ref_val/Asound());
           t0 = T/TAU();
           p0 = p/PF();
           res_p = fabs((p0-p/PF())/p0);
           res_t = fabs((t0-T/TAU())/t0);
           Wg(ref_val);
           iter++;
         } while ((res_p > 0.0001 || res_t > 0.0001) && iter < 100);
    }
}

extern "C" void
magma_ctrdtype1cbHLsym_withQ(
    magma_int_t N, magma_int_t NB,
    magmaFloatComplex *A, magma_int_t LDA,
    magmaFloatComplex *V, magmaFloatComplex *TAU,
    magma_int_t st, magma_int_t ed, magma_int_t sweep, magma_int_t Vblksiz)
{
    //magma_int_t    J1, J2, J3, i, j;
    magma_int_t    len, LDX;
    magma_int_t    IONE=1;
    magma_int_t    blkid, vpos, taupos, tpos;
    //magmaFloatComplex conjtmp;
    magmaFloatComplex Z_ONE  =  MAGMA_C_ONE;
    magmaFloatComplex *WORK;
    magma_cmalloc_cpu( &WORK, N );
    
    
    findVTpos(N,NB,Vblksiz,sweep-1,st-1, &vpos, &taupos, &tpos, &blkid);
    //printf("voici vpos %d taupos %d  tpos %d  blkid %d \n", vpos, taupos, tpos, blkid);
    LDX     = LDA-1;
    len     = ed-st+1;
    *V(vpos)  = Z_ONE;
    memcpy(V(vpos+1), A(st+1, st-1), (len-1)*sizeof(magmaFloatComplex));
    memset(A(st+1, st-1), 0, (len-1)*sizeof(magmaFloatComplex));
    /* Eliminate the col  at st-1 */
    lapackf77_clarfg( &len, A(st, st-1), V(vpos+1), &IONE, TAU(taupos) );
    /* apply left and right on A(st:ed,st:ed)*/
    magma_clarfxsym(len,A(st,st),LDX,V(vpos),TAU(taupos));
    //conjtmp = MAGMA_C_CONJ(*TAU(taupos));
    //lapackf77_clarfy("L", &len, V(vpos), &IONE, &conjtmp, A(st,st), &LDX, WORK); //&(MAGMA_C_CONJ(*TAU(taupos)))
    magma_free_cpu(WORK);
}

extern "C" void
magma_strdtype1cbHLsym_withQ_v2(magma_int_t n, magma_int_t nb,
                                float *A, magma_int_t lda,
                                float *V, magma_int_t ldv,
                                float *TAU,
                                magma_int_t st, magma_int_t ed,
                                magma_int_t sweep, magma_int_t Vblksiz,
                                float *work)
{
/*
    WORK (workspace) float real array, dimension N
*/

    magma_int_t ione = 1;
    magma_int_t vpos, taupos, len, len2;

    float c_one    =  MAGMA_S_ONE;

    magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, st-1, ldv, &vpos, &taupos);
    //printf("voici vpos %d taupos %d  tpos %d  blkid %d \n", vpos, taupos, tpos, blkid);

    len     = ed-st+1;
    *V(vpos)  = c_one;

    len2 = len-1;
    blasf77_scopy( &len2, A(st+1, st-1), &ione, V(vpos+1), &ione );
    //memcpy(V(vpos+1), A(st+1, st-1), (len-1)*sizeof(float));
    memset(A(st+1, st-1), 0, (len-1)*sizeof(float));

    /* Eliminate the col  at st-1 */
    lapackf77_slarfg( &len, A(st, st-1), V(vpos+1), &ione, TAU(taupos) );

    /* apply left and right on A(st:ed,st:ed)*/
    magma_slarfxsym_v2(len, A(st,st), lda-1, V(vpos), TAU(taupos), work);
}

extern "C" void
magma_ctrdtype2cbHLsym_withQ_v2(
    magma_int_t n, magma_int_t nb,
    magmaFloatComplex *A, magma_int_t lda,
    magmaFloatComplex *V, magma_int_t ldv,
    magmaFloatComplex *TAU,
    magma_int_t st, magma_int_t ed,
    magma_int_t sweep, magma_int_t Vblksiz,
    magmaFloatComplex *work)
{
    /*
     WORK (workspace) float complex array, dimension NB
    */

    magma_int_t ione = 1;
    magma_int_t vpos, taupos;

    magmaFloatComplex conjtmp;

    magmaFloatComplex c_one = MAGMA_C_ONE;

    magma_int_t ldx = lda-1;
    magma_int_t len = ed - st + 1;
    magma_int_t lem = min(ed+nb, n) - ed;
    magma_int_t lem2;
    
    if (lem > 0) {
        magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, st-1, ldv, &vpos, &taupos);
        /* apply remaining right coming from the top block */
        lapackf77_clarfx("R", &lem, &len, V(vpos), TAU(taupos), A(ed+1, st), &ldx, work);
    }
    if (lem > 1) {
        magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, ed, ldv, &vpos, &taupos);

        /* remove the first column of the created bulge */
        *V(vpos)  = c_one;
        //memcpy(V(vpos+1), A(ed+2, st), (lem-1)*sizeof(magmaFloatComplex));
        lem2 = lem-1;
        blasf77_ccopy( &lem2, A(ed+2, st), &ione, V(vpos+1), &ione );
        memset(A(ed+2, st),0,(lem-1)*sizeof(magmaFloatComplex));

        /* Eliminate the col at st */
        lapackf77_clarfg( &lem, A(ed+1, st), V(vpos+1), &ione, TAU(taupos) );

        /* apply left on A(J1:J2,st+1:ed) */
        len = len-1; /* because we start at col st+1 instead of st. col st is the col that has been removed; */
        conjtmp = MAGMA_C_CNJG(*TAU(taupos));
        lapackf77_clarfx("L", &lem, &len, V(vpos),  &conjtmp, A(ed+1, st+1), &ldx, work);
    }
}

extern "C" void
magma_ctrdtype3cbHLsym_withQ(
    magma_int_t N, magma_int_t NB,
    magmaFloatComplex *A, magma_int_t LDA,
    magmaFloatComplex *V, magmaFloatComplex *TAU,
    magma_int_t st, magma_int_t ed, magma_int_t sweep, magma_int_t Vblksiz)
{
    //magma_int_t    J1, J2, J3, i, j;
    magma_int_t    len, LDX;
    //magma_int_t    IONE=1;
    magma_int_t    blkid, vpos, taupos, tpos;
    //magmaFloatComplex conjtmp;
    magmaFloatComplex *WORK;
    magma_cmalloc_cpu( &WORK, N );
    
    
    findVTpos(N,NB,Vblksiz,sweep-1,st-1, &vpos, &taupos, &tpos, &blkid);
    LDX    = LDA-1;
    len    = ed-st+1;
    
    /* apply left and right on A(st:ed,st:ed)*/
    magma_clarfxsym(len,A(st,st),LDX,V(vpos),TAU(taupos));
    //conjtmp = MAGMA_C_CONJ(*TAU(taupos));
    //lapackf77_clarfy("L", &len, V(vpos), &IONE,  &(MAGMA_C_CONJ(*TAU(taupos))), A(st,st), &LDX, WORK);
    magma_free_cpu(WORK);
}

/***************************************************************************//**
 *          TYPE 3-BAND Lower-columnwise-Householder
 ***************************************************************************/
extern "C" void
magma_dsbtype3cb(magma_int_t n, magma_int_t nb,
                double *A, magma_int_t lda,
                double *V, magma_int_t ldv, 
                double *TAU,
                magma_int_t st, magma_int_t ed, magma_int_t sweep, 
                magma_int_t Vblksiz, magma_int_t wantz,
                double *work)
{
    magma_int_t len;
    magma_int_t vpos, taupos;
    //magma_int_t blkid, tpos;


    if ( wantz == 0 ) {
        vpos   = (sweep%2)*n + st;
        taupos = (sweep%2)*n + st;
    } else {
        magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep, st, ldv, &vpos, &taupos);
        //findVTpos(n, nb, Vblksiz, sweep, st, &vpos, &taupos, &tpos, &blkid);
    }

    len = ed-st+1;
    /* Apply left and right on A(st:ed,st:ed)*/
    magma_dlarfy(len, A(st,st), lda-1, V(vpos), TAU(taupos), work);
    return;
}

static void magma_stile_bulge_computeT_parallel(magma_int_t my_core_id, magma_int_t cores_num, float *V, magma_int_t ldv, float *TAU,
                                                float *T, magma_int_t ldt, magma_int_t n, magma_int_t nb, magma_int_t Vblksiz)
{
    //%===========================
    //%   local variables
    //%===========================
    magma_int_t firstcolj;
    magma_int_t rownbm;
    magma_int_t st,ed,fst,vlen,vnb,colj;
    magma_int_t blkid,vpos,taupos,tpos;
    magma_int_t blkpercore, myid;
    
    if(n<=0)
        return ;
    
    magma_int_t blkcnt = magma_bulge_get_blkcnt(n, nb, Vblksiz);
    
    blkpercore = blkcnt/cores_num;
    
    magma_int_t nbGblk  = magma_ceildiv(n-1, Vblksiz);
    
    if(my_core_id==0) printf("  COMPUTE T parallel threads %d with  N %d   NB %d   Vblksiz %d \n",cores_num,n,nb,Vblksiz);
    
    for (magma_int_t bg = nbGblk; bg>0; bg--)
    {
        firstcolj = (bg-1)*Vblksiz + 1;
        rownbm    = magma_ceildiv(n-(firstcolj+1), nb);
        if(bg==nbGblk) 
            rownbm    = magma_ceildiv(n-firstcolj ,nb);  // last blk has size=1 used for real to handle A(N,N-1)

        for (magma_int_t m = rownbm; m>0; m--)
        {
            vlen = 0;
            vnb  = 0;
            colj      = (bg-1)*Vblksiz; // for k=0;I compute the fst and then can remove it from the loop
            fst       = (rownbm -m)*nb+colj +1;
            for (magma_int_t k=0; k<Vblksiz; k++)
            {
                colj     = (bg-1)*Vblksiz + k;
                st       = (rownbm -m)*nb+colj +1;
                ed       = min(st+nb-1,n-1);
                if(st>ed)
                    break;
                if((st==ed)&&(colj!=n-2))
                    break;
                
                vlen=ed-fst+1;
                vnb=k+1;
            }        
            colj     = (bg-1)*Vblksiz;
            magma_bulge_findVTAUTpos(n, nb, Vblksiz, colj, fst, ldv, ldt, &vpos, &taupos, &tpos, &blkid);
            myid = blkid/blkpercore;
            if(my_core_id==(myid%cores_num)){
                if((vlen>0)&&(vnb>0))
                    lapackf77_slarft( "F", "C", &vlen, &vnb, V(vpos), &ldv, TAU(taupos), T(tpos), &ldt);
            }
        }
    }
}

extern "C" void
magma_ctrdtype2cbHLsym_withQ(
    magma_int_t N, magma_int_t NB,
    magmaFloatComplex *A, magma_int_t LDA,
    magmaFloatComplex *V, magmaFloatComplex *TAU,
    magma_int_t st, magma_int_t ed, magma_int_t sweep, magma_int_t Vblksiz)
{
    magma_int_t    J1, J2, len, lem, LDX;
    //magma_int_t    i, j;
    magma_int_t    IONE=1;
    magma_int_t    blkid, vpos, taupos, tpos;
    magmaFloatComplex conjtmp;
    magmaFloatComplex Z_ONE  =  MAGMA_C_ONE;
    //magmaFloatComplex WORK[NB];
    magmaFloatComplex *WORK;
    magma_cmalloc_cpu( &WORK, NB );
    
    
    findVTpos(N,NB,Vblksiz,sweep-1,st-1, &vpos, &taupos, &tpos, &blkid);
    LDX    = LDA-1;
    J1     = ed+1;
    J2     = min(ed+NB,N);
    len    = ed-st+1;
    lem    = J2-J1+1;
    if (lem > 0) {
        /* apply remaining right commming from the top block */
        lapackf77_clarfx("R", &lem, &len, V(vpos), TAU(taupos), A(J1, st), &LDX, WORK);
    }
    if (lem > 1) {
        findVTpos(N,NB,Vblksiz,sweep-1,J1-1, &vpos, &taupos, &tpos, &blkid);
        /* remove the first column of the created bulge */
        *V(vpos)  = Z_ONE;
        memcpy(V(vpos+1), A(J1+1, st), (lem-1)*sizeof(magmaFloatComplex));
        memset(A(J1+1, st),0,(lem-1)*sizeof(magmaFloatComplex));
        /* Eliminate the col at st */
        lapackf77_clarfg( &lem, A(J1, st), V(vpos+1), &IONE, TAU(taupos) );
        /* apply left on A(J1:J2,st+1:ed) */
        len = len-1; /* because we start at col st+1 instead of st. col st is the col that has been revomved; */
        conjtmp = MAGMA_C_CONJ(*TAU(taupos));
        lapackf77_clarfx("L", &lem, &len, V(vpos),  &conjtmp, A(J1, st+1), &LDX, WORK);
    }
    magma_free_cpu(WORK);
}

extern "C" void
magma_dsbtype1cb(magma_int_t n, magma_int_t nb,
                double *A, magma_int_t lda,
                double *V, magma_int_t ldv, 
                double *TAU,
                magma_int_t st, magma_int_t ed, magma_int_t sweep, 
                magma_int_t Vblksiz, magma_int_t wantz,
                double *work)
{
    magma_int_t len;
    magma_int_t vpos, taupos;
    //magma_int_t blkid, tpos;
    
    magma_int_t ione = 1;
    double c_one    =  MAGMA_D_ONE;

    /* find the pointer to the Vs and Ts as stored by the bulgechasing
     * note that in case no eigenvector required V and T are stored
     * on a vector of size n
     * */
     if ( wantz == 0 ) {
         vpos   = (sweep%2)*n + st;
         taupos = (sweep%2)*n + st;
     } else {
         magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep, st, ldv, &vpos, &taupos);
         //findVTpos(n, nb, Vblksiz, sweep, st, &vpos, &taupos, &tpos, &blkid);
     }

    len = ed-st+1;
    *(V(vpos)) = c_one;

    //magma_int_t len2 = len-1;
    //blasf77_dcopy( &len2, A(st+1, st-1), &ione, V(vpos+1), &ione );
    memcpy( V(vpos+1), A(st+1, st-1), (len-1)*sizeof(double) );
    memset( A(st+1, st-1), 0, (len-1)*sizeof(double) );

    /* Eliminate the col  at st-1 */
    lapackf77_dlarfg( &len, A(st, st-1), V(vpos+1), &ione, TAU(taupos) );

    /* Apply left and right on A(st:ed,st:ed) */
    magma_dlarfy(len, A(st,st), lda-1, V(vpos), TAU(taupos), work);
}

extern "C" void
magma_ztrdtype1cbHLsym_withQ_v2(magma_int_t n, magma_int_t nb, 
                                magmaDoubleComplex *A, magma_int_t lda, 
                                magmaDoubleComplex *V, magma_int_t ldv, 
                                magmaDoubleComplex *TAU,
                                magma_int_t st, magma_int_t ed, 
                                magma_int_t sweep, magma_int_t Vblksiz, 
                                magmaDoubleComplex *work) 
{

/*
    WORK (workspace) double complex array, dimension N
*/

    magma_int_t ione = 1;
    magma_int_t vpos, taupos, len;

    magmaDoubleComplex c_one    =  MAGMA_Z_ONE;

    magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, st-1, ldv, &vpos, &taupos);
    //printf("voici vpos %d taupos %d  tpos %d  blkid %d \n", vpos, taupos, tpos, blkid);

    len     = ed-st+1;
    *V(vpos)  = c_one;

    cblas_zcopy(len-1, A(st+1, st-1), ione, V(vpos+1), ione);
    //memcpy(V(vpos+1), A(st+1, st-1), (len-1)*sizeof(magmaDoubleComplex));
    memset(A(st+1, st-1), 0, (len-1)*sizeof(magmaDoubleComplex));

    /* Eliminate the col  at st-1 */
    lapackf77_zlarfg( &len, A(st, st-1), V(vpos+1), &ione, TAU(taupos) );

    /* apply left and right on A(st:ed,st:ed)*/
    magma_zlarfxsym_v2(len, A(st,st), lda-1, V(vpos), TAU(taupos), work);

}

int SpIO_H5ReadTau(hid_t h5f_id, Zone *zone)
/* Write data of all children to an HDF5 table */
{
	SpPhys *pp = zone->data;

	/* Just in case the programmer did something stupid */
	Deb_ASSERT(pp->mol != NULL);
	Deb_ASSERT(pp->tau != NULL);

	int status = 0;
	size_t
		i, j,
		nrad = pp->mol->nrad,
		record_size =  sizeof(double) * nrad,
		field_sizes[nrad],
		field_offsets[nrad];
	double *tau = Mem_CALLOC(zone->nchildren * nrad, tau);
	herr_t hstatus;

	/* Init fields */
	for(i = 0; i < nrad; i++) {
		field_offsets[i] = i * sizeof(double);
		field_sizes[i] = sizeof(double);
	}

	/* Read tau table */
	hstatus = H5TBread_table(h5f_id, "TAU", record_size, field_offsets, field_sizes, tau);
	if(hstatus < 0)
		status = Err_SETSTRING("Error reading HDF5 `%s' table", "TAU");

	#define TAU(i, j)\
		tau[(j) + nrad * (i)]

	for(i = 0; i < zone->nchildren; i++) {
		pp = zone->children[i]->data;

		for(j = 0; j < nrad; j++) {
			pp->tau[j] = TAU(i, j);
		}
	}

	#undef TAU

	free(tau);

	return status;
}

extern "C" void
magma_strdtype3cbHLsym_withQ_v2(magma_int_t n, magma_int_t nb,
                                float *A, magma_int_t lda,
                                float *V, magma_int_t ldv,
                                float *TAU,
                                magma_int_t st, magma_int_t ed,
                                magma_int_t sweep, magma_int_t Vblksiz,
                                float *work)
{
    /*
     WORK (workspace) float real array, dimension N
     */

    magma_int_t vpos, taupos;

    magma_bulge_findVTAUpos(n, nb, Vblksiz, sweep-1, st-1, ldv, &vpos, &taupos);

    magma_int_t len = ed-st+1;

    /* apply left and right on A(st:ed,st:ed)*/
    magma_slarfxsym_v2(len, A(st,st), lda-1, V(vpos), TAU(taupos), work);
}

/* Subroutine */ int dormbr_(char *vect, char *side, char *trans, integer *m, 
	integer *n, integer *k, doublereal *a, integer *lda, doublereal *tau, 
	doublereal *c, integer *ldc, doublereal *work, integer *lwork, 
	integer *info)
{
/*  -- LAPACK routine (version 2.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       September 30, 1994   


    Purpose   
    =======   

    If VECT = 'Q', DORMBR overwrites the general real M-by-N matrix C   
    with   
                    SIDE = 'L'     SIDE = 'R'   
    TRANS = 'N':      Q * C          C * Q   
    TRANS = 'T':      Q**T * C       C * Q**T   

    If VECT = 'P', DORMBR overwrites the general real M-by-N matrix C   
    with   
                    SIDE = 'L'     SIDE = 'R'   
    TRANS = 'N':      P * C          C * P   
    TRANS = 'T':      P**T * C       C * P**T   

    Here Q and P**T are the orthogonal matrices determined by DGEBRD when 
  
    reducing a real matrix A to bidiagonal form: A = Q * B * P**T. Q and 
  
    P**T are defined as products of elementary reflectors H(i) and G(i)   
    respectively.   

    Let nq = m if SIDE = 'L' and nq = n if SIDE = 'R'. Thus nq is the   
    order of the orthogonal matrix Q or P**T that is applied.   

    If VECT = 'Q', A is assumed to have been an NQ-by-K matrix:   
    if nq >= k, Q = H(1) H(2) . . . H(k);   
    if nq < k, Q = H(1) H(2) . . . H(nq-1).   

    If VECT = 'P', A is assumed to have been a K-by-NQ matrix:   
    if k < nq, P = G(1) G(2) . . . G(k);   
    if k >= nq, P = G(1) G(2) . . . G(nq-1).   

    Arguments   
    =========   

    VECT    (input) CHARACTER*1   
            = 'Q': apply Q or Q**T;   
            = 'P': apply P or P**T.   

    SIDE    (input) CHARACTER*1   
            = 'L': apply Q, Q**T, P or P**T from the Left;   
            = 'R': apply Q, Q**T, P or P**T from the Right.   

    TRANS   (input) CHARACTER*1   
            = 'N':  No transpose, apply Q  or P;   
            = 'T':  Transpose, apply Q**T or P**T.   

    M       (input) INTEGER   
            The number of rows of the matrix C. M >= 0.   

    N       (input) INTEGER   
            The number of columns of the matrix C. N >= 0.   

    K       (input) INTEGER   
            If VECT = 'Q', the number of columns in the original   
            matrix reduced by DGEBRD.   
            If VECT = 'P', the number of rows in the original   
            matrix reduced by DGEBRD.   
            K >= 0.   

    A       (input) DOUBLE PRECISION array, dimension   
                                  (LDA,min(nq,K)) if VECT = 'Q'   
                                  (LDA,nq)        if VECT = 'P'   
            The vectors which define the elementary reflectors H(i) and   
            G(i), whose products determine the matrices Q and P, as   
            returned by DGEBRD.   

    LDA     (input) INTEGER   
            The leading dimension of the array A.   
            If VECT = 'Q', LDA >= max(1,nq);   
            if VECT = 'P', LDA >= max(1,min(nq,K)).   

    TAU     (input) DOUBLE PRECISION array, dimension (min(nq,K))   
            TAU(i) must contain the scalar factor of the elementary   
            reflector H(i) or G(i) which determines Q or P, as returned   
            by DGEBRD in the array argument TAUQ or TAUP.   

    C       (input/output) DOUBLE PRECISION array, dimension (LDC,N)   
            On entry, the M-by-N matrix C.   
            On exit, C is overwritten by Q*C or Q**T*C or C*Q**T or C*Q   
            or P*C or P**T*C or C*P or C*P**T.   

    LDC     (input) INTEGER   
            The leading dimension of the array C. LDC >= max(1,M).   

    WORK    (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 
  
            On exit, if INFO = 0, WORK(1) returns the optimal LWORK.   
//.........这里部分代码省略.........

/* Subroutine */ int ssytd2_(char *uplo, integer *n, real *a, integer *lda, 
	real *d, real *e, real *tau, integer *info)
{
/*  -- LAPACK routine (version 2.0) --   
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
       Courant Institute, Argonne National Lab, and Rice University   
       October 31, 1992   


    Purpose   
    =======   

    SSYTD2 reduces a real symmetric matrix A to symmetric tridiagonal   
    form T by an orthogonal similarity transformation: Q' * A * Q = T.   

    Arguments   
    =========   

    UPLO    (input) CHARACTER*1   
            Specifies whether the upper or lower triangular part of the   
            symmetric matrix A is stored:   
            = 'U':  Upper triangular   
            = 'L':  Lower triangular   

    N       (input) INTEGER   
            The order of the matrix A.  N >= 0.   

    A       (input/output) REAL array, dimension (LDA,N)   
            On entry, the symmetric matrix A.  If UPLO = 'U', the leading 
  
            n-by-n upper triangular part of A contains the upper   
            triangular part of the matrix A, and the strictly lower   
            triangular part of A is not referenced.  If UPLO = 'L', the   
            leading n-by-n lower triangular part of A contains the lower 
  
            triangular part of the matrix A, and the strictly upper   
            triangular part of A is not referenced.   
            On exit, if UPLO = 'U', the diagonal and first superdiagonal 
  
            of A are overwritten by the corresponding elements of the   
            tridiagonal matrix T, and the elements above the first   
            superdiagonal, with the array TAU, represent the orthogonal   
            matrix Q as a product of elementary reflectors; if UPLO   
            = 'L', the diagonal and first subdiagonal of A are over-   
            written by the corresponding elements of the tridiagonal   
            matrix T, and the elements below the first subdiagonal, with 
  
            the array TAU, represent the orthogonal matrix Q as a product 
  
            of elementary reflectors. See Further Details.   

    LDA     (input) INTEGER   
            The leading dimension of the array A.  LDA >= max(1,N).   

    D       (output) REAL array, dimension (N)   
            The diagonal elements of the tridiagonal matrix T:   
            D(i) = A(i,i).   

    E       (output) REAL array, dimension (N-1)   
            The off-diagonal elements of the tridiagonal matrix T:   
            E(i) = A(i,i+1) if UPLO = 'U', E(i) = A(i+1,i) if UPLO = 'L'. 
  

    TAU     (output) REAL array, dimension (N-1)   
            The scalar factors of the elementary reflectors (see Further 
  
            Details).   

    INFO    (output) INTEGER   
            = 0:  successful exit   
            < 0:  if INFO = -i, the i-th argument had an illegal value.   

    Further Details   
    ===============   

    If UPLO = 'U', the matrix Q is represented as a product of elementary 
  
    reflectors   

       Q = H(n-1) . . . H(2) H(1).   

    Each H(i) has the form   

       H(i) = I - tau * v * v'   

    where tau is a real scalar, and v is a real vector with   
    v(i+1:n) = 0 and v(i) = 1; v(1:i-1) is stored on exit in   
    A(1:i-1,i+1), and tau in TAU(i).   

    If UPLO = 'L', the matrix Q is represented as a product of elementary 
  
    reflectors   

       Q = H(1) H(2) . . . H(n-1).   

    Each H(i) has the form   

       H(i) = I - tau * v * v'   

    where tau is a real scalar, and v is a real vector with   
//.........这里部分代码省略.........

int 
CORE_shbrce(int uplo, int N,
            PLASMA_desc *A,
            float *V,
            float *TAU,
            int st,
            int ed,
            int eltsize)
{
    int    NB, J1, J2, J3, KDM2, len, pt;
    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(6, "Illegal value of st and ed (internal)");
        return -6;
    }

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

    NB = A->mb;
    KDM2 = A->mb-2;
    if( uplo == PlasmaLower ) {
        /* ========================
         *       LOWER CASE
         * ========================*/
        for (i = ed; i >= st+1 ; i--){
            /* apply Householder from the right. and create newnnz outside the band if J3 < N */
            J1  = ed+1;
            J2  = min((i+1+KDM2), N);
            J3  = min((J2+1), N);
            len = J3-J1+1;
            if(J3>J2)*A(J3,(i-1))=zzero;/* could be removed because A is supposed to be band.*/

            t1ed  = (J3/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1; /* can be negative*/
            len2  = J3-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));
            /* if nonzero element a(j+kd,j-1) has been created outside the band (if index < N) then eliminate it.*/
            len    = J3-J2; // soit 1 soit 0
            if(len>0){
                /* generate Householder to annihilate a(j+kd,j-1) within the band */
                *V(J3)         = *A(J3,i-1);
                *A(J3,i-1)   = 0.0;
                LAPACKE_slarfg_work( 2, A(J2,i-1), V(J3), 1, TAU(J3));
            }
        }
        /* APPLY LEFT ON THE REMAINING ELEMENT OF KERNEL 2 */
        for (i = ed; i >= st+1 ; i--){
            /* find if there was a nnz created. if yes apply left else nothing to be done.*/
            J2  = min((i+1+KDM2), N);
            J3  = min((J2+1), N);
            len    = J3-J2;
            if(len>0){
                pt    = J2;
                J1    = i;
                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(J3), (*TAU(J3)), A(pt, i   ), ELTLDD(vA, pt),  A((pt+1),  i  ), ELTLDD(vA, pt+1) );
                if(len2>0)CORE_slarfx2(PlasmaLeft, len2 , *V(J3), (*TAU(J3)), A(pt, t2st), ELTLDD(vA, pt),  A((pt+1), t2st), ELTLDD(vA, pt+1) );
            }
        }
    } else {
        /* ========================
         *       UPPER CASE
         * ========================*/
        for (i = ed; i >= st+1 ; i--){
            /* apply Householder from the right. and create newnnz outside the band if J3 < N */
            J1  = ed+1;
            J2  = min((i+1+KDM2), N);
            J3  = min((J2+1), N);
            len = J3-J1+1;
            if(J3>J2)*A((i-1), J3)=zzero;/* could be removed because A is supposed to be band.*/

            t1ed  = (J3/NB)*NB;
            t2st  = max(t1ed+1,J1);
            len1  = t1ed-J1+1;  /* can be negative*/
            len2  = J3-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) );
            /* if nonzero element a(j+kd,j-1) has been created outside the band (if index < N) then eliminate it.*/
            len    = J3-J2; /* either 1 soit 0*/
            if(len>0){
                /* generate Householder to annihilate a(j+kd,j-1) within the band*/
                *V(J3)         = *A((i-1), J3);
//.........这里部分代码省略.........