From efdae14faeb47829c77e1726765107f0bd3315fa Mon Sep 17 00:00:00 2001 From: Antonio De Lucreziis Date: Thu, 13 Mar 2025 16:57:15 +0100 Subject: [PATCH] fewpjfiewjipfwe --- .vscode/settings.json | 2 + arnoldi.c | 66 ++++++---- ex10.c | 82 ------------- experiment-1.sh | 12 ++ main.c | 276 ------------------------------------------ 5 files changed, 54 insertions(+), 384 deletions(-) delete mode 100644 ex10.c create mode 100644 experiment-1.sh delete mode 100644 main.c diff --git a/.vscode/settings.json b/.vscode/settings.json index 00f9e85..7e29fcf 100644 --- a/.vscode/settings.json +++ b/.vscode/settings.json @@ -11,9 +11,11 @@ "MATDENSE", "MATMPIAIJ", "matwrite", + "mpirun", "PETSC", "PETSCVIEWERHDF", "Reorthogonalization", + "SBATCH", "spdiagm", "VECMPI" ] diff --git a/arnoldi.c b/arnoldi.c index ca91c20..f8f0968 100644 --- a/arnoldi.c +++ b/arnoldi.c @@ -8,6 +8,7 @@ #include #include #include +#include #include #include #include @@ -59,6 +60,9 @@ int main(int argc, char **argv) { if (!flg) snprintf(matrix_name, sizeof(matrix_name), "./matrices/laplacian/laplacian-discretization-3d.mat"); + PetscLogDouble load_start_time; + PetscCall(PetscTime(&load_start_time)); + Mat A; MatCreate(PETSC_COMM_WORLD, &A); PetscViewer v; @@ -77,7 +81,11 @@ int main(int argc, char **argv) { MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); - // TIME: Assembly construction + PetscLogDouble load_start_end; + PetscCall(PetscTime(&load_start_end)); + + PetscLogDouble load_time = load_start_end - load_start_time; + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Load time: %f seconds\n", load_time)); if (n == -1) { MatGetSize(A, &n, NULL); @@ -94,7 +102,7 @@ int main(int argc, char **argv) { // 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")); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Allocating memory for Krylov subspace basis\n")); Vec *Q; PetscMalloc1(l, &Q); @@ -103,7 +111,7 @@ int main(int argc, char **argv) { VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]); } - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Constructing Hessenberg matrix\n")); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Constructing Hessenberg matrix\n")); // Mat H; // PetscCall(MatCreate(PETSC_COMM_SELF, &H)); @@ -117,9 +125,11 @@ int main(int argc, char **argv) { double *H = (double *)malloc((l + 1) * l * sizeof(double)); - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Starting iteration\n")); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Starting iteration\n")); - // TIME:START + // ARNOLDI TIME START + PetscLogDouble arnoldi_start_time; + PetscCall(PetscTime(&arnoldi_start_time)); PetscCall(ArnoldiIteration(A, b, l, n, Q, H)); @@ -127,22 +137,12 @@ int main(int argc, char **argv) { // PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY)); // PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY)); - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Done\n")); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Done\n")); int rank; PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank)); 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"); - } - - // 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)); @@ -150,9 +150,23 @@ 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); - // TIME:END + PetscLogDouble arnoldi_end_time; + PetscCall(PetscTime(&arnoldi_end_time)); + PetscLogDouble arnoldi_time = arnoldi_end_time - arnoldi_start_time; + PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Arnoldi time: %f seconds\n", arnoldi_time)); + // ARNOLDI TIME END + + // 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"); + } // sort eigenvalues for (int i = 0; i < l; i++) { @@ -194,13 +208,13 @@ int main(int argc, char **argv) { PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, double *h) { PetscFunctionBeginUser; - int rank, total; - MPI_Comm_rank(PETSC_COMM_WORLD, &rank); - MPI_Comm_size(PETSC_COMM_WORLD, &total); + // int rank, total; + // MPI_Comm_rank(PETSC_COMM_WORLD, &rank); + // MPI_Comm_size(PETSC_COMM_WORLD, &total); - char hostname[64]; - PetscCall(PetscGetHostName(hostname, sizeof(hostname))); - printf("Process %s: %d of %d\n", hostname, rank, total); + // char hostname[64]; + // PetscCall(PetscGetHostName(hostname, sizeof(hostname))); + // printf("Process %s: %d of %d\n", hostname, rank, total); PetscScalar eps = 1e-12; @@ -220,7 +234,7 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, do for (PetscInt k = 1; k < n + 1; k++) { - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Iteration %d\n", k)); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Iteration %d\n", k)); Vec v; PetscCall(VecDuplicate(b, &v)); @@ -228,7 +242,7 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, do // Reorthogonalization using modified Gram-Schmidt for (PetscInt j = 0; j < k; j++) { // anche solo 3 - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Reorthogonalization %d %d\n", k, j)); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Reorthogonalization %d %d\n", k, j)); PetscScalar h_ij; PetscCall(VecDot(Q[j], v, &h_ij)); @@ -251,7 +265,7 @@ PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, PetscInt m, Vec *Q, do PetscCall(VecNormalize(v, NULL)); Q[k] = v; } else { - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Early breakdown")); + // PetscCall(PetscPrintf(PETSC_COMM_WORLD, "[Arnoldi] Early breakdown")); break; } } diff --git a/ex10.c b/ex10.c deleted file mode 100644 index 0d0d9f2..0000000 --- a/ex10.c +++ /dev/null @@ -1,82 +0,0 @@ -#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/experiment-1.sh b/experiment-1.sh new file mode 100644 index 0000000..0937b17 --- /dev/null +++ b/experiment-1.sh @@ -0,0 +1,12 @@ +#!/bin/bash +#SBATCH --job-name=experiment-1 +#SBATCH --nodes=20 +#SBATCH --output=%x_%j.log + +for i in {1..30} +do + echo "Iteration $i" + mpirun -np $i ./build/arnoldi -l 100 + sleep 1 +done + diff --git a/main.c b/main.c deleted file mode 100644 index ad2f629..0000000 --- a/main.c +++ /dev/null @@ -1,276 +0,0 @@ -static char help[] = "Basic vector routines.\n\n"; - -/* - Include "petscvec.h" so that we can use vectors. Note that this file - automatically includes: - petscsys.h - base PETSc routines petscis.h - index sets - petscviewer.h - viewers -*/ - -#include -#include -#include - -#include - -int example_lapack() { - int n = 2; // Dimension of the matrix A - int nrhs = 1; // Number of right-hand sides - int lda = 2; // Leading dimension of A - int ldb = 2; // Leading dimension of B - int info; // Status output - int ipiv[2]; // Pivot indices for the LU decomposition - - // Define matrix A (row-major order) - double A[4] = {3, 1, 1, 2}; - - // Define right-hand side vector b - double b[2] = {9, 8}; - - // Call LAPACK dgesv to solve Ax = b - info = LAPACKE_dgesv(LAPACK_ROW_MAJOR, n, nrhs, A, lda, ipiv, b, ldb); - - // Check for success - if (info == 0) { - printf("Solution:\n"); - for (int i = 0; i < n; i++) { - printf("x[%d] = %f\n", i, b[i]); - } - } else if (info < 0) { - printf("Argument %d had an illegal value.\n", -info); - } else { - printf("Matrix is singular. Solution could not be computed.\n"); - } - - return info; -} - -int main(int argc, char **argv) { - printf("LAPACKE Example:\n"); - example_lapack(); - - printf("PETSC Example:\n"); - Vec x, y, w; /* vectors */ - Vec *z; /* array of vectors */ - PetscReal norm, v, v1, v2, maxval; - PetscInt n = 20, maxind; - PetscScalar one = 1.0, two = 2.0, three = 3.0, dots[3], dot; - - PetscFunctionBeginUser; - PetscCall(PetscInitialize(&argc, &argv, NULL, help)); - PetscCall(PetscOptionsGetInt(NULL, NULL, "-n", &n, NULL)); - - /* - Create a vector, specifying only its global dimension. - When using VecCreate(), VecSetSizes() and VecSetFromOptions(), the vector - format (currently parallel, shared, or sequential) is determined at - runtime. Also, the parallel partitioning of the vector is determined by - PETSc at runtime. - */ - PetscCall(VecCreate(PETSC_COMM_WORLD, &x)); - PetscCall(VecSetSizes(x, PETSC_DECIDE, n)); - PetscCall(VecSetFromOptions(x)); - - /* - Duplicate some work vectors (of the same format and - partitioning as the initial vector). - */ - PetscCall(VecDuplicate(x, &y)); - PetscCall(VecDuplicate(x, &w)); - - /* - Duplicate more work vectors (of the same format and - partitioning as the initial vector). Here we duplicate - an array of vectors, which is often more convenient than - duplicating individual ones. - */ - PetscCall(VecDuplicateVecs(x, 3, &z)); - /* - Set the vectors to entries to a constant value. - */ - PetscCall(VecSet(x, one)); - PetscCall(VecSet(y, two)); - PetscCall(VecSet(z[0], one)); - PetscCall(VecSet(z[1], two)); - PetscCall(VecSet(z[2], three)); - /* - Demonstrate various basic vector routines. - */ - PetscCall(VecDot(x, y, &dot)); - PetscCall(VecMDot(x, 3, z, dots)); - - /* - Note: If using a complex numbers version of PETSc, then - PETSC_USE_COMPLEX is defined in the makefiles; otherwise, - (when using real numbers) it is undefined. - */ - - 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", - (double)maxval, maxind)); - - PetscCall(VecMin(x, &maxind, &maxval)); - PetscCall(PetscPrintf(PETSC_COMM_WORLD, - "VecMin %g, VecInd %" PetscInt_FMT "\n", - (double)maxval, maxind)); - PetscCall(PetscPrintf(PETSC_COMM_WORLD, - "All other values should be near zero\n")); - - PetscCall(VecScale(x, two)); - PetscCall(VecNorm(x, NORM_2, &norm)); - v = norm - 2.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecScale %g\n", (double)v)); - - PetscCall(VecCopy(x, w)); - PetscCall(VecNorm(w, NORM_2, &norm)); - v = norm - 2.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecCopy %g\n", (double)v)); - - PetscCall(VecAXPY(y, three, x)); - PetscCall(VecNorm(y, NORM_2, &norm)); - v = norm - 8.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAXPY %g\n", (double)v)); - - PetscCall(VecAYPX(y, two, x)); - PetscCall(VecNorm(y, NORM_2, &norm)); - v = norm - 18.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecAYPX %g\n", (double)v)); - - PetscCall(VecSwap(x, y)); - PetscCall(VecNorm(y, NORM_2, &norm)); - v = norm - 2.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecSwap %g\n", (double)v)); - PetscCall(VecNorm(x, NORM_2, &norm)); - v = norm - 18.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecSwap %g\n", (double)v)); - - PetscCall(VecWAXPY(w, two, x, y)); - PetscCall(VecNorm(w, NORM_2, &norm)); - v = norm - 38.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecWAXPY %g\n", (double)v)); - - PetscCall(VecPointwiseMult(w, y, x)); - PetscCall(VecNorm(w, NORM_2, &norm)); - v = norm - 36.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall( - PetscPrintf(PETSC_COMM_WORLD, "VecPointwiseMult %g\n", (double)v)); - - PetscCall(VecPointwiseDivide(w, x, y)); - PetscCall(VecNorm(w, NORM_2, &norm)); - v = norm - 9.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall( - PetscPrintf(PETSC_COMM_WORLD, "VecPointwiseDivide %g\n", (double)v)); - - PetscCall(VecSetValue(y, 0, 0.0, INSERT_VALUES)); - PetscCall(VecAssemblyBegin(y)); - PetscCall(VecAssemblyEnd(y)); - PetscCall(VecPointwiseDivide(w, x, y)); - PetscCall(VecNorm(w, NORM_2, &norm)); - v = norm - 9.0 * PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall( - PetscPrintf(PETSC_COMM_WORLD, "VecPointwiseDivide %g\n", (double)v)); - - dots[0] = one; - dots[1] = three; - dots[2] = two; - - PetscCall(VecSet(x, one)); - PetscCall(VecMAXPY(x, 3, dots, z)); - PetscCall(VecNorm(z[0], NORM_2, &norm)); - v = norm - PetscSqrtReal((PetscReal)n); - if (v > -PETSC_SMALL && v < PETSC_SMALL) - v = 0.0; - PetscCall(VecNorm(z[1], NORM_2, &norm)); - v1 = norm - 2.0 * PetscSqrtReal((PetscReal)n); - if (v1 > -PETSC_SMALL && v1 < PETSC_SMALL) - v1 = 0.0; - PetscCall(VecNorm(z[2], NORM_2, &norm)); - v2 = norm - 3.0 * PetscSqrtReal((PetscReal)n); - if (v2 > -PETSC_SMALL && v2 < PETSC_SMALL) - v2 = 0.0; - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "VecMAXPY %g %g %g \n", (double)v, - (double)v1, (double)v2)); - - /* - Free work space. All PETSc objects should be destroyed when they - are no longer needed. - */ - PetscCall(VecDestroy(&x)); - PetscCall(VecDestroy(&y)); - PetscCall(VecDestroy(&w)); - PetscCall(VecDestroyVecs(3, &z)); - PetscCall(PetscFinalize()); - return 0; -} - -/*TEST - - testset: - output_file: output/ex1_1.out - # This is a test where the exact numbers are critical - diff_args: -j - - test: - - test: - suffix: cuda - args: -vec_type cuda - requires: cuda - - test: - suffix: kokkos - args: -vec_type kokkos - requires: kokkos_kernels - - test: - suffix: hip - args: -vec_type hip - requires: hip - - test: - suffix: 2 - nsize: 2 - - test: - suffix: 2_cuda - nsize: 2 - args: -vec_type cuda - requires: cuda - - test: - suffix: 2_kokkos - nsize: 2 - args: -vec_type kokkos - requires: kokkos_kernels - - test: - suffix: 2_hip - nsize: 2 - args: -vec_type hip - requires: hip - -TEST*/