reverse.c

cholesky decomposition for dsp c6

#include <stdio.h>

#define N 5

void cholesky(short input[][N], short output[], short d[], int n);

//short A[N][N]= {1,1,1,1,1,1,2,3,4,5,1,3,6,10,15,1,4,10,20,35,1,5,15,35,70}; 
//short A[N][N]= {1,2,3,4,6,2,13,18,23,39,3,18,106,104,117,4,23,104,121,129,6,39,117,129,216};
//short A[N][N]= {2,2,2,2,2,2,5,8,11,14,2,8,21,41,68,2,11,41,104,212,2,14,68,212,514};
short A[N][N]= {22,22,22,22,22,22,35,48,61,74,22,48,91,151,228,22,61,151,322,604,22,74,228,604,1345};
short B[N][N] = {0};
short index[N-1] = {0,1,3,6};
short d[N];
short L[N*(N-1)/2] = {0};

int main(void)
{ 
	//cholesky decomposed
	int i,j;
	for (i = 0; i < N; i++)
	{
		for (j = 0; j <N; j++)
		{
			B[i][j] = A[i][j] << 4;
		}
	}
	//B is Q(16,12);
	cholesky(B, L, d, N);
	return 1;
}


void cholesky(short input[][N], short output[], short d[], int n)//input, output and d are all Q(16,12)
{

	int i,j,k;
	int temp1, temp2;
	short t[N*(N-1)/2];
	d[0] = input[0][0];		//Q(16,12)
	for (i = 1; i < N; i++)
	{
		temp1 = 0;
		for (j = 0; j < i; j++)
		{
			temp2 = 0;
			for (k = 0; k < j; k++)
			{
				temp2 += (t[index[i-1]+k] * output[index[j-1]+k] + 0x0008)>>4;	//rounding, temp2 is Q(32,28)
			} 
			t[index[i-1]+j] = input[i][j] - (short)temp2;//Q(16,12)
			output[index[i-1]+j] = (t[index[i-1]+j] * ((1 << 15)/ d[j]) + 0x0400)>> 11; //rounding, output is Q(16,12)		
			temp1 += (t[index[i-1]+j] * output[index[i-1]+j] + 0x0008) >> 4; //temp2 is Q(32,28)			
		}	
		d[i] = input[i][i] - (short)temp1; //Q(16,12)
	}	
}