From a1937f7eddce07bd6afbfcaa3a2234a443f8f2bb Mon Sep 17 00:00:00 2001
From: Antonio De Lucreziis <antonio.delucreziis@gmail.com>
Date: Sat, 4 May 2024 15:06:12 +0200
Subject: [PATCH] feat: this maybe works

---
 .gitignore                                    |   3 +
 main.c                                        | 160 +++++++++++-------
 matrices/laplacian/julia-laplaciano.jl        |  13 ++
 .../laplacian/laplacian-discretization-3d.mat | Bin 0 -> 19240 bytes
 4 files changed, 116 insertions(+), 60 deletions(-)
 create mode 100644 matrices/laplacian/julia-laplaciano.jl
 create mode 100644 matrices/laplacian/laplacian-discretization-3d.mat

diff --git a/.gitignore b/.gitignore
index 57fd4da..08d45fa 100644
--- a/.gitignore
+++ b/.gitignore
@@ -4,3 +4,6 @@
 
 # Binaries
 build/
+
+# Clang
+.cache/
\ No newline at end of file
diff --git a/main.c b/main.c
index 9fca88c..a17041e 100644
--- a/main.c
+++ b/main.c
@@ -17,7 +17,6 @@ static char help[] = "Example PETSc program\n\n";
 PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H);
 
 int main(int argc, char **argv) {
-    Mat A;
     Vec b;
     PetscInt n, l;
 
@@ -27,7 +26,7 @@ int main(int argc, char **argv) {
     PetscBool flg;
     PetscOptionsGetInt(NULL, NULL, "-n", &n, &flg);
     if (!flg)
-        n = 10;
+        n = 176;
 
     PetscOptionsGetInt(NULL, NULL, "-l", &l, &flg);
     if (!flg)
@@ -40,42 +39,65 @@ int main(int argc, char **argv) {
     VecSet(b, 1.0);
     // VecSetValue(b, 0, 1.0, INSERT_VALUES);
 
+    Mat A;
     MatCreate(PETSC_COMM_WORLD, &A);
-    MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
-
-    MatSetType(A, MATMPIAIJ);
-
-    // A := diag(-1, 2, -1)
-    for (PetscInt i = 0; i < n; i++) {
-        PetscScalar v[3] = {-1.0, 2.0, -1.0};
-        PetscInt col[3] = {i - 1, i, i + 1};
-        PetscInt ncol = 0;
-        if (i > 0) {
-            col[ncol] = i - 1;
-            v[ncol] = -1.0;
-            ncol++;
-        }
-        col[ncol] = i;
-        v[ncol] = 2.0;
-        ncol++;
-        if (i < n - 1) {
-            col[ncol] = i + 1;
-            v[ncol] = -1.0;
-            ncol++;
-        }
-        MatSetValues(A, 1, &i, ncol, col, v, INSERT_VALUES);
-        // MatSetValue(A, i, i, (PetscScalar)(i + 1), INSERT_VALUES);
-    }
+    // MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n);
+    // MatSetType(A, MATMPIAIJ);
+    PetscViewer v;
+    PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &v));
+    PetscCall(PetscViewerSetType(v, PETSCVIEWERHDF5));
+    PetscCall(PetscViewerPushFormat(v, PETSC_VIEWER_HDF5_MAT));
+    PetscCall(PetscViewerSetFromOptions(v));
+    PetscCall(PetscViewerFileSetMode(v, FILE_MODE_READ));
+    PetscCall(PetscViewerFileSetName(
+        v, "../matrices/laplacian/laplacian-discretization-3d.mat"));
+
+    PetscCall(MatSetOptionsPrefix(A, "a_"));
+    PetscCall(PetscObjectSetName((PetscObject)A, "A"));
+
+    // PetscCall(
+    //     PetscOptionsGetString(NULL, NULL, "-f", name, sizeof(name), &flg));
+    // PetscCheck(flg, PETSC_COMM_WORLD, PETSC_ERR_SUP,
+    //            "Must provide a binary file for the matrix");
+
+    // // PetscCall(MatLoad(A, v));
+    // PetscCall(PetscViewerBinaryOpen(
+    //     PETSC_COMM_WORLD,
+    //     "../matrices/laplacian/laplacian-discretization-3d.mat",
+    //     FILE_MODE_READ, &v));
+    PetscCall(MatLoad(A, v));
+
+    // // A := diag(-1, 2, -1)
+    // for (PetscInt i = 0; i < n; i++) {
+    //     PetscScalar v[3] = {-1.0, 2.0, -1.0};
+    //     PetscInt col[3] = {i - 1, i, i + 1};
+    //     PetscInt ncol = 0;
+    //     if (i > 0) {
+    //         col[ncol] = i - 1;
+    //         v[ncol] = -1.0;
+    //         ncol++;
+    //     }
+    //     col[ncol] = i;
+    //     v[ncol] = 2.0;
+    //     ncol++;
+    //     if (i < n - 1) {
+    //         col[ncol] = i + 1;
+    //         v[ncol] = -1.0;
+    //         ncol++;
+    //     }
+    //     MatSetValues(A, 1, &i, ncol, col, v, INSERT_VALUES);
+    //     // MatSetValue(A, i, i, (PetscScalar)(i + 1), INSERT_VALUES);
+    // }
 
     MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
     MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
 
     // MatView(A, PETSC_VIEWER_DRAW_WORLD);
-    MatView(A, PETSC_VIEWER_STDOUT_WORLD);
+    // MatView(A, PETSC_VIEWER_STDOUT_WORLD);
     // VecView(b, PETSC_VIEWER_DRAW_WORLD);
-    VecView(b, PETSC_VIEWER_STDOUT_WORLD);
+    // VecView(b, PETSC_VIEWER_STDOUT_WORLD);
 
-    printf("Allocating memory for Krylov subspace basis\n");
+    printf("[Arnoldi] Allocating memory for Krylov subspace basis\n");
 
     Vec *Q;
     PetscMalloc1(l, &Q);
@@ -84,34 +106,48 @@ int main(int argc, char **argv) {
         VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]);
     }
 
