diff --git a/.clang-format b/.clang-format index 321c6de..36e5e57 100644 --- a/.clang-format +++ b/.clang-format @@ -1,2 +1,3 @@ # .clang-format IndentWidth: 4 +ColumnLimit: 120 diff --git a/CMakeLists.txt b/CMakeLists.txt index 6243ed8..9216cc9 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -24,3 +24,7 @@ target_link_libraries(main PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConf 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) + +add_executable(ex10 ex10.c) +target_include_directories(ex10 PUBLIC ${PETSc_INCLUDE_DIRS} ${OpenBLAS_INCLUDE_DIRS} ${OpenMPI_INCLUDE_DIRS}) +target_link_libraries(ex10 PUBLIC m PkgConfig::PETSc PkgConfig::OpenBLAS PkgConfig::OpenMPI) diff --git a/arnoldi.c b/arnoldi.c index 45728f3..3e26a0b 100644 --- a/arnoldi.c +++ b/arnoldi.c @@ -1,20 +1,41 @@ +#include + +#include +#include +#include #include #include #include #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); +void swap(double *a, double *b) { + double t = *a; + *a = *b; + *b = t; +} + +// int compare_double(const void *a, const void *b) { +// int aa = *(int *)a; +// int bb = *(int *)b; -PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H); +// if (aa < bb) +// return -1; +// if (aa > bb) +// return 1; +// return 0; +// } + +PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, double *h); int main(int argc, char **argv) { Vec b; @@ -26,78 +47,47 @@ int main(int argc, char **argv) { PetscBool flg; PetscOptionsGetInt(NULL, NULL, "-n", &n, &flg); if (!flg) - n = 176; + n = -1; 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); + l = 10; 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(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); + // TIME: Assembly construction - printf("[Arnoldi] Allocating memory for Krylov subspace basis\n"); + if (n == -1) { + MatGetSize(A, &n, NULL); + } + + VecCreate(PETSC_COMM_WORLD, &b); + VecSetSizes(b, PETSC_DECIDE, n); + VecSetType(b, VECMPI); + + VecSet(b, 1.0); + // VecSetValue(b, 0, 1.0, INSERT_VALUES); + + // MatView(A, PETSC_VIEWER_STDOUT_WORLD); + // VecView(b, PETSC_VIEWER_STDOUT_WORLD); + + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Allocating memory for Krylov subspace basis\n")); Vec *Q; PetscMalloc1(l, &Q); @@ -106,65 +96,110 @@ int main(int argc, char **argv) { VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]); } - printf("[Arnoldi] Constructing Hessenberg matrix\n"); + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Constructing Hessenberg matrix\n")); + + // Mat H; + // PetscCall(MatCreate(PETSC_COMM_SELF, &H)); + // PetscCall(PetscObjectSetName((PetscObject)H, "hessenberg")); + // PetscCall(MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l)); + // PetscCall(MatSetType(H, MATDENSE)); + + // PetscCall(MatCreate(PETSC_COMM_WORLD, &H)); - 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"); + double *H = (double *)malloc((l + 1) * l * sizeof(double)); - PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Starting iteration\n")); - PetscCall(ArnoldiIteration(A, b, l, Q, H)); + // TIME:START - PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); + PetscCall(ArnoldiIteration(A, b, l, n, Q, H)); - printf("[Arnoldi] Done\n"); + // PetscCall(MatSetValue(H, 2, 3, -1, INSERT_VALUES)); + // PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); + // PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); - MatView(H, PETSC_VIEWER_DRAW_WORLD); - MatView(H, PETSC_VIEWER_STDOUT_WORLD); + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Done\n")); - // print Hessenberg matrix to file + int rank; + PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); - 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)); + if (rank == 0) { + // print Hessenberg matrix + printf("H = \n"); + for (int i = 0; i < l + 1; i++) { + for (int j = 0; j < l; j++) { + printf("%.3f ", H[i * l + j]); + } + printf("\n"); + } - // for (PetscInt i = 0; i < l + 1; i++) { - // VecView(Q[i], PETSC_VIEWER_STDOUT_WORLD); - // } + // call LAPACK function "DHSEQR" to compute the eigenvalues of the Hessenberg matrix + + double *wr = (double *)malloc(l * sizeof(double)); + double *wi = (double *)malloc(l * sizeof(double)); + double *z = (double *)malloc(l * l * sizeof(double)); + double *work = (double *)malloc(3 * l * sizeof(double)); + int info; + + LAPACKE_dhseqr(LAPACK_ROW_MAJOR, 'E', 'I', l, 1, l, H, l, wr, wi, z, l); + + // TIME:END + + // sort eigenvalues + for (int i = 0; i < l; i++) { + for (int j = i + 1; j < l; j++) { + if (wr[i] > wr[j]) { + swap(&wr[i], &wr[j]); + swap(&wi[i], &wi[j]); + for (int k = 0; k < l; k++) { + swap(&z[i * l + k], &z[j * l + k]); + } + } + } + } - // MatView(H, PETSC_VIEWER_STDOUT_WORLD); + // print eigenvalues + printf("Eigenvalues = \n"); + for (int i = 0; i < l; i++) { + printf("%.3f + %.3f i\n", wr[i], wi[i]); + } + } - // for (PetscInt i = 0; i < l + 1; i++) { - // PetscCall(VecDestroy(&Q[i])); + // print eigenvectors + // printf("Eigenvectors = \n"); + // for (int i = 0; i < l; i++) { + // for (int j = 0; j < l; j++) { + // printf("%f ", z[i * l + j]); + // } + // printf("\n"); // } - // PetscCall(PetscFree(Q)); - PetscCall(MatDestroy(&A)); PetscCall(VecDestroy(&b)); + // PetscCall(MatDestroy(&H)); + PetscFinalize(); return 0; } -PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) { +PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, double *h) { PetscFunctionBeginUser; + int rank; + MPI_Comm_rank(PETSC_COMM_WORLD, &rank); + printf("Rank %d\n", rank); + PetscScalar eps = 1e-12; - PetscInt m; - PetscCall(VecGetSize(b, &m)); Vec q; - PetscCall(MatZeroEntries(H)); + for (PetscInt i = 0; i < n + 1; i++) { + for (PetscInt j = 0; j < n; j++) { + h[i * n + j] = 0.0; + } + } PetscCall(VecDuplicate(b, &q)); PetscCall(VecCopy(b, q)); @@ -173,32 +208,39 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) { Q[0] = q; for (PetscInt k = 1; k < n + 1; k++) { - // printf("[Arnoldi] Iteration %d\n", k); + + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[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); + for (PetscInt j = 0; j < k; j++) { // anche solo 3 + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Reorthogonalization %d %d\n", k, j)); + + PetscScalar h_ij; + PetscCall(VecDot(Q[j], v, &h_ij)); - PetscScalar h; - PetscCall(VecDot(Q[j], v, &h)); - PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES)); - PetscCall(VecAXPY(v, -h, Q[j])); + // PetscCall(MatSetValue(H, j, k - 1, h, INSERT_VALUES)); + h[j * n + k - 1] = h_ij; + + PetscCall(VecAXPY(v, -h_ij, Q[j])); } // Normalize - PetscScalar h; - PetscCall(VecNorm(v, NORM_2, &h)); - PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES)); + PetscScalar h_ij; + PetscCall(VecNorm(v, NORM_2, &h_ij)); + + // PetscCall(MatSetValue(H, k, k - 1, h, INSERT_VALUES)); + h[k * n + k - 1] = h_ij; // Check for convergence - if (h > eps) { + if (h_ij > eps) { PetscCall(VecNormalize(v, NULL)); Q[k] = v; } else { + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Early breakdown")); break; } } diff --git a/ex10.c b/ex10.c new file mode 100644 index 0000000..0d0d9f2 --- /dev/null +++ b/ex10.c @@ -0,0 +1,82 @@ +#include +#include +#include +#include + +static char help[] = "Test to write HDF5 file from PETSc DMDA Vec.\n\n"; + +int main(int argc, char **argv) { + DM da2D; + PetscInt i, j, ixs, ixm, iys, iym; + PetscViewer H5viewer; + PetscScalar xm = -1.0, xp = 1.0; + PetscScalar ym = -1.0, yp = 1.0; + PetscScalar value = 1.0, dx, dy; + PetscInt Nx = 40, Ny = 40; + Vec gauss, input; + PetscScalar **gauss_ptr; + PetscReal norm; + const char *vecname; + + dx = (xp - xm) / (Nx - 1); + dy = (yp - ym) / (Ny - 1); + + /* Initialize the Petsc context */ + PetscFunctionBeginUser; + PetscCall(PetscInitialize(&argc, &argv, NULL, help)); + PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, Nx, Ny, + PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da2D)); + PetscCall(DMSetFromOptions(da2D)); + PetscCall(DMSetUp(da2D)); + + /* Set the coordinates */ + PetscCall(DMDASetUniformCoordinates(da2D, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0)); + + /* Declare gauss as a DMDA component */ + PetscCall(DMCreateGlobalVector(da2D, &gauss)); + PetscCall(PetscObjectSetName((PetscObject)gauss, "pressure")); + + /* Initialize vector gauss with a constant value (=1) */ + PetscCall(VecSet(gauss, value)); + + /* Get the coordinates of the corners for each process */ + PetscCall(DMDAGetCorners(da2D, &ixs, &iys, 0, &ixm, &iym, 0)); + + /* Build the gaussian profile (exp(-x^2-y^2)) */ + PetscCall(DMDAVecGetArray(da2D, gauss, &gauss_ptr)); + for (j = iys; j < iys + iym; j++) { + for (i = ixs; i < ixs + ixm; i++) + gauss_ptr[j][i] = PetscExpScalar(-(xm + i * dx) * (xm + i * dx) - (ym + j * dy) * (ym + j * dy)); + } + PetscCall(DMDAVecRestoreArray(da2D, gauss, &gauss_ptr)); + + /* Create the HDF5 viewer */ + PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_WRITE, &H5viewer)); + PetscCall(PetscViewerSetFromOptions(H5viewer)); + + /* Write the H5 file */ + PetscCall(VecView(gauss, H5viewer)); + + /* Close the viewer */ + PetscCall(PetscViewerDestroy(&H5viewer)); + + PetscCall(VecDuplicate(gauss, &input)); + PetscCall(PetscObjectGetName((PetscObject)gauss, &vecname)); + PetscCall(PetscObjectSetName((PetscObject)input, vecname)); + + /* Create the HDF5 viewer for reading */ + PetscCall(PetscViewerHDF5Open(PETSC_COMM_WORLD, "gauss.h5", FILE_MODE_READ, &H5viewer)); + PetscCall(PetscViewerSetFromOptions(H5viewer)); + PetscCall(VecLoad(input, H5viewer)); + PetscCall(PetscViewerDestroy(&H5viewer)); + + PetscCall(VecAXPY(input, -1.0, gauss)); + PetscCall(VecNorm(input, NORM_2, &norm)); + PetscCheck(norm <= 1.e-6, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "Vec read in does not match vector written out"); + + PetscCall(VecDestroy(&input)); + PetscCall(VecDestroy(&gauss)); + PetscCall(DMDestroy(&da2D)); + PetscCall(PetscFinalize()); + return 0; +} diff --git a/main.c b/main.c index 2e81cc9..ad2f629 100644 --- a/main.c +++ b/main.c @@ -107,6 +107,7 @@ int main(int argc, char **argv) { PetscCall( PetscPrintf(PETSC_COMM_WORLD, "Vector length %" PetscInt_FMT "\n", n)); + PetscCall(VecMax(x, &maxind, &maxval)); PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMax %g, VecInd %" PetscInt_FMT "\n", diff --git a/matrices/laplacian/eigs.jl b/matrices/laplacian/eigs.jl new file mode 100644 index 0000000..d033f12 --- /dev/null +++ b/matrices/laplacian/eigs.jl @@ -0,0 +1,17 @@ +using SparseArrays +using LinearAlgebra +using MAT + +# 11 x 16 +nx = 10 +ny = 15 +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); + +println("Laplacian matrix L:") + +display(sparse(L)) +display(eigvals(L)) \ No newline at end of file