The problem is template-mesh deformation to match with silhouettes.  In this 
process, there are two kinds of linear systems to solve.  This system        
(Yoshiyasu/image_interp) is a smooth vector field construction from images,  
which is harmonic interpolation (minimizing laplacian: Lx=0) of intensity    
gradient field p.  This can be solved by normal equation and cholesky        
factorization, x=(A1'*A1)/(A1'*b1), where A1=[L;C] and                       
b1=[zeros(size(length(L),1);1);C*p]. C is a square diagonal matrix containing
weights.  This is for a 400x300 image, so Ix=reshape(x,400,300) must be done 
to get the vector field. After solving y direction for Iy, the result is     
visualized with quiver(Ix,Iy).   At each iteration the both C submatrix and  
the right-hand-side change but L remains unchanged.  [Note by T. Davis:      
since C is of high rank, update/downdate will not be effective, since it is  
meant for low-rank changes.]