-    printf("Constructing Hessenberg matrix\n");
+    printf("[Arnoldi] Constructing Hessenberg matrix\n");
 
     Mat H;
-    MatCreate(PETSC_COMM_WORLD, &H);
-    MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l);
+    PetscCall(MatCreate(PETSC_COMM_WORLD, &H));
+    PetscCall(MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l));
+    PetscCall(MatSetType(H, MATDENSE));
     // MatSetType(H, MATMPIAIJ);
-    MatSetType(H, MATDENSE);
 
-    printf("Starting Arnoldi iteration\n");
+    printf("[Arnoldi] Starting iteration\n");
+
+    PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
 
-    MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
     PetscCall(ArnoldiIteration(A, b, l, Q, H));
-    MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
 
-    for (PetscInt i = 0; i < l + 1; i++) {
-        VecView(Q[i], PETSC_VIEWER_STDOUT_WORLD);
-    }
+    PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
 
-    MatView(H, PETSC_VIEWER_STDOUT_WORLD);
+    printf("[Arnoldi] Done\n");
 
-    for (PetscInt i = 0; i < l + 1; i++) {
-        VecDestroy(&Q[i]);
-    }
+    // print Hessenberg matrix to file
+
+    PetscViewer v2;
+    PetscCall(PetscViewerCreate(PETSC_COMM_WORLD, &v2));
+    PetscCall(PetscViewerSetType(v2, PETSCVIEWERHDF5));
+    PetscCall(PetscViewerPushFormat(v2, PETSC_VIEWER_HDF5_MAT));
+    PetscCall(PetscViewerFileSetMode(v2, FILE_MODE_WRITE));
+    PetscCall(PetscViewerFileSetName(v2, "hessenberg.mat"));
+    PetscCall(MatView(H, v2));
+
+    // for (PetscInt i = 0; i < l + 1; i++) {
+    //     VecView(Q[i], PETSC_VIEWER_STDOUT_WORLD);
+    // }
+
+    // MatView(H, PETSC_VIEWER_STDOUT_WORLD);
 
-    // PetscFree(Q);
+    // for (PetscInt i = 0; i < l + 1; i++) {
+    //     PetscCall(VecDestroy(&Q[i]));
+    // }
 
-    MatDestroy(&A);
-    VecDestroy(&b);
+    // PetscCall(PetscFree(Q));
+
+    PetscCall(MatDestroy(&A));
+    PetscCall(VecDestroy(&b));
     PetscFinalize();
     return 0;
 }
