/* dposl.c --- solve symmetric positive definite matrix.
 *	       call dpofa or dpoco prior to calling dposl.
 */

void dposl(double **a,int n,double *b)
{
	double t;
	int i,k,kb;

/*	
 *	solve trans(r)*y = b
 */
 	for (k=0;k<n;k++) {
/*		t = ddot(k-1,a(1,k),1,b(1),1)		*/
		t = 0;
		for (i = 0;i < k;i++)
			t += (a[i][k] * b[i]);
		b[k] = (b[k]-t)/a[k][k];
	}

/*	
 *	solve r*x = y
 */
	for (kb=0;kb<n;kb++) {
		k = n - 1 - kb;
		b[k] /= a[k][k];
		t = -b[k];
/*		daxpy(k-1,t,a(1,k),1,b(1),1)		*/
		for (i=0;i<k;i++)
			b[i] += (t * a[i][k]);
	}
}