solchol.c

speech signal process tools

solchol(R,b,x,t,n)
/* R upper triangular from cholesky, b constants, x unknowns */
/* t is a temporary vector */
/* Mx+b=0   M=Rt*R */
register double R[],b[],*x,t[];
register int n;
{
	register int k,np1,i,iM;
	register double *M,dtemp;
	register double S;
	double dot();

	np1=n+1;

/* forward substitution */

	for(k=0,M=R;k<n;k++,M+=np1)
		{
		dtemp= -(b[k]+dot(R+k,n,t,1,k))/(*M);
		t[k]=dtemp;
		}
	/*
	for(k=0,kM=0;k<n;k++,kM+=n)
		{
		S=b[k];
		for(i=0,iM=0;i<k;i++,iM+=n)
			S+=R[iM+k]*t[i];
		dtemp= -S/R[kM+k];
		t[k]=dtemp;
		}
	*/

/* backward substitution */

	/*
	for(k=0,x=x+n,M=R+n*n;k<n;k++,M-=np1,x--)
		{
		dtemp=(t[n-1-k]-dot(M,1,x,1,k))/(*(M-1));
		x[k]=dtemp;
		}
	*/
	/*
	for(i=n-1,iM=(n-1)*n;i>=0;i--,iM-=n)
		{
		S=t[i];
		for(k=i+1;k<n;k++)
			S-=R[iM+k]*x[k];
		dtemp=S/R[iM+i];
		x[i]=dtemp;
		}
	*/
	for(i=n-1,iM=(n-1)*n;i>=0;i--,iM-=n)
		{
		S=t[i]-dot(R+iM+i+1,1,x+i+1,1,n-i-1);
		dtemp=S/R[iM+i];
		x[i]=dtemp;
		}
	return;
}