cholesky update for RA

This commit is contained in:
James Hensman 2013-05-03 14:00:22 +01:00
parent 7561c4c232
commit d9252d0e36
3 changed files with 35 additions and 2 deletions

View file

@ -276,3 +276,30 @@ def symmetrify_murray(A):
nn = A.shape[0]
A[[range(nn),range(nn)]] /= 2.0
def cholupdate(L,x):
"""
update the LOWER cholesky factor of a pd matrix IN PLACE
if L is the lower chol. of K, then this function computes L_
where L_ is the lower chol of K + x*x^T
"""
support_code = """
#include <math.h>
"""
code="""
double r,c,s;
int j,i;
for(j=0; j<N; j++){
r = sqrt(L(j,j)*L(j,j) + x(j)*x(j));
c = r / L(j,j);
s = x(j) / L(j,j);
L(j,j) = r;
for (i=j+1; i<N; i++){
L(i,j) = (L(i,j) + s*x(i))/c;
x(i) = c*x(i) - s*L(i,j);
}
}
"""
x = x.copy()
N = x.size
weave.inline(code, support_code=support_code, arg_names=['N','L','x'], type_converters=weave.converters.blitz)