changed Krylov routine using Lanczos

scikit_tt/solvers/ode.py

@@ -1311,38 +1311,35 @@ def krylov(operator: 'TT', initial_value: 'TT', dimension: int, step_size: float
           Annals of Physics, 411, 167998, 2019
     """
 
-    # construct Krylov subspace basis
-    krylov_tensors = [(1/initial_value.norm())*initial_value]
-    v_tmp = operator@krylov_tensors[0]
-    v_tmp = v_tmp - (v_tmp.transpose(conjugate=True)@krylov_tensors[0])*krylov_tensors[0]
-    v_tmp = v_tmp.ortho(threshold=threshold, max_rank=max_rank)
-    v_tmp = (1/v_tmp.norm())*v_tmp
-    krylov_tensors.append(v_tmp)
-    for i in range(2,dimension):
-        v_tmp = operator@krylov_tensors[-1]
-        v_tmp = v_tmp - (v_tmp.transpose(conjugate=True)@krylov_tensors[-2])*krylov_tensors[-2] - (v_tmp.transpose(conjugate=True)@krylov_tensors[-1])*krylov_tensors[-1]
-        v_tmp = (1/v_tmp.norm())*v_tmp
-        v_tmp = v_tmp.ortho(threshold=threshold, max_rank=max_rank)
-        krylov_tensors.append(v_tmp)
-
-    # compute effective H
+    # construct Krylov subspace basis and effective H
     T = np.zeros([dimension, dimension], dtype=complex)
-    for i in range(dimension):
-        for j in range(dimension):
-            T[i,j] = krylov_tensors[i].transpose(conjugate=True)@operator@krylov_tensors[j]
-
+    krylov_tensors = [initial_value]
+    w_tmp = operator@krylov_tensors[-1]
+    alpha = (w_tmp.transpose(conjugate=True)@krylov_tensors[-1])
+    T[0,0] = alpha
+    w_tmp = w_tmp - alpha*krylov_tensors[-1]
+    w_tmp = w_tmp.ortho(threshold=threshold, max_rank=max_rank)
+    for i in range(1,dimension):
+        beta = w_tmp.norm()
+        T[i,i-1] = beta
+        T[i-1,i] = beta
+        krylov_tensors.append((1/beta)*w_tmp)
+        w_tmp = operator@krylov_tensors[-1]
+        alpha = (w_tmp.transpose(conjugate=True)@krylov_tensors[-1])
+        T[i,i] = alpha
+        w_tmp = w_tmp - alpha*krylov_tensors[-1] - beta*krylov_tensors[-2]
+        w_tmp = w_tmp.ortho(threshold=threshold, max_rank=max_rank)
+
     # compute time-evolved state
     w_tmp = np.zeros([dimension], dtype=complex)
-    for j in range(dimension):
-        w_tmp[j] = krylov_tensors[j].transpose(conjugate=True)@initial_value
+    w_tmp[0] = 1
     w_tmp = expm_multiply(-1j*T*step_size, w_tmp) 
     solution = w_tmp[0]*krylov_tensors[0]
     for j in range(1,dimension):
         solution = solution + w_tmp[j]*krylov_tensors[j]
     solution = solution.ortho(threshold=threshold, max_rank=max_rank)
     if normalize > 0:
         solution = (1 / solution.norm(p=normalize)) * solution
-
     return solution