Commit 9903873c authored by thomas dettbarn's avatar thomas dettbarn
Browse files

second iteration is working to a point.

parent 704ccb7f
Loading
Loading
Loading
Loading
+58 −58
Original line number Diff line number Diff line
@@ -902,9 +902,25 @@ static void HouseholderReduction_fx(
			nCh
		);

		{
			int i,j;
			printf("\ncompare%d, (%d)\x1b[1;32mstart\x1b[0m\n",nCh,singularVectors_Left_e);
			for (i=0;i<nChannelsL;i++)
			{
				printf("compare%d,%d: \x1b[1;32m",nCh,i);
				for (j=0;j<nChannelsC;j++)
				{
					Word16 n;
					n=W_norm(singularVectors_Left_64[i][j]);
					printf("%08X ",W_extract_h(W_shl(singularVectors_Left_64[i][j],n)));
				}
				printf("\x1b[0m\n");
			}
			printf("compare%d \x1b[1;32mend\x1b[0m\n",nCh);
		}

		biDiagonalReductionRight_64(
			singularVectors_Left_64,0,
			singularVectors_Left_64,nCh,
			secDiag_fx,secDiag_fx_e,
			nChannelsL,
			nChannelsC,
@@ -946,6 +962,21 @@ static void HouseholderReduction_fx(
    {
#ifdef	MYCHANGES
        biDiagonalReductionLeft_fx( singularVectors_Left_fx, singularValues_fx, singularVectors_Left_fx_e, singularValues_fx_e, nChannelsL, nChannelsC, nCh );
       		{
			int i,j;
			printf("compare%d, start\n",nCh);
			for (i=0;i<nChannelsL;i++)
			{
				printf("compare%d,%d: ",nCh,i);
				for (j=0;j<nChannelsC;j++)
				{
					Word16 n;
					printf("%08X ",singularVectors_Left_fx[i][j]);
				}
				printf("\n");
			}
			printf("compare%d, end\n",nCh);
		}
	biDiagonalReductionRight_fx( singularVectors_Left_fx, secDiag_fx, singularVectors_Left_fx_e, secDiag_fx_e, nChannelsL, nChannelsC, nCh, &g_fx, &g_e );
		{
			int i,j;
@@ -1087,7 +1118,6 @@ static void biDiagonalReductionLeft_64(
            norm_x=W_extract_h(W_shl(norm_64,norm_x_e));
            f=Mpy_32_32(norm_x,invVal);
            f_e=add(invVal_e, sub(norm_x_e, r_e ));
		//printf("MYADD%d,%d: [%08X<%2x]    (norm_x:%d r_e:%d  invVal_e:%d f_e:%d -->%d)  \x1b[1;33m",iCh,currChannel,f,f_e,norm_x_e,r_e,invVal_e,f_e,shifter[iCh]);
            for (jCh=currChannel;jCh<nChannelsL; jCh++)
            {

@@ -1097,10 +1127,8 @@ static void biDiagonalReductionLeft_64(
		{
			Word16 n;
			n=W_norm(singularVectors_Left_64[jCh][iCh]);
			//printf("%08X ",W_extract_h(W_shl(singularVectors_Left_64[jCh][iCh],n)));
		}
            }
		//printf("\x1b[0m\n");
        }
    }
}
@@ -1149,7 +1177,6 @@ static void biDiagonalReductionRight_64(
        move64();
        move64();
        tmpe=add(singularVectors_e,1);
	printf("NORM64: \x1b[1;33m");
        for (jCh=idx;jCh<nChannelsC; jCh++)
        {
            tmp=W_extract_l(W_shr(singularVectors_Left_64[currChannel][jCh],tmpe));
@@ -1164,7 +1191,10 @@ static void biDiagonalReductionRight_64(
        tmpe=W_norm(abs_64);
        abs_x=W_extract_h(W_shl(abs_64, tmpe ));
        abs_x_e = add(sub(shl(singularVectors_e, 1), tmpe), 3 );
	printf("%016llX --> %08X<%2X   abs_x:%016llX --> %08X<%2X\x1b[0m\n",norm_64,norm_x,norm_x_e,abs_64,abs_x,abs_x_e);
	printf("\n");
	printf("START%d: \x1b[1;37m",currChannel);
	printf("norm:%016llx  %08X<%2X   abs:%016llx %08X<%2X",norm_64,norm_x,norm_x_e, abs_64,abs_x,abs_x_e);
	printf("\x1b[0m\n");
        IF (norm_x)
        {
            Word16 invVal_e;
@@ -1197,23 +1227,30 @@ static void biDiagonalReductionRight_64(
            tmpmul=W_mult0_32_32(*g,tmp);
            tmpmul=W_shl(tmpmul,tmpe);
            r_64=W_sub(W_shr(tmpmul,1), W_shr(norm_64,1) );
		printf("MUL: \x1b[1;33ms:%016llX -> tmp:%08X -> g:%08X tmpmul:%016llX -> norm_64:%016llX --> R:%016llX\x1b[0m\n",singularVectors_Left_64[currChannel][idx],tmp,*g,tmpmul,norm_64,r_64);
            
            tmpe2=W_norm(r_64);
            r=W_extract_h(W_shl(r_64,tmpe2));
            r_e = sub(sub(add(add(1,singularVectors_e),*g_e),tmpe),tmpe2);
		printf("G%d: \x1b[1;37m %08X<%2X  R: %016llx %08X<%2X\x1b[0m\n",currChannel,*g,*g_e, r_64,r,r_e);

		printf("R: \x1b[1;33m%016llX --> %08X<%2X   g:%08X<%2X\x1b[0m\n",r_64,r,r_e,*g,*g_e);
    
            invVal_e = r_e;
            invVal = BASOP_Util_Divide3232_Scale_newton( MAXVAL_WORD32, maxWithSign_fx( r ), &invVal_e);
    
//            tmpe=add(31,sub(singularVectors_e,g_e));	// TODO: maybe the other way around??
		tmpe=sub(30,singularVectors_e);
		switch(currChannel)
		{
			case 0: tmpe=30;break;
			case 1: tmpe=33;break;
			default: tmpe=32-(*g_e);break;
		}
            tmp64=W_shr(W_deposit32_h(*g),tmpe);
		printf("SING: \x1b[1;33m%016llX-g:%08X (%016llX) = ",singularVectors_Left_64[currChannel][idx],*g,tmp64);
		printf("SING%d: \x1b[1;37m",currChannel);
		printf("%016llx-%016llx=",singularVectors_Left_64[currChannel][idx],tmp64);
    	    singularVectors_Left_64[currChannel][idx]=W_sub(singularVectors_Left_64[currChannel][idx],tmp64);	// exponent +1
		printf("%016llX\x1b[0m\n",singularVectors_Left_64[currChannel][idx]);
		printf("%016llx",singularVectors_Left_64[currChannel][idx]);
		printf("\x1b[0m   (tmpe:%d invval_e:%d norm_x_e:%d g_e:%d abs_x_e:%d\n",tmpe,invVal_e,norm_x_e, *g_e, abs_x_e);


            tmpe=add(singularVectors_e,2);
@@ -1225,58 +1262,36 @@ static void biDiagonalReductionRight_64(
                tmp2=W_extract_l(W_shr(singularVectors_Left_64[iCh][currChannel],tmpe));
		norm_64=0;
		move64();
		printf("norm%d,%d: \x1b[1;33m",currChannel,iCh);
                for (jCh=idx;jCh<nChannelsC; jCh++)
                {
                    tmp2=W_extract_l(W_shr(singularVectors_Left_64[iCh][jCh],tmpe));
                    tmp3=W_extract_l(W_shr(singularVectors_Left_64[currChannel][jCh],tmpe));
                    norm_64=W_add(norm_64, W_mult0_32_32(tmp2,tmp3));
		}
		printf("%016llx\x1b[0m\n",norm_64);
                norm_x_e=W_norm(norm_64);
                norm_x=W_extract_h(W_shl(norm_64,norm_x_e));
                f=Mpy_32_32(norm_x,invVal);
                f_e=add(invVal_e, sub(norm_x_e, r_e ));
		printf("F: \x1b[1;33mnorm:%016llX --> %08X<%2X   %08X<%2X\x1b[0m\n",norm_64,norm_x,norm_x_e,f,f_e);
   		printf("singularVectors_Left%d,%d: \x1b[1;35m",currChannel,iCh); 
		printf("f_e:d norm_x_e:%d norm_x_e2:%d tmpe:%d ",f_e,norm_x_e,norm_x_e2,tmpe);

			switch(iCh)
			{
				case  1: magic_shift=37;break;
				case  2: magic_shift=30;break;
//				case  3: magic_shift=35;break;
//				case  4: magic_shift=33;break;
				default: magic_shift=norm_x_e+27;break;
			}	
		printf("magic_shift:%d  ",magic_shift);
		magic_shift=norm_x_e+27;
                for (jCh=idx;jCh<nChannelsC; jCh++)
                {
                    tmp2=W_extract_l(W_shr(singularVectors_Left_64[currChannel][jCh],tmpe));
			printf("s=0x%016llxll;f=0x%08xll;tmp2=%016llx --> 0x%08xll;",singularVectors_Left_64[iCh][jCh],f,singularVectors_Left_64[currChannel][jCh],tmp2);
		    singularVectors_Left_64[iCh][jCh]=W_add(singularVectors_Left_64[iCh][jCh],W_shr(W_mult0_32_32(f,tmp2),magic_shift));
			printf("\x1b[45m%016llX\x1b[0m\x1b[1;35m ",singularVectors_Left_64[iCh][jCh]);
                }
		printf("\x1b[0m\n");
            }
            invVal_e = 0;
            move16();
            tmpe = W_norm(abs_x);
            invVal = BASOP_Util_Divide3232_Scale_newton( abs_x, maxWithSign_fx( r ), &invVal_e);
            invVal_e = add(invVal_e, sub(tmpe, r_e));
		printf("invVal: \x1b[1;36m[%016llx] %08X*%08X = %08X\x1b[0m\n",abs_64,abs_x,r,invVal);
            tmpe = add(1,singularVectors_e);
		printf("SECDIAG: \x1b[1;36m");
            for ( jCh = idx; jCh < nChannelsC ; jCh++)
            {
		printf("[%016llX] ",singularVectors_Left_64[currChannel][jCh]);
                secDiag[jCh] = Mpy_32_32( W_extract_l(W_shr(singularVectors_Left_64[currChannel][jCh],tmpe)), invVal);
                move32();
                secDiag_e[jCh]=add(invVal_e, tmpe);
                move16();
		printf("%08X*%08X= %08X<%2x   ",W_extract_l(W_shr(singularVectors_Left_64[currChannel][jCh],tmpe)), invVal, secDiag[jCh], secDiag_e[jCh]);
            }
		printf("\x1b[0m\n");
        }
    }
    return;
@@ -1424,14 +1439,11 @@ static void biDiagonalReductionLeft_fx(
                f = Mpy_32_32( norm_x, invVal ); /* invVal_e + (norm_x_e - r_e) */
                f_e = add( invVal_e, sub( norm_x_e, r_e ) );

		//printf("MYADD%d,%d: [%08X<%2x]  ",iCh,currChannel,f,f_e);
                FOR( jCh = currChannel; jCh < nChannelsL; jCh++ ) /* nChannelsL */
                {
                    singularVectors[jCh][iCh] = BASOP_Util_Add_Mant32Exp( singularVectors[jCh][iCh], singularVectors2_e[jCh][iCh], Mpy_32_32( f, singularVectors[jCh][currChannel] ), add( f_e, singularVectors2_e[jCh][currChannel] ), &singularVectors2_e[jCh][iCh] );
			//printf("%08X ",singularVectors[jCh][iCh]);
                    move32();
                }
		//printf("\n");
            }
        }
        // rescaling block
@@ -1439,7 +1451,6 @@ static void biDiagonalReductionLeft_fx(
        move32();
        singularValues_e[currChannel] = g_e;
        move16();
	//printf("MYADD done\n");
    }

    return;
@@ -1621,14 +1632,16 @@ static void biDiagonalReductionRight_fx(
        move32();
        abs_x_e = 0;
        move16();
	printf("NORM64: ");
        FOR( jCh = idx; jCh < nChannelsC; jCh++ ) /*nChannelsC */
        {
            norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[currChannel][jCh], singularVectors[currChannel][jCh] ), shl( singularVectors2_e[currChannel][jCh], 1 ), &norm_x_e ); /* exp(norm_x_e) */
            abs_x = BASOP_Util_Add_Mant32Exp( abs_x, abs_x_e, L_abs( singularVectors[currChannel][jCh] ), singularVectors2_e[currChannel][jCh], &abs_x_e );
        }
	printf("%08X<%2X   abs_x:%08X<%2X \x1b[0m\n",norm_x,norm_x_e,abs_x,abs_x_e);

	printf("\n");
	printf("START%d: ",currChannel);
	printf("norm: %08X<%2X   abs: %08X<%2X",norm_x,norm_x_e, abs_x,abs_x_e);
	printf("\n");
        IF( norm_x ) /*(fabsf(*sig_x) > EPSILON * fabsf(*sig_x)) { */
        {
   
@@ -1652,13 +1665,13 @@ static void biDiagonalReductionRight_fx(


            r = BASOP_Util_Add_Mant32Exp( Mpy_32_32( ( *g ), singularVectors[currChannel][idx] ), add( *g_e, singularVectors2_e[currChannel][idx] ), -norm_x, norm_x_e, &r_e );                                      /* exp(r_e) */
		printf("MUL: s:%08X -> g:%08X ->tmpmul:%08X norm:%08X --> r:%08X\n",singularVectors[currChannel][idx],*g,Mpy_32_32( ( *g ), singularVectors[currChannel][idx]),norm_x,r);
		printf("SING: %08X-%08X = ",singularVectors[currChannel][idx],*g);
		printf("G%d: %08X<%2X  R:  %08X<%2X\x1b[0m\n",currChannel,*g,*g_e, r,r_e);
		printf("SING%d: ",currChannel);
		printf("%08X<%02X-%08X<%02X ",singularVectors[currChannel][idx],singularVectors2_e[currChannel][idx], *g,*g_e);
             singularVectors[currChannel][idx] = BASOP_Util_Add_Mant32Exp( singularVectors[currChannel][idx], singularVectors2_e[currChannel][idx], -( *g ), ( *g_e ), &singularVectors2_e[currChannel][idx] ); /* exp(sing_exp) */
		printf("%08X\n",singularVectors[currChannel][idx]);
		printf("%08X<%02X\n",singularVectors[currChannel][idx],singularVectors2_e[currChannel][idx]);
            move32();

		printf("R: --> %08X<%2X g:%08X<%2x\n",r,r_e,*g,*g_e);

            FOR( iCh = currChannel + 1; iCh < nChannelsL; iCh++ ) /*  nChannelsL */
            {
@@ -1666,43 +1679,30 @@ static void biDiagonalReductionRight_fx(
                move32();
                norm_x_e = 0;
                move16();
		printf("norm%d,%d: ",currChannel,iCh);
                FOR( jCh = idx; jCh < nChannelsC; jCh++ ) /* nChannelsC */
                {
                    norm_x = BASOP_Util_Add_Mant32Exp( norm_x, norm_x_e, Mpy_32_32( singularVectors[iCh][jCh], singularVectors[currChannel][jCh] ), add( singularVectors2_e[iCh][jCh], singularVectors2_e[currChannel][jCh] ), &norm_x_e ); /* exp(norm_x_e) */
                }
		printf("%08X<%02X  ",norm_x,norm_x_e);
		printf("\n");
		printf("F: norm %08X<%2X  ",norm_x,norm_x_e);
                norm_x = BASOP_Util_Divide3232_Scale_newton( norm_x, r, &invVal_e );
                norm_x_e = add( invVal_e, sub(norm_x_e, r_e ) );
		printf("-> %08X<%2X\n",norm_x,norm_x_e);
   		printf("singularVectors_Left%d,%d: ",currChannel,iCh); 
                FOR( jCh = idx; jCh < nChannelsC; jCh++ ) /*  nChannelsC */
                {
                    singularVectors[iCh][jCh] = BASOP_Util_Add_Mant32Exp( singularVectors[iCh][jCh], singularVectors2_e[iCh][jCh], Mpy_32_32( norm_x, singularVectors[currChannel][jCh] ), add( norm_x_e, singularVectors2_e[currChannel][jCh] ), &singularVectors2_e[iCh][jCh] ); /* exp(sing_exp2) */
			printf("|%08X*%08X|expected=0x%08Xll;",norm_x, singularVectors[currChannel][jCh],singularVectors[iCh][jCh]);
			printf("%08X<%2X  ",singularVectors[iCh][jCh], singularVectors2_e[iCh][jCh]);
                    move32();
                }
		printf("\n");
            }
            invVal_e = 0;
            move16();
            invVal = BASOP_Util_Divide3232_Scale_newton( abs_x, maxWithSign_fx( r ), &invVal_e );
            invVal_e = add(invVal_e, sub( abs_x_e,r_e ) );
		printf("invVal: %08X*%08X = %08X\n",abs_x,r,invVal);
		printf("SECDIAG: ");
            FOR( jCh = idx; jCh < nChannelsC; jCh++ ) /* nChannelsC */
            {
                secDiag[jCh] = Mpy_32_32( singularVectors[currChannel][jCh], invVal ); /* exp(sing_exp + (singularVectors_e - sig_x_e) */
                move32();
                secDiag_exp[jCh] = add( invVal_e, singularVectors2_e[currChannel][jCh] );
                move16();
		printf("%08X*%08X= %08X<%2x   ",singularVectors[currChannel][jCh], invVal, secDiag[jCh], secDiag_exp[jCh]);
			
            }
		printf("\n");
        }
    }
    return;