function x = pd_solve( L, b);
% solve positive def. system LL^tx=b
% input: L ... lower triangle ( cholesky factor)
%        b ... right hand side
% out:   x solution
% call: x = pd_solve( L, b);

[n, n1] = size( L);
% first solve Ly = b;
b( 1) = b( 1) / L( 1, 1);
for i = 2:n
        b( i) = (b( i) - L( i, 1:i-1) * b( 1:i-1) ) / L( i, i);
        end;

% now solve L^t x = b = Ly
b = b';
b( n) = b( n) / L( n,n);
for i = n-1 : - 1 : 1
        b( i) = (b( i) - b( i+1: n) * L( i+1:n, i) ) / L( i,i);
        end;
x = b';
