def symorth(S):
    "Symmetric orthogonalization"
    E,U = np.linalg.eigh(S)
    n = len(E)
    Shalf = np.identity(n,'d')
    for i in range(n):
        Shalf[i,i] /= np.sqrt(E[i])
    return simx(Shalf,U,True)