@@ -121,39 +157,43 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) {
 
     PetscScalar eps = 1e-12;
     PetscInt m;
-    VecGetSize(b, &m);
+    PetscCall(VecGetSize(b, &m));
 
     Vec q;
 
-    MatZeroEntries(H);
+    PetscCall(MatZeroEntries(H));
 
-    VecDuplicate(b, &q);
-    VecCopy(b, q);
-    VecNormalize(q, NULL);
+    PetscCall(VecDuplicate(b, &q));
+    PetscCall(VecCopy(b, q));
+    PetscCall(VecNormalize(q, NULL));
 
     Q[0] = q;
 
     for (PetscInt k = 1; k < n + 1; k++) {
+        // printf("[Arnoldi] Iteration %d\n", k);
+
         Vec v;
-        VecDuplicate(b, &v);
-        MatMult(A, Q[k - 1], v);
+        PetscCall(VecDuplicate(b, &v));
+        PetscCall(MatMult(A, Q[k - 1], v));
 
         // Reorthogonalization using modified Gram-Schmidt
         for (PetscInt j = 0; j < k; j++) {
+            // printf("[Arnoldi] Reorthogonalization %d\n", j);
+
             PetscScalar h;
-            VecDot(Q[j], v, &h);
-            MatSetValue(H, j, k - 1, h, INSERT_VALUES);
-            VecAXPY(v, -h, Q[j]);
+            PetscCall(VecDot(Q[j], v, &h));
+            PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES));
+            PetscCall(VecAXPY(v, -h, Q[j]));
         }
 
         // Normalize
         PetscScalar h;
