From 4f1adee70d88c1df9f024e7da71d093d6829bc71 Mon Sep 17 00:00:00 2001 From: Antonio De Lucreziis Date: Sat, 15 Feb 2025 19:17:32 +0100 Subject: [PATCH] copied old arnoldi code --- CMakeLists.txt | 9 +- README.md | 3 + arnoldi.c | 207 ++++++++++++++++++ matrices/laplacian/julia-laplaciano.jl | 13 ++ .../laplacian/laplacian-discretization-3d.mat | Bin 0 -> 19240 bytes 5 files changed, 229 insertions(+), 3 deletions(-) create mode 100644 arnoldi.c create mode 100644 matrices/laplacian/julia-laplaciano.jl create mode 100644 matrices/laplacian/laplacian-discretization-3d.mat diff --git a/CMakeLists.txt b/CMakeLists.txt index 1d14f4c..6243ed8 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,13 +11,16 @@ set(CMAKE_Fortran_COMPILER gfortran) project(main) -add_executable(main main.c) - find_package(PkgConfig REQUIRED) pkg_search_module(OpenBLAS REQUIRED IMPORTED_TARGET openblas) pkg_search_module(PETSc REQUIRED IMPORTED_TARGET petsc) pkg_search_module(OpenMPI REQUIRED IMPORTED_TARGET ompi) +add_executable(main main.c) target_include_directories(main PUBLIC ${PETSc_INCLUDE_DIRS} ${OpenBLAS_INCLUDE_DIRS} ${OpenMPI_INCLUDE_DIRS}) -target_link_libraries(main PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConfig::OpenMPI) \ No newline at end of file +target_link_libraries(main PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConfig::OpenMPI) + +add_executable(arnoldi arnoldi.c) +target_include_directories(arnoldi PUBLIC ${PETSc_INCLUDE_DIRS} ${OpenBLAS_INCLUDE_DIRS} ${OpenMPI_INCLUDE_DIRS}) +target_link_libraries(arnoldi PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConfig::OpenMPI) diff --git a/README.md b/README.md index 9e1f53c..f3e970a 100644 --- a/README.md +++ b/README.md @@ -4,6 +4,9 @@ ### Spack +> I recently wrote +> [an article about this](https://aziis98.com/articles/using-spack/) + This assumes that you have spack installed and sourced. If not, you can install it using the following commands diff --git a/arnoldi.c b/arnoldi.c new file mode 100644 index 0000000..45728f3 --- /dev/null +++ b/arnoldi.c @@ -0,0 +1,207 @@ +#include +#include +#include +#include +#include + +#include +#include +#include + +static char help[] = "Example PETSc program\n\n"; + +// extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *); +// extern PetscErrorCode ComputeRHS(KSP, Vec, void *); +// extern PetscErrorCode ComputeInitialSolution(DM, Vec); + +PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H); + +int main(int argc, char **argv) { + Vec b; + PetscInt n, l; + + PetscFunctionBeginUser; + PetscInitialize(&argc, &argv, (char *)0, help); + + PetscBool flg; + PetscOptionsGetInt(NULL, NULL, "-n", &n, &flg); + if (!flg) + n = 176; + + PetscOptionsGetInt(NULL, NULL, "-l", &l, &flg); + if (!flg) + l = 4; + + VecCreate(PETSC_COMM_WORLD, &b); + VecSetSizes(b, PETSC_DECIDE, n); + VecSetType(b, VECMPI); + + 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); + 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); + VecView(b, PETSC_VIEWER_DRAW_WORLD); + VecView(b, PETSC_VIEWER_STDOUT_WORLD); + + printf("[Arnoldi] Allocating memory for Krylov subspace basis\n"); + + Vec *Q; + PetscMalloc1(l, &Q); + + for (PetscInt i = 0; i < l; i++) { + VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]); + } + + printf("[Arnoldi] Constructing Hessenberg matrix\n"); + + Mat H; + PetscCall(MatCreate(PETSC_COMM_SELF, &H)); + PetscCall(MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l)); + PetscCall(MatSetType(H, MATDENSE)); + // MatSetType(H, MATMPIAIJ); + + printf("[Arnoldi] Starting iteration\n"); + + PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); + + PetscCall(ArnoldiIteration(A, b, l, Q, H)); + + PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); + + printf("[Arnoldi] Done\n"); + + MatView(H, PETSC_VIEWER_DRAW_WORLD); + MatView(H, PETSC_VIEWER_STDOUT_WORLD); + + // print Hessenberg matrix to file + + PetscViewer v2; + PetscCall(PetscViewerCreate(PETSC_COMM_SELF, &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); + + // for (PetscInt i = 0; i < l + 1; i++) { + // PetscCall(VecDestroy(&Q[i])); + // } + + // PetscCall(PetscFree(Q)); + + PetscCall(MatDestroy(&A)); + PetscCall(VecDestroy(&b)); + PetscFinalize(); + return 0; +} + +PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) { + PetscFunctionBeginUser; + + PetscScalar eps = 1e-12; + PetscInt m; + PetscCall(VecGetSize(b, &m)); + + Vec q; + + PetscCall(MatZeroEntries(H)); + + 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; + 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; + PetscCall(VecDot(Q[j], v, &h)); + PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES)); + PetscCall(VecAXPY(v, -h, Q[j])); + } + + // Normalize + PetscScalar h; + PetscCall(VecNorm(v, NORM_2, &h)); + PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES)); + + // Check for convergence + if (h > eps) { + PetscCall(VecNormalize(v, NULL)); + Q[k] = v; + } else { + break; + } + } + + PetscFunctionReturn(PETSC_SUCCESS); +} 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 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?-)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-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!7fy&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_wPLcy*n@0`!srbk9;EhL{jh0%k_Q|J-tV+|5keL zi}Uj7aQ(Y4dOsPi_sI2r3O&6$FW392^z$~ z>D_s`-p{6|_sI2r4n4g)FW37#dU}ss@9p&T?z~*@=hD-A)4TI>yu`{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^gPhF=Y!rK zWBy6H$CHPWN0AfAN#s2w6+kkt4`@vXR`MJdkW6k0g&HTga*8baDncn>>@8M_xcKAg>^=C2t}Zl6R4d z$p^_L