def glsfit(psr,renormalize=True):
    """Solve local GLS problem using scipy.linalg.cholesky.
    Update psr[...].val and psr[...].err from solution.
    If renormalize=True, normalize each design-matrix column by its norm."""
    
    mask = psr.deleted == 0
    res, err = psr.residuals(removemean=False)[mask], psr.toaerrs[mask]
    M = psr.designmatrix(updatebats=False,fixunits=True,fixsigns=True)[mask,:]
        
    C = numpy.diag((err * 1e-6)**2)
    
    if renormalize:
        norm = numpy.sqrt(numpy.sum(M**2,axis=0))
        M /= norm
    else:
        norm = np.ones_like(M[0,:])
    
    mtcm = numpy.dot(M.T,numpy.dot(numpy.linalg.inv(C),M))
    mtcy = numpy.dot(M.T,numpy.dot(numpy.linalg.inv(C),res))
    
    xvar = numpy.linalg.inv(mtcm)
        
    c = scipy.linalg.cho_factor(mtcm)
    xhat = scipy.linalg.cho_solve(c,mtcy)

    sol = psr.vals()
    psr.vals(sol + xhat[1:] / norm[1:])
    psr.errs(numpy.sqrt(numpy.diag(xvar)[1:]) / norm[1:])
    
    return chisq(psr)
Solve local GLS problem using scipy.linalg.cholesky.
    Update psr[...].val and psr[...].err from solution.
    If renormalize=True, normalize each design-matrix column by its norm.