From 4b6295d0dfb7fdf6bc2de99462e0f7bdb47d896d Mon Sep 17 00:00:00 2001
From: Antonio De Lucreziis <antonio.delucreziis@gmail.com>
Date: Sun, 16 Mar 2025 17:32:43 +0100
Subject: [PATCH] fowfewpfeopwjpfwe

---
 .gitignore                 |  6 +++-
 .vscode/settings.json      |  1 +
 arnoldi.c                  | 30 ++++++++++-------
 matrices/laplacian/eigs.jl | 68 ++++++++++++++++++++++++++++++++++----
 4 files changed, 85 insertions(+), 20 deletions(-)

diff --git a/.gitignore b/.gitignore
index 9fcd6a7..73c33da 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,2 +1,6 @@
 build*/
-.spack-env/
\ No newline at end of file
+.spack-env/
+
+*.log
+*.out
+*.lock?
\ No newline at end of file
diff --git a/.vscode/settings.json b/.vscode/settings.json
index 76ebf34..f8cdad8 100644
--- a/.vscode/settings.json
+++ b/.vscode/settings.json
@@ -4,6 +4,7 @@
         "Arnoldi",
         "AXPY",
         "DHSEQR",
+        "eigvals",
         "Hessenberg",
         "Krylov",
         "LAPACK",
diff --git a/arnoldi.c b/arnoldi.c
index 4392dd9..ef5df8d 100644
--- a/arnoldi.c
+++ b/arnoldi.c
@@ -126,11 +126,16 @@ int main(int argc, char **argv) {
 
     // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Starting iteration\n"));
 
-    // ARNOLDI TIME START
     PetscLogDouble arnoldi_start_time;
     PetscCall(PetscTime(&arnoldi_start_time));
-
-    PetscCall(ArnoldiIteration(A, b, l, n, Q, H));
+    {
+        PetscCall(ArnoldiIteration(A, b, l, n, Q, H));
+    }
+    PetscLogDouble arnoldi_end_time;
+    PetscCall(PetscTime(&arnoldi_end_time));
+    PetscLogDouble arnoldi_time = arnoldi_end_time - arnoldi_start_time;
+    if (rank == 0)
+        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Arnoldi time: %f seconds\n", arnoldi_time));
 
     // PetscCall(MatSetValue(H, 2, 3, -1, INSERT_VALUES));
     // PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
@@ -149,15 +154,16 @@ int main(int argc, char **argv) {
         double *work = (double *)malloc(3 * l * sizeof(double));
         int info;
 
-        // call LAPACK function "DHSEQR" to compute the eigenvalues of the Hessenberg matrix
-        LAPACKE_dhseqr(LAPACK_ROW_MAJOR, 'E', 'I', l, 1, l, H, l, wr, wi, z, l);
-
-        PetscLogDouble arnoldi_end_time;
-        PetscCall(PetscTime(&arnoldi_end_time));
-        PetscLogDouble arnoldi_time = arnoldi_end_time - arnoldi_start_time;
-        if (rank == 0)
-            PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Arnoldi time: %f seconds\n", arnoldi_time));
-        // ARNOLDI TIME END
+        PetscLogDouble lapack_start_time;
+        PetscCall(PetscTime(&lapack_start_time));
+        {
+            // call LAPACK function "DHSEQR" to compute the eigenvalues of the Hessenberg matrix
+            LAPACKE_dhseqr(LAPACK_ROW_MAJOR, 'E', 'I', l, 1, l, H, l, wr, wi, z, l);
+        }
+        PetscLogDouble lapack_end_time;
+        PetscCall(PetscTime(&lapack_end_time));
+        PetscLogDouble lapack_time = lapack_end_time - lapack_start_time;
+        PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] LAPACK time: %f seconds\n", lapack_time));
 
         // // print Hessenberg matrix
         // printf("H = \n");
diff --git a/matrices/laplacian/eigs.jl b/matrices/laplacian/eigs.jl
index d033f12..73cd409 100644
--- a/matrices/laplacian/eigs.jl
+++ b/matrices/laplacian/eigs.jl
@@ -2,16 +2,70 @@ using SparseArrays
 using LinearAlgebra
 using MAT
 
-# 11 x 16
-nx = 10
-ny = 15
+# 20 x 20 x 20 grid
+nx = 20 - 1
+ny = 20 - 1
+nz = 20 - 1
+
 ex = fill(1, nx)
 ey = fill(1, ny)
-Dxx = diagm(-1 => ex, 0 => -2 * ex, +1 => ex)
-Dyy = diagm(-1 => ey, 0 => -2 * ey, +1 => ey)
-L = kron(Dyy, diagm(0 => [ex; 1])) + kron(diagm(0 => [ey; 1]), Dxx);
+ez = fill(1, nz)
+
+Dxx = spdiagm(-1 => ex, 0 => -2 * ex, +1 => ex)
+Dyy = spdiagm(-1 => ey, 0 => -2 * ey, +1 => ey)
+Dzz = spdiagm(-1 => ez, 0 => -2 * ez, +1 => ez)
+
+Ix = spdiagm(0 => [ex; 1])
+Iy = spdiagm(0 => [ey; 1])
+Iz = spdiagm(0 => [ez; 1])
+
+# L = kron(Dxx, Iy, Iz) + kron(Ix, Dyy, Iz) + kron(Ix, Iy, Dzz)
+
+# # 10 x 17 grid
+# nx = 11 - 1
+# ny = 16 - 1
+
+# ex = fill(1, nx)
+# ey = fill(1, ny)
+
+# Dxx = spdiagm(-1 => ex, 0 => -2 * ex, +1 => ex)
+# Dyy = spdiagm(-1 => ey, 0 => -2 * ey, +1 => ey)
+
+# Ix = spdiagm(0 => [ex; 1])
+# Iy = spdiagm(0 => [ey; 1])
+
+# L = kron(Dxx, Iy) + kron(Ix, Dyy)
 
 println("Laplacian matrix L:")
 
 display(sparse(L))
-display(eigvals(L))
\ No newline at end of file
+
+l = 100
+
+# arnoldi iteration
+Q = Matrix{Float64}(undef, size(L, 1), l)
+Q[:, 1] = ones(size(L, 1))
+H = zeros(l, l)
+
+for j in 2:l
+    # Arnoldi iteration
+    v = L * Q[:, j - 1]
+    
+    # Reorthogonalization
+    for i in 1:(j - 1)
+        H[i, j - 1] = dot(Q[:, i], v)
+        v -= H[i, j - 1] * Q[:, i]
+    end
+    H[j, j - 1] = norm(v)
+    if H[j, j - 1] < 1e-10
+        break
+    end
+    Q[:, j] = v / H[j, j - 1]
+end
+
+H = H[1:(l-1), 1:(l-1)]
+
+println("Hessenberg matrix H:")
+display(H)
+
+display(eigvals(H))