-        VecNorm(v, NORM_2, &h);
-        MatSetValue(H, k, k - 1, h, INSERT_VALUES);
+        PetscCall(VecNorm(v, NORM_2, &h));
+        PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES));
 
         // Check for convergence
         if (h > eps) {
-            VecNormalize(v, NULL);
+            PetscCall(VecNormalize(v, NULL));
             Q[k] = v;
         } else {
             break;
diff --git a/matrices/laplacian/julia-laplaciano.jl b/matrices/laplacian/julia-laplaciano.jl
new file mode 100644
index 0000000..dcc29bd
--- /dev/null
+++ b/matrices/laplacian/julia-laplaciano.jl
@@ -0,0 +1,13 @@
+using SparseArrays
+using MAT
+
+# 11 x 16
+nx = 10
+ny = 15
+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)
+L = kron(Dyy, spdiagm(0 => [ex; 1])) + kron(spdiagm(0 => [ey; 1]), Dxx);
+
+matwrite("laplacian-discretization-3d.mat", Dict("A" => L))
diff --git a/matrices/laplacian/laplacian-discretization-3d.mat b/matrices/laplacian/laplacian-discretization-3d.mat
new file mode 100644
index 0000000000000000000000000000000000000000..bad09f42c3a122e61ba35e339fac6d86e60ba805
GIT binary patch
literal 19240
zcmeI4X>eV07>CcjDTy{E6s1K&(I{GkHYk#kDs3rZU)o4pEh*ZxYAKDqweNeWt#+*$
zjAaHhn5hqZ=m*R&gBgrv1~Ztkw(h;>c{1@%OA!0XhyIfP^LtNjp8q*V>K?{5G>vT-
zSv`D6ZS}Z@rtPOpYn`0@Wgq$d#{FgAEBi*PG_YXwKBIc~?ANmdD=Ly6Nmbg*r?HZc
znlV5BI{hk1Ch|#+oKN1*e6DcBBHM`dDvtl|8tKAl)|(D)Y%G52_cOs+{Lc-HZfqD|
zRB!|yD#&G=XM9;Y*NUH)$}awa%q6{6dH&X{nI7~ij-8#qe?ap8V03lo?|klqD)2tq
zUFUDU<Id|H)l9!*&!nQ`nklZV$^AK3mCH3$m-(^^t5v1XmzAeW-nT209ZCyjN$<ju
z3dbt^7E$<xVnR!6b6Z=|r!Z9f0sX0_ojhyU&ZT~3QpeodNoA+Gwo{vDwoNX5v-A?L
z?-*WM8sLvy;gyz-%As@ooVi@Ouq!X!EGrzcNU91)^}xkLtx^2mh0N#krNIi0Q<f_>
ziyzB2PBZ;;J*JE0?`wulNs{)?_;36Cx9+q1`Oy~XOHOw&Nq*q;D9I?Gm>Ew6MK6B6
zVm1DL&71RCFFhaiVRpBdY@brk-R-rzkC*LJwogR^f7X4VyU*8u`th>sDZ8Gs>uKdm
z?<cFvu3pF4^^{#t+4YoNPhGuV{8xI}eJ{K3W%s@8zV~0bpLO-;$nJaDeJ{K3{rUI3
z)m>-)eJ$g^{656*et)cRa(C_hu{XC#?*#vj_s4v{t7gmeRr`Kdr}sPa(`LT8c}k1_
zzVVmeUwn+B+425jVZCeq7VVuL;^*%G|5M*z?BuD(q~=-8Px1Uir}g}aEgKEgUzVAH
z%s^%!Gw}Bs@b{%&;bpk)<8zhpxb71`z2Wv<@4@=oWL)?0rgh-<UFRPoJYAGNEmW+!
z&f`YOy3EIQ9@0tHgWGp~E!O*xaXo$THA&Wo+jo5f)(t)CYeRtEBG+4Adirx-u0NlX
z)nnv(^r0kuIxjCgcvbr7N9sFrz4oW4U+3le-GrW=BiHk$^z`k#T;H3~(|hFlA3)E2
zab8{x*T4Is_s!vYk6iCO4pn%XE&Xs_uJ?g(y+^Khe%(y;?z~*@TX8(SN3Qp+>FM2h
zx!woS(|hE4uc4=R=jD3mm;H|GB=TyGGnjN=oR@D4*T4Is_wC?%k6iED)6={2a=j0s
zr}xP9K9ru`otNv~zc=VTa=q`!@$~MzT<<&4(|hE4ucfDV=jD1&pRP<2y+^M9T{vDf
z>AZXx>-u+J^gbM}_sI2LM^Eq0%k{o1J-tV+_uc5}-Fdm*cc-WK$n`#gp5C38>wOP;
zdXHT1d(zXp^K!lKMNjXM>%X3!`{KO30j_`dMeifwdXHT1d(+dq^K!lKLr?FK>%Ebl
z-kq20eH1;tN3QpM>FM2hx!y<9(|hE4-;bW&otM|c_b2rpx&Ft{b6=d7kA>^sebM_k
zxZWez`vLUy?z~*@2h!7f<a!@ZPw&pl^?ndNy+^M1gX!trdAZ&Xp{MuA_1;8J@6OBh
zeki?qGIISN26tbammdz-zx$&1Bj9?ET<=HH)4TI>y&pwS?~&{MXnJ~gUat3J=;=Lj
zy&p?Y@6OBhejGi$N3Qn?^z`n$T<;U<={<7&H`8-poR=RD*T4Is_ZGO`BiH*RdU|(W
zuJ_6G^d7n1r_j^8^K!kPKu_<H>wPLcy*n@0`!srbk9;EhL{jh0%k_Q|J-tV+|5keL
zi}Uj7aQ(Y4dOsPi_sI2r3O&6$FW392^z<IN-e=I$yYq6r&!ngK$o1YvPw&pl^*)Q9
z-XquhY<d$(=jD1o4X*dd^*@K6`{KO(bh!TA7rmbW*L&o8Ka-x`otNwVEP8s6T<>$~
z>D_s`-p{6|_sI2r4n4g)FW37#dU}ss@9p&T?z~*@=hD-A<oZ94p8Mjw{Cv3n-50%I
z0M~oudcTmK-kq20{UUmLk6iB;)6={2a=p)|r}xP9zJQ+IotL-6FCq0Fx!y0Or+4S&
zdcTaG-Xqul<@DSa=jB(x_3ysu{YtprBiH*?^z`n$T<=%Y(|hE4zlNUPotNwVT6%hq
zT<_P>)4TI>y<bnSos3-XH^BAoyj<@$($jn7`oD>u`{KO(X1M;{7royC*L&o8zm=Ze
zotNwVHhOxGT<;6%>D_s`-fyR;_sI2r2R*$zFW37bdU}ss?|0JEyYq6r-$hUFk?a3%
zdhUz!@_XR=cVF~=FI?}D>wPgjy*n@0`+fBE9=YD{r>A%4<$8aBp57y01b>j!yYq6r
zKSWROk?Z|odU|(WuJ=di={<7&KT6MiabCUzu7CGM@A{7G`qs03*Y&L5xUOHl+IL;A
z`i$%PT*T|Q@46oK7uWTtH~X&ZO<!?cZ+f)P%Sp=kq@3(Uu0!@A`;z_10pvh(5IL9}
zLe`RX<R0Wma$j;Rc@TLRc{DkJY$2zT)5)3S9C9vsF6sUAQoi<j_PX)>^gPhF=Y!rK
zWBy6<S@H$)W%3pBRdPA`2KhGm9{C~pDftEYHTfO6f=vIb@+2wGr3;nhy5xprf3lj~
ziX2Q1C2Psy<Os5X97T>H$CHPWN0AfAN#s<rm7GDIMxI5slNXYgkXMq|kvEeI$ver#
z<U{0A@(J=8@_F(l@)hzmayj`r`4;&u`2qO}`8oL&`7ODEEX$|+Q%?3GtH=$=e&hgh
zOL7}>2w6+kkt4`@vXR`MJdkW6k0g&HTga*8baDncn>>@8M_xcKAg>^=C2t}Zl6R4d
z$p^_L<m2Si<a6YUr2jvbSNVE5`3CtG`40I$`7!w!`6c-c`6HR<Mp#Z(l2v3MvM<@6
V96$~v2a(&7L&;ilIO+dC;wP~XXzc(1

literal 0
HcmV?d00001