code: def glsfit(psr,renormalize=True):
      
      
      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)
docs: 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.