From bdf39a41f7bb38efa24d16ace7859351a707669a Mon Sep 17 00:00:00 2001
From: alberto <albedef@protonmail.com>
Date: Sun, 12 Apr 2026 17:34:37 +0200
Subject: [PATCH] Try demmel's implementation.

---
 src/afgl/util/lanczos.py | 25 ++++++++++++-------------
 1 file changed, 12 insertions(+), 13 deletions(-)

diff --git a/src/afgl/util/lanczos.py b/src/afgl/util/lanczos.py
index 82d5cea..30548a9 100644
--- a/src/afgl/util/lanczos.py
+++ b/src/afgl/util/lanczos.py
@@ -3,6 +3,7 @@ import numpy.linalg as LA
 
 """
 Classic Lanczos method (without re-orthogonalization)
+Using Demmel's book version.
 
 Arguments
 L : Real valued NxN symmetric matrix
@@ -23,22 +24,20 @@ beta : ndarray
 def lanczos(L, s, M):
     N = len(s)
     alp = np.zeros(M)
-    beta = np.zeros(M - 1)
-    V = np.zeros((N, M))
-    V[:, 0] = s / LA.norm(s)
+    beta = np.zeros(M)
+    V = np.zeros((N, M + 1))
+    V[:, 1] = s / LA.norm(s)
 
-    for j in range(M):
+    for j in range(1, M):
         w = L @ V[:, j]
         alp[j] = np.dot(V[:, j], w)
 
-        v_tilde = w - V[:, j] * alp[j]
-        if j > 0:
-            v_tilde = v_tilde - V[:, j - 1] * beta[j - 1]
+        w = w - V[:, j] * alp[j] - V[:, j - 1] * beta[j - 1]
 
-        if j < M - 1:
-            beta[j] = LA.norm(v_tilde)
-            if beta[j] == 0:
-                break
-            V[:, j + 1] = v_tilde / beta[j]
+        beta[j] = LA.norm(w)
+        if beta[j] == 0:
+            print("Breakdown")
+            break
+        V[:, j + 1] = w / beta[j]
 
-    return [V, alp, beta]
+    return [V[:, 1:], alp, beta[1:]]