From 01bc6b3297c62cab785ceb4faef8dd753e80a378 Mon Sep 17 00:00:00 2001
From: singularitti <singularitti@outlook.com>
Date: Tue, 12 May 2026 21:01:49 -0500
Subject: [PATCH] Vectorize Rayleigh-Ritz projection fallback

Add rayleighRitzProduct to replace repeated MATLAB nested loops for
building the Rayleigh-Ritz projection matrix with a BLAS-backed matrix
multiply. Preserve the old upper-triangle-then-mirror behavior so the
MATLAB fallback remains consistent with the existing MEX implementation.
---
 RSDFT/DiagonalizationFiles/chefsi1.m             | 7 +------
 RSDFT/DiagonalizationFiles/chsubsp.m             | 9 +--------
 RSDFT/DiagonalizationFiles/first_filt.m          | 9 +--------
 RSDFT/DiagonalizationFiles/lanczos.m             | 7 +------
 RSDFT/DiagonalizationFiles/lanczosForChsubsp.m   | 7 +------
 RSDFT/DiagonalizationFiles/rayleighRitzProduct.m | 6 ++++++
 6 files changed, 11 insertions(+), 34 deletions(-)
 create mode 100644 RSDFT/DiagonalizationFiles/rayleighRitzProduct.m

diff --git a/RSDFT/DiagonalizationFiles/chefsi1.m b/RSDFT/DiagonalizationFiles/chefsi1.m
index fbd8cce..0f8c6c0 100644
--- a/RSDFT/DiagonalizationFiles/chefsi1.m
+++ b/RSDFT/DiagonalizationFiles/chefsi1.m
@@ -47,12 +47,7 @@
 if (OPTIMIZATIONLEVEL~=0)
 	G=Rayleighritz(Vin,W,n2);
 else
-    for j=1:n2
-        for i=1:j
-            G(i,j) = Vin(:,i)'*W(:,j);
-            G(j,i) = G(i,j);
-        end
-    end
+    G = rayleighRitzProduct(Vin,W);
     if (enableMexFilesTest==1)
         G2=Rayleighritz(Vin,W,n2);
         if (any(abs(G-G2)>0.000001))
diff --git a/RSDFT/DiagonalizationFiles/chsubsp.m b/RSDFT/DiagonalizationFiles/chsubsp.m
index eb862b3..bbd46ce 100644
--- a/RSDFT/DiagonalizationFiles/chsubsp.m
+++ b/RSDFT/DiagonalizationFiles/chsubsp.m
@@ -47,14 +47,7 @@
     if (OPTIMIZATIONLEVEL~=0)
         G=Rayleighritz(Vin,W,n2);
     else
-        %preallocating G
-        G=zeros(n2,n2);
-        for j=1:n2
-            for i=1:j
-                G(i,j) = Vin(:,i)'*W(:,j);
-                G(j,i) = G(i,j);
-            end
-        end
+        G = rayleighRitzProduct(Vin,W);
         if (enableMexFilesTest==1)
             G2=Rayleighritz(Vin,W,n2);
             if (any(abs(G-G2)>0.000001))
diff --git a/RSDFT/DiagonalizationFiles/first_filt.m b/RSDFT/DiagonalizationFiles/first_filt.m
index 309bc2b..26f78d4 100644
--- a/RSDFT/DiagonalizationFiles/first_filt.m
+++ b/RSDFT/DiagonalizationFiles/first_filt.m
@@ -52,14 +52,7 @@
     if (OPTIMIZATIONLEVEL~=0)
         G=Rayleighritz(Vin,W,n2);
     else
-        % reuse G and overwrite it
-        G=zeros(n2,n2);
-        for j=1:n2
-            for i=1:j
-                G(i,j) = Vin(:,i)'*W(:,j);
-                G(j,i) = G(i,j);
-            end
-        end
+        G = rayleighRitzProduct(Vin,W);
     end
 
     [Q, D] = eig(G);
diff --git a/RSDFT/DiagonalizationFiles/lanczos.m b/RSDFT/DiagonalizationFiles/lanczos.m
index 5375958..34d3ada 100644
--- a/RSDFT/DiagonalizationFiles/lanczos.m
+++ b/RSDFT/DiagonalizationFiles/lanczos.m
@@ -197,12 +197,7 @@
     if (OPTIMIZATIONLEVEL~=0)
         G=Rayleighritz(Vin,W,mev);
     else
-        for j=1:mev
-            for i=1:j
-                G(i,j) = Vin(:,i)'*W(:,j);
-                G(j,i) = G(i,j);
-            end
-        end
+        G = rayleighRitzProduct(Vin,W);
         if (enableMexFilesTest==1)
             G2=Rayleighritz(Vin,W,mev);
             if (any(abs(G-G2)>0.000001))
diff --git a/RSDFT/DiagonalizationFiles/lanczosForChsubsp.m b/RSDFT/DiagonalizationFiles/lanczosForChsubsp.m
index 9d2a1d3..19eb447 100644
--- a/RSDFT/DiagonalizationFiles/lanczosForChsubsp.m
+++ b/RSDFT/DiagonalizationFiles/lanczosForChsubsp.m
@@ -172,12 +172,7 @@
     if (OPTIMIZATIONLEVEL~=0)
         G=Rayleighritz(Vin,W,mev);
     else
-        for j=1:mev
-            for i=1:j
-                G(i,j) = Vin(:,i)'*W(:,j);
-                G(j,i) = G(i,j);
-            end
-        end
+        G = rayleighRitzProduct(Vin,W);
         if (enableMexFilesTest==1)
             G2=Rayleighritz(Vin,W,mev);
             if (any(abs(G-G2)>0.000001))
diff --git a/RSDFT/DiagonalizationFiles/rayleighRitzProduct.m b/RSDFT/DiagonalizationFiles/rayleighRitzProduct.m
new file mode 100644
index 0000000..4a59822
--- /dev/null
+++ b/RSDFT/DiagonalizationFiles/rayleighRitzProduct.m
@@ -0,0 +1,6 @@
+function G = rayleighRitzProduct(Vin,W)
+%% G = rayleighRitzProduct(Vin,W)
+%% computes the same Rayleigh-Ritz projection as the old MATLAB loops.
+
+G = Vin' * W;
+G = triu(G) + triu(G,1)';