diff --git a/CMakeLists.txt b/CMakeLists.txt index d59429d..e65a87a 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,29 +1,57 @@ +# cmake_minimum_required(VERSION 3.5.0) + +# set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) + +# # set root of location to find PETSc's pkg-config +# set(PETSC $ENV{PETSC_DIR}/$ENV{PETSC_ARCH}) +# set(ENV{PKG_CONFIG_PATH} ${PETSC}/lib/pkgconfig) + +# # Remove the lines below if you do not wish to have PETSc determine the compilers +# execute_process ( COMMAND pkg-config PETSc --variable=ccompiler COMMAND tr -d '\n' OUTPUT_VARIABLE C_COMPILER) +# SET(CMAKE_C_COMPILER ${C_COMPILER}) +# execute_process ( COMMAND pkg-config PETSc --variable=cxxcompiler COMMAND tr -d '\n' OUTPUT_VARIABLE CXX_COMPILER) +# if (CXX_COMPILER) +# SET(CMAKE_CXX_COMPILER ${CXX_COMPILER}) +# endif (CXX_COMPILER) +# execute_process ( COMMAND pkg-config PETSc --variable=fcompiler COMMAND tr -d '\n' OUTPUT_VARIABLE FORTRAN_COMPILER) +# if (FORTRAN_COMPILER) +# SET(CMAKE_Fortran_COMPILER ${FORTRAN_COMPILER}) +# enable_language(Fortran) +# endif (FORTRAN_COMPILER) + +# # tells CMake to build the application ex1 from the source file ex1.c +# # this must appear AFTER the compilers are set +# project(main) +# add_executable(main main.c) + +# find_package(MPI REQUIRED) +# include_directories(${MPI_INCLUDE_PATH}) + +# find_package(PkgConfig REQUIRED) +# pkg_search_module(PETSC REQUIRED IMPORTED_TARGET PETSc) +# target_link_libraries(main PkgConfig::PETSC) + cmake_minimum_required(VERSION 3.5.0) set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) -# set root of location to find PETSc's pkg-config +# Set root of location to find PETSc's pkg-config set(PETSC $ENV{PETSC_DIR}/$ENV{PETSC_ARCH}) set(ENV{PKG_CONFIG_PATH} ${PETSC}/lib/pkgconfig) # Remove the lines below if you do not wish to have PETSc determine the compilers -execute_process ( COMMAND pkg-config PETSc --variable=ccompiler COMMAND tr -d '\n' OUTPUT_VARIABLE C_COMPILER) -SET(CMAKE_C_COMPILER ${C_COMPILER}) -execute_process ( COMMAND pkg-config PETSc --variable=cxxcompiler COMMAND tr -d '\n' OUTPUT_VARIABLE CXX_COMPILER) -if (CXX_COMPILER) - SET(CMAKE_CXX_COMPILER ${CXX_COMPILER}) -endif (CXX_COMPILER) -execute_process ( COMMAND pkg-config PETSc --variable=fcompiler COMMAND tr -d '\n' OUTPUT_VARIABLE FORTRAN_COMPILER) -if (FORTRAN_COMPILER) - SET(CMAKE_Fortran_COMPILER ${FORTRAN_COMPILER}) - enable_language(Fortran) -endif (FORTRAN_COMPILER) - -# tells CMake to build the application ex1 from the source file ex1.c -# this must appear AFTER the compilers are set +SET(CMAKE_C_COMPILER gcc) +SET(CMAKE_CXX_COMPILER g++) +SET(CMAKE_Fortran_COMPILER gfortran) +enable_language(Fortran) + +# Tells CMake to build the application ex1 from the source file ex1.c +# This must appear AFTER the compilers are set project(main) add_executable(main main.c) +find_package(MPI REQUIRED) find_package(PkgConfig REQUIRED) + pkg_search_module(PETSC REQUIRED IMPORTED_TARGET PETSc) -target_link_libraries(main PkgConfig::PETSC) +target_link_libraries(main MPI::MPI_C MPI::MPI_CXX MPI::MPI_Fortran PkgConfig::PETSC) diff --git a/main.c b/main.c index dfc1152..3112af0 100644 --- a/main.c +++ b/main.c @@ -1,173 +1,161 @@ -/* -Laplacian in 3D. Modeled by the partial differential equation +#include +#include +#include +#include +#include - - Laplacian u = 1,0 < x,y,z < 1, +#include +#include +#include -with boundary conditions +static char help[] = "Example PETSc program\n\n"; - u = 1 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1. +// extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *); +// extern PetscErrorCode ComputeRHS(KSP, Vec, void *); +// extern PetscErrorCode ComputeInitialSolution(DM, Vec); - This uses multigrid to solve the linear system +PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H); - See src/snes/tutorials/ex50.c +int main(int argc, char **argv) { + Mat A; + Vec b; + PetscInt n, l; - Can also be run with -pc_type exotic -ksp_type fgmres + PetscFunctionBeginUser; + PetscInitialize(&argc, &argv, (char *)0, help); + + PetscBool flg; + PetscOptionsGetInt(NULL, NULL, "-n", &n, &flg); + if (!flg) + n = 10; + + 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); + + 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); + } -*/ + MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY); + MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY); -static char help[] = "Solves 3D Laplacian using multigrid.\n\n"; + // MatView(A, PETSC_VIEWER_DRAW_WORLD); + MatView(A, PETSC_VIEWER_STDOUT_WORLD); + // VecView(b, PETSC_VIEWER_DRAW_WORLD); + VecView(b, PETSC_VIEWER_STDOUT_WORLD); -#include -#include -#include + printf("Allocating memory for Krylov subspace basis\n"); -extern PetscErrorCode ComputeMatrix(KSP, Mat, Mat, void *); -extern PetscErrorCode ComputeRHS(KSP, Vec, void *); -extern PetscErrorCode ComputeInitialGuess(KSP, Vec, void *); + Vec *Q; + PetscMalloc1(l, &Q); -int main(int argc, char **argv) { - KSP ksp; - PetscReal norm; - DM da; - Vec x, b, r; - Mat A; + for (PetscInt i = 0; i < l; i++) { + VecCreateMPI(PETSC_COMM_WORLD, PETSC_DECIDE, n, &Q[i]); + } - PetscFunctionBeginUser; - PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); - - PetscCall(KSPCreate(PETSC_COMM_WORLD, &ksp)); - PetscCall(DMDACreate3d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, - DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 7, 7, 7, - PETSC_DECIDE, PETSC_DECIDE, PETSC_DECIDE, 1, 1, 0, 0, - 0, &da)); - PetscCall(DMSetFromOptions(da)); - PetscCall(DMSetUp(da)); - PetscCall(KSPSetDM(ksp, da)); - PetscCall(KSPSetComputeInitialGuess(ksp, ComputeInitialGuess, NULL)); - PetscCall(KSPSetComputeRHS(ksp, ComputeRHS, NULL)); - PetscCall(KSPSetComputeOperators(ksp, ComputeMatrix, NULL)); - PetscCall(DMDestroy(&da)); - - PetscCall(KSPSetFromOptions(ksp)); - PetscCall(KSPSolve(ksp, NULL, NULL)); - PetscCall(KSPGetSolution(ksp, &x)); - PetscCall(KSPGetRhs(ksp, &b)); - PetscCall(VecDuplicate(b, &r)); - PetscCall(KSPGetOperators(ksp, &A, NULL)); - - PetscCall(MatMult(A, x, r)); - PetscCall(VecAXPY(r, -1.0, b)); - PetscCall(VecNorm(r, NORM_2, &norm)); - PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Residual norm %g\n", (double)norm)); - - PetscCall(VecDestroy(&r)); - PetscCall(KSPDestroy(&ksp)); - PetscCall(PetscFinalize()); - return 0; -} + printf("Constructing Hessenberg matrix\n"); -PetscErrorCode ComputeRHS(KSP ksp, Vec b, void *ctx) { - PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs; - DM dm; - PetscScalar Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy; - PetscScalar ***barray; + Mat H; + MatCreate(PETSC_COMM_WORLD, &H); + MatSetSizes(H, PETSC_DECIDE, PETSC_DECIDE, l + 1, l); + // MatSetType(H, MATMPIAIJ); + MatSetType(H, MATDENSE); - PetscFunctionBeginUser; - PetscCall(KSPGetDM(ksp, &dm)); - PetscCall(DMDAGetInfo(dm, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0)); - Hx = 1.0 / (PetscReal)(mx - 1); - Hy = 1.0 / (PetscReal)(my - 1); - Hz = 1.0 / (PetscReal)(mz - 1); - HxHydHz = Hx * Hy / Hz; - HxHzdHy = Hx * Hz / Hy; - HyHzdHx = Hy * Hz / Hx; - PetscCall(DMDAGetCorners(dm, &xs, &ys, &zs, &xm, &ym, &zm)); - PetscCall(DMDAVecGetArray(dm, b, &barray)); - - for (k = zs; k < zs + zm; k++) { - for (j = ys; j < ys + ym; j++) { - for (i = xs; i < xs + xm; i++) { - if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || - k == mz - 1) { - barray[k][j][i] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx); - } else { - barray[k][j][i] = Hx * Hy * Hz; - } - } - } + printf("Starting Arnoldi iteration\n"); + + 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(DMDAVecRestoreArray(dm, b, &barray)); - PetscFunctionReturn(PETSC_SUCCESS); + + MatView(H, PETSC_VIEWER_STDOUT_WORLD); + + for (PetscInt i = 0; i < l + 1; i++) { + VecDestroy(&Q[i]); + } + + // PetscFree(Q); + + MatDestroy(&A); + VecDestroy(&b); + PetscFinalize(); + return 0; } -PetscErrorCode ComputeInitialGuess(KSP ksp, Vec b, void *ctx) { +PetscErrorCode ArnoldiIteration(Mat A, Vec b, PetscInt n, Vec *Q, Mat H) { PetscFunctionBeginUser; - PetscCall(VecSet(b, 0)); - PetscFunctionReturn(PETSC_SUCCESS); -} -PetscErrorCode ComputeMatrix(KSP ksp, Mat jac, Mat B, void *ctx) { - DM da; - PetscInt i, j, k, mx, my, mz, xm, ym, zm, xs, ys, zs; - PetscScalar v[7], Hx, Hy, Hz, HxHydHz, HyHzdHx, HxHzdHy; - MatStencil row, col[7]; + PetscScalar eps = 1e-12; + PetscInt m; + VecGetSize(b, &m); - PetscFunctionBeginUser; - PetscCall(KSPGetDM(ksp, &da)); - PetscCall(DMDAGetInfo(da, 0, &mx, &my, &mz, 0, 0, 0, 0, 0, 0, 0, 0, 0)); - Hx = 1.0 / (PetscReal)(mx - 1); - Hy = 1.0 / (PetscReal)(my - 1); - Hz = 1.0 / (PetscReal)(mz - 1); - HxHydHz = Hx * Hy / Hz; - HxHzdHy = Hx * Hz / Hy; - HyHzdHx = Hy * Hz / Hx; - PetscCall(DMDAGetCorners(da, &xs, &ys, &zs, &xm, &ym, &zm)); - - for (k = zs; k < zs + zm; k++) { - for (j = ys; j < ys + ym; j++) { - for (i = xs; i < xs + xm; i++) { - row.i = i; - row.j = j; - row.k = k; - if (i == 0 || j == 0 || k == 0 || i == mx - 1 || j == my - 1 || - k == mz - 1) { - v[0] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx); - PetscCall(MatSetValuesStencil(B, 1, &row, 1, &row, v, INSERT_VALUES)); - } else { - v[0] = -HxHydHz; - col[0].i = i; - col[0].j = j; - col[0].k = k - 1; - v[1] = -HxHzdHy; - col[1].i = i; - col[1].j = j - 1; - col[1].k = k; - v[2] = -HyHzdHx; - col[2].i = i - 1; - col[2].j = j; - col[2].k = k; - v[3] = 2.0 * (HxHydHz + HxHzdHy + HyHzdHx); - col[3].i = row.i; - col[3].j = row.j; - col[3].k = row.k; - v[4] = -HyHzdHx; - col[4].i = i + 1; - col[4].j = j; - col[4].k = k; - v[5] = -HxHzdHy; - col[5].i = i; - col[5].j = j + 1; - col[5].k = k; - v[6] = -HxHydHz; - col[6].i = i; - col[6].j = j; - col[6].k = k + 1; - PetscCall(MatSetValuesStencil(B, 1, &row, 7, col, v, INSERT_VALUES)); - } - } + Vec q; + + MatZeroEntries(H); + + VecDuplicate(b, &q); + VecCopy(b, q); + VecNormalize(q, NULL); + + Q[0] = q; + + for (PetscInt k = 1; k < n + 1; k++) { + Vec v; + VecDuplicate(b, &v); + MatMult(A, Q[k - 1], v); + + for (PetscInt j = 0; j < k; j++) { + PetscScalar h; + VecDot(Q[j], v, &h); + MatSetValue(H, j, k - 1, h, INSERT_VALUES); + VecAXPY(v, -h, Q[j]); + } + + PetscScalar h; + VecNorm(v, NORM_2, &h); + MatSetValue(H, k, k - 1, h, INSERT_VALUES); + + if (h > eps) { + VecNormalize(v, NULL); + Q[k] = v; + } else { + break; } } - PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); - PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); + PetscFunctionReturn(PETSC_SUCCESS); }