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 0000000..bad09f4 Binary files /dev/null and b/matrices/laplacian/laplacian-discretization-3d.mat differ