From 30b1861a461ce17eea82cb4bdfb64ab18e4c372b Mon Sep 17 00:00:00 2001 From: kadir Date: Wed, 20 Jan 2016 09:00:17 +0200 Subject: [PATCH] smaples for krylov methods are added --- README.md | 9 + krylov/CMakeLists.txt | 80 ++++++++ krylov/README.md | 60 ++++++ krylov/ksp_solver_multi_rhs.c | 221 ++++++++++++++++++++++ krylov/ksp_solver_simple.c | 259 ++++++++++++++++++++++++++ krylov/ksp_solver_two_sys.c | 335 ++++++++++++++++++++++++++++++++++ krylov/run.sh | 11 ++ 7 files changed, 975 insertions(+) create mode 100644 krylov/CMakeLists.txt create mode 100644 krylov/README.md create mode 100644 krylov/ksp_solver_multi_rhs.c create mode 100644 krylov/ksp_solver_simple.c create mode 100644 krylov/ksp_solver_two_sys.c create mode 100755 krylov/run.sh diff --git a/README.md b/README.md index 9e3f567..0605bb3 100644 --- a/README.md +++ b/README.md @@ -12,3 +12,12 @@ In this Sparse Linear Algebra folder, there are sample codes for beginners, as w - spgemm: Multiplication of two sparse matrices. - mkl_shmem: Using MKL's routine mkl_dcsrmultcsr() on a multicore processor - mkl_xphi: Using MKL's routine mkl_dcsrmultcsr() via offloading to a Xeon Phi coprocessor +- Krylov Subspace Methods + - Linear system solution in parallel + - 2D Laplacian (2D mesh) + - Repeatedly solving two linear systems + - Same preconditioner + - Two different matrices with the same nonzero pattern + - Solving multiple linear systems + - Same cofficient matrix + - Different right-hand-side vectors diff --git a/krylov/CMakeLists.txt b/krylov/CMakeLists.txt new file mode 100644 index 0000000..5e93f55 --- /dev/null +++ b/krylov/CMakeLists.txt @@ -0,0 +1,80 @@ + +# ================================================================================================== +# This file is part of the CodeVault project. The project is licensed under Apache Version 2.0. +# CodeVault is part of the EU-project PRACE-4IP (WP7.3.C). +# +# Author(s): +# Valeriu Codreanu +# +# ================================================================================================== +cmake_minimum_required(VERSION 2.8.10 FATAL_ERROR) +list(APPEND CMAKE_MODULE_PATH "${CMAKE_CURRENT_SOURCE_DIR}/../../../../cmake/Modules") + +set(CMAKE_VERBOSE_MAKEFILE ON) + + +## OGUZ: Edit accordingly +include(${CMAKE_CURRENT_SOURCE_DIR}/../../../cmake/common.cmake) + +# ================================================================================================== + +if ("${DWARF_PREFIX}" STREQUAL "") + set(DWARF_PREFIX 2_sparse) +endif() + +## OGUZ: Add executables. There must be more convenient way :-) +set(NAME ${DWARF_PREFIX}_krylov_simple) +set(NAMEB ${DWARF_PREFIX}_krylov_twosys) +set(NAMEC ${DWARF_PREFIX}_krylov_multirhs) + +# ================================================================================================== +# C++ compiler settings + +find_package(Common) +find_package(PETSc REQUIRED) +find_package(MPI) + +select_compiler_flags(cxx_flags + GNU "-march=native" # I suggest remove "-O3" as this is controlled by the CMAKE_BUILD_TYPE + CLANG "-march=native" # same here + Intel "-axavx2,avx") +set(CXX_FLAGS ${cxx_flags}) +if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") + set(CXX_FLAGS "${CXX_FLAGS} -Wall -Wno-comment") + if(APPLE) + set(CXX_FLAGS "${CXX_FLAGS} -Wa,-q") + endif() +endif() +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${CXX_FLAGS}") + +# ================================================================================================== + +# LUD with the MKL library +if (PETSc_FOUND AND MPI_FOUND) + message("Petsc dir: ${PETSC_DIR}") + message("Petsc arch: ${PETSC_ARCH}") + message("Petsc includes: ${PETSC_INCLUDES}") + message("Petsc libs: ${PETSC_LIBRARIES}") + include_directories(${PETSC_INCLUDES}) + add_definitions (${PETSC_DEFINITIONS}) + link_directories(${petsc_lib_dir}) + + add_executable(${NAME} ksp_solver_simple.c) + target_link_libraries(${NAME} ${PETSC_LIBRARIES}) + install(TARGETS ${NAME} DESTINATION bin) + + add_executable(${NAMEB} ksp_solver_two_sys.c) + target_link_libraries(${NAMEB} ${PETSC_LIBRARIES} mpi) + install(TARGETS ${NAMEB} DESTINATION bin) + + add_executable(${NAMEC} ksp_solver_multi_rhs.c) + target_link_libraries(${NAMEC} ${PETSC_LIBRARIES} mpi) + install(TARGETS ${NAMEC} DESTINATION bin) +else () + message("## Skipping '${NAME}': no PETSC support found") + install(CODE "MESSAGE(\"${NAME} can only be built with PETSC.\")") +endif() + +unset(NAME) + +# ================================================================================================== diff --git a/krylov/README.md b/krylov/README.md new file mode 100644 index 0000000..1ea40a3 --- /dev/null +++ b/krylov/README.md @@ -0,0 +1,60 @@ +======= +README +======= +The samples in this folder demonstrate the usage of PETSc (implicit parallelism) +- Basic structures (Mat, Vec, KSP, PC, etc.) +- Profiling stages of the code +- Examples can make use of different preconditioners and solvers (provided via command line arguments) + +- 1. Code sample name +ksp_solver_simple.c +ksp_solver_two_sys.c +ksp_colver_multi_rhs.c + +- 2. Description of the code sample package + Samples for introduction to Krylov Subspace Methods by using PETSc: + - Linear system solution in parallel + - 2D Laplacian (2D mesh) + - Repeatedly solving two linear systems + - Same preconditioner + - Two different matrices with the same nonzero pattern + - Solving multiple linear systems + - Same cofficient matrix + - Different right-hand-side vectors + +- 3. Release date +20 January 2016 + +- 4. Version history +1.0: initial version + +- 5. Contributor (s) / Maintainer(s) +Cevdet Aykanat (aykanat@cs.bilkent.edu.tr) +Kadir Akbudak (kadir.cs@gmail.com) +Reha Oguz Selvitopi(reha@cs.bilkent.edu.tr) + +- 6. Copyright / License of the code sample + +- 7. Language(s) +C + +- 8. Parallelisation Implementation(s) +MPI + +- 9. Level of the code sample complexity +new starters + +- 10. Instructions on how to compile the code +*Note that PETSc and MPI libraries must be available.* + +$ rm -rf CMakeFiles CMakeCache.txt; cmake .;make + +- 11. Instructions on how to run the code +$./run.sh + +- 12. Sample input(s) +No input is required. + +- 13. Sample output(s) +No file is produced. norm are printed to standard output. + diff --git a/krylov/ksp_solver_multi_rhs.c b/krylov/ksp_solver_multi_rhs.c new file mode 100644 index 0000000..8b65695 --- /dev/null +++ b/krylov/ksp_solver_multi_rhs.c @@ -0,0 +1,221 @@ +/* +2-clause BSD license + +Copyright (c) 1991-2014, UChicago Argonne, LLC and the PETSc Development Team +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above copyright notice, this + list of conditions and the following disclaimer in the documentation and/or + other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +static char help[] = "Solves a sequence of linear systems with different right-hand-side vectors.\n\ +Input parameters include:\n\ + -ntimes : number of linear systems to solve\n\ + -view_exact_sol : write exact solution vector to stdout\n\ + -m : number of mesh points in x-direction\n\ + -n : number of mesh points in y-direction\n\n"; + +/*T + Concepts: KSP^repeatedly solving linear systems; + Concepts: KSP^Laplacian, 2d + Concepts: Laplacian, 2d + Processors: n +T*/ + +/* + Include "petscksp.h" so that we can use KSP solvers. Note that this file + automatically includes: + petscsys.h - base PETSc routines petscvec.h - vectors + petscmat.h - matrices + petscis.h - index sets petscksp.h - Krylov subspace methods + petscviewer.h - viewers petscpc.h - preconditioners +*/ +#include + +#undef __FUNCT__ +#define __FUNCT__ "main" +int main(int argc,char **args) +{ + Vec x,b,u; /* approx solution, RHS, exact solution */ + Mat A; /* linear system matrix */ + KSP ksp; /* linear solver context */ + PetscReal norm; /* norm of solution error */ + PetscErrorCode ierr; + PetscInt ntimes,i,j,k,Ii,J,Istart,Iend; + PetscInt m = 8,n = 7,its; + PetscBool flg = PETSC_FALSE; + PetscScalar v,one = 1.0,neg_one = -1.0,rhs; + + PetscInitialize(&argc,&args,(char*)0,help); + ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr); + ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr); + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Compute the matrix for use in solving a series of + linear systems of the form, A x_i = b_i, for i=1,2,... + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* + Create parallel matrix, specifying only its global dimensions. + When using MatCreate(), the matrix format can be specified at + runtime. Also, the parallel partitioning of the matrix is + determined by PETSc at runtime. + */ + ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); + ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); + ierr = MatSetFromOptions(A);CHKERRQ(ierr); + ierr = MatSetUp(A);CHKERRQ(ierr); + + /* + Currently, all PETSc parallel matrix formats are partitioned by + contiguous chunks of rows across the processors. Determine which + rows of the matrix are locally owned. + */ + ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); + + /* + Set matrix elements for the 2-D, five-point stencil in parallel. + - Each processor needs to insert only elements that it owns + locally (but any non-local elements will be sent to the + appropriate processor during matrix assembly). + - Always specify global rows and columns of matrix entries. + */ + for (Ii=Istart; Ii0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} + if (i0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,INSERT_VALUES);CHKERRQ(ierr);} + if (j -pc_type -ksp_monitor -ksp_rtol + These options will override those specified above as long as + KSPSetFromOptions() is called _after_ any other customization + routines. + */ + ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Solve several linear systems of the form A x_i = b_i + I.e., we retain the same matrix (A) for all systems, but + change the right-hand-side vector (b_i) at each step. + + In this case, we simply call KSPSolve() multiple times. The + preconditioner setup operations (e.g., factorization for ILU) + be done during the first call to KSPSolve() only; such operations + will NOT be repeated for successive solves. + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + ntimes = 2; + ierr = PetscOptionsGetInt(NULL,"-ntimes",&ntimes,NULL);CHKERRQ(ierr); + for (k=1; k + +#undef __FUNCT__ +#define __FUNCT__ "main" +int main(int argc,char **args) +{ + Vec x,b,u; /* approx solution, RHS, exact solution */ + Mat A; /* linear system matrix */ + KSP ksp; /* linear solver context */ + PetscRandom rctx; /* random number generator context */ + PetscReal norm; /* norm of solution error */ + PetscInt i,j,Ii,J,Istart,Iend,m = 8,n = 7,its; + PetscErrorCode ierr; + PetscBool flg = PETSC_FALSE; + PetscScalar v; +#if defined(PETSC_USE_LOG) + PetscLogStage stage; +#endif + + PetscInitialize(&argc,&args,(char*)0,help); + ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr); + ierr = PetscOptionsGetInt(NULL,"-n",&n,NULL);CHKERRQ(ierr); + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Compute the matrix and right-hand-side vector that define + the linear system, Ax = b. + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + /* + Create parallel matrix, specifying only its global dimensions. + When using MatCreate(), the matrix format can be specified at + runtime. Also, the parallel partitioning of the matrix is + determined by PETSc at runtime. + + Performance tuning note: For problems of substantial size, + preallocation of matrix memory is crucial for attaining good + performance. See the matrix chapter of the users manual for details. + */ + ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); + ierr = MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); + ierr = MatSetFromOptions(A);CHKERRQ(ierr); + ierr = MatMPIAIJSetPreallocation(A,5,NULL,5,NULL);CHKERRQ(ierr); + ierr = MatSeqAIJSetPreallocation(A,5,NULL);CHKERRQ(ierr); + + /* + Currently, all PETSc parallel matrix formats are partitioned by + contiguous chunks of rows across the processors. Determine which + rows of the matrix are locally owned. + */ + ierr = MatGetOwnershipRange(A,&Istart,&Iend);CHKERRQ(ierr); + + /* + Set matrix elements for the 2-D, five-point stencil in parallel. + - Each processor needs to insert only elements that it owns + locally (but any non-local elements will be sent to the + appropriate processor during matrix assembly). + - Always specify global rows and columns of matrix entries. + + Note: this uses the less common natural ordering that orders first + all the unknowns for x = h then for x = 2h etc; Hence you see J = Ii +- n + instead of J = I +- m as you might expect. The more standard ordering + would first do all variables for y = h, then y = 2h etc. + + */ + ierr = PetscLogStageRegister("Assembly", &stage);CHKERRQ(ierr); + ierr = PetscLogStagePush(stage);CHKERRQ(ierr); + for (Ii=Istart; Ii0) {J = Ii - n; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + if (i0) {J = Ii - 1; ierr = MatSetValues(A,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + if (j -pc_type -ksp_monitor -ksp_rtol + These options will override those specified above as long as + KSPSetFromOptions() is called _after_ any other customization + routines. + */ + ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Solve the linear system + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); + + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + Check solution and clean up + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + + /* + Check the error + */ + ierr = VecAXPY(x,-1.0,u);CHKERRQ(ierr); + ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); + ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); + + /* + Print convergence information. PetscPrintf() produces a single + print statement from all processes that share a communicator. + An alternative is PetscFPrintf(), which prints to a file. + */ + ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g iterations %D\n",(double)norm,its);CHKERRQ(ierr); + + /* + Free work space. All PETSc objects should be destroyed when they + are no longer needed. + */ + ierr = KSPDestroy(&ksp);CHKERRQ(ierr); + ierr = VecDestroy(&u);CHKERRQ(ierr); ierr = VecDestroy(&x);CHKERRQ(ierr); + ierr = VecDestroy(&b);CHKERRQ(ierr); ierr = MatDestroy(&A);CHKERRQ(ierr); + + /* + Always call PetscFinalize() before exiting a program. This routine + - finalizes the PETSc libraries as well as MPI + - provides summary and diagnostic information if certain runtime + options are chosen (e.g., -log_summary). + */ + ierr = PetscFinalize(); + return 0; +} diff --git a/krylov/ksp_solver_two_sys.c b/krylov/ksp_solver_two_sys.c new file mode 100644 index 0000000..de5edd6 --- /dev/null +++ b/krylov/ksp_solver_two_sys.c @@ -0,0 +1,335 @@ +/* +2-clause BSD license + +Copyright (c) 1991-2014, UChicago Argonne, LLC and the PETSc Development Team +All rights reserved. + +Redistribution and use in source and binary forms, with or without modification, +are permitted provided that the following conditions are met: + +* Redistributions of source code must retain the above copyright notice, this + list of conditions and the following disclaimer. +* Redistributions in binary form must reproduce the above copyright notice, this + list of conditions and the following disclaimer in the documentation and/or + other materials provided with the distribution. + + THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND + ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED + WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE + DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR + ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES + (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; + LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON + ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT + (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS + SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +*/ + +static char help[] = "Solves two linear systems in parallel with KSP. The code\n\ +illustrates repeated solution of linear systems with the same preconditioner\n\ +method but different matrices (having the same nonzero structure). The code\n\ +also uses multiple profiling stages. Input arguments are\n\ + -m : problem size\n\ + -mat_nonsym : use nonsymmetric matrix (default is symmetric)\n\n"; + +/*T + Concepts: KSP^repeatedly solving linear systems; + Concepts: PetscLog^profiling multiple stages of code; + Processors: n +T*/ + +/* + Include "petscksp.h" so that we can use KSP solvers. Note that this file + automatically includes: + petscsys.h - base PETSc routines petscvec.h - vectors + petscmat.h - matrices + petscis.h - index sets petscksp.h - Krylov subspace methods + petscviewer.h - viewers petscpc.h - preconditioners +*/ +#include + +#undef __FUNCT__ +#define __FUNCT__ "main" +int main(int argc,char **args) +{ + KSP ksp; /* linear solver context */ + Mat C; /* matrix */ + Vec x,u,b; /* approx solution, RHS, exact solution */ + PetscReal norm; /* norm of solution error */ + PetscScalar v,none = -1.0; + PetscInt Ii,J,ldim,low,high,iglobal,Istart,Iend; + PetscErrorCode ierr; + PetscInt i,j,m = 3,n = 2,its; + PetscMPIInt size,rank; + PetscBool mat_nonsymmetric = PETSC_FALSE; + PetscBool testnewC = PETSC_FALSE; +#if defined(PETSC_USE_LOG) + PetscLogStage stages[2]; +#endif + + PetscInitialize(&argc,&args,(char*)0,help); + ierr = PetscOptionsGetInt(NULL,"-m",&m,NULL);CHKERRQ(ierr); + ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&rank);CHKERRQ(ierr); + ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); + n = 2*size; + + /* + Set flag if we are doing a nonsymmetric problem; the default is symmetric. + */ + ierr = PetscOptionsGetBool(NULL,"-mat_nonsym",&mat_nonsymmetric,NULL);CHKERRQ(ierr); + + /* + Register two stages for separate profiling of the two linear solves. + Use the runtime option -log_summary for a printout of performance + statistics at the program's conlusion. + */ + ierr = PetscLogStageRegister("Original Solve",&stages[0]);CHKERRQ(ierr); + ierr = PetscLogStageRegister("Second Solve",&stages[1]);CHKERRQ(ierr); + + /* -------------- Stage 0: Solve Original System ---------------------- */ + /* + Indicate to PETSc profiling that we're beginning the first stage + */ + ierr = PetscLogStagePush(stages[0]);CHKERRQ(ierr); + + /* + Create parallel matrix, specifying only its global dimensions. + When using MatCreate(), the matrix format can be specified at + runtime. Also, the parallel partitioning of the matrix is + determined by PETSc at runtime. + */ + ierr = MatCreate(PETSC_COMM_WORLD,&C);CHKERRQ(ierr); + ierr = MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n);CHKERRQ(ierr); + ierr = MatSetFromOptions(C);CHKERRQ(ierr); + ierr = MatSetUp(C);CHKERRQ(ierr); + + /* + Currently, all PETSc parallel matrix formats are partitioned by + contiguous chunks of rows across the processors. Determine which + rows of the matrix are locally owned. + */ + ierr = MatGetOwnershipRange(C,&Istart,&Iend);CHKERRQ(ierr); + + /* + Set matrix entries matrix in parallel. + - Each processor needs to insert only elements that it owns + locally (but any non-local elements will be sent to the + appropriate processor during matrix assembly). + - Always specify global row and columns of matrix entries. + */ + for (Ii=Istart; Ii0) {J = Ii - n; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + if (i0) {J = Ii - 1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + if (j1) {J = Ii-n-1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + } + } else { + ierr = MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE);CHKERRQ(ierr); + ierr = MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE);CHKERRQ(ierr); + } + + /* + Assemble matrix, using the 2-step process: + MatAssemblyBegin(), MatAssemblyEnd() + Computations can be done while messages are in transition + by placing code between these two statements. + */ + ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); + ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); + + /* + Create parallel vectors. + - When using VecSetSizes(), we specify only the vector's global + dimension; the parallel partitioning is determined at runtime. + - Note: We form 1 vector from scratch and then duplicate as needed. + */ + ierr = VecCreate(PETSC_COMM_WORLD,&u);CHKERRQ(ierr); + ierr = VecSetSizes(u,PETSC_DECIDE,m*n);CHKERRQ(ierr); + ierr = VecSetFromOptions(u);CHKERRQ(ierr); + ierr = VecDuplicate(u,&b);CHKERRQ(ierr); + ierr = VecDuplicate(b,&x);CHKERRQ(ierr); + + /* + Currently, all parallel PETSc vectors are partitioned by + contiguous chunks across the processors. Determine which + range of entries are locally owned. + */ + ierr = VecGetOwnershipRange(x,&low,&high);CHKERRQ(ierr); + + /* + Set elements within the exact solution vector in parallel. + - Each processor needs to insert only elements that it owns + locally (but any non-local entries will be sent to the + appropriate processor during vector assembly). + - Always specify global locations of vector entries. + */ + ierr = VecGetLocalSize(x,&ldim);CHKERRQ(ierr); + for (i=0; i -pc_type ) + */ + ierr = KSPSetFromOptions(ksp);CHKERRQ(ierr); + + /* + Solve linear system. Here we explicitly call KSPSetUp() for more + detailed performance monitoring of certain preconditioners, such + as ICC and ILU. This call is optional, as KSPSetUp() will + automatically be called within KSPSolve() if it hasn't been + called already. + */ + ierr = KSPSetUp(ksp);CHKERRQ(ierr); + ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); + + /* + Check the error + */ + ierr = VecAXPY(x,none,u);CHKERRQ(ierr); + ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); + ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); + if (norm > 1.e-13) { + ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr); + } + + /* -------------- Stage 1: Solve Second System ---------------------- */ + /* + Solve another linear system with the same method. We reuse the KSP + context, matrix and vector data structures, and hence save the + overhead of creating new ones. + + Indicate to PETSc profiling that we're concluding the first + stage with PetscLogStagePop(), and beginning the second stage with + PetscLogStagePush(). + */ + ierr = PetscLogStagePop();CHKERRQ(ierr); + ierr = PetscLogStagePush(stages[1]);CHKERRQ(ierr); + + /* + Initialize all matrix entries to zero. MatZeroEntries() retains the + nonzero structure of the matrix for sparse formats. + */ + ierr = MatZeroEntries(C);CHKERRQ(ierr); + + /* + Assemble matrix again. Note that we retain the same matrix data + structure and the same nonzero pattern; we just change the values + of the matrix entries. + */ + for (i=0; i0) {J = Ii - n; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + if (i0) {J = Ii - 1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + if (j1) {J = Ii-n-1; ierr = MatSetValues(C,1,&Ii,1,&J,&v,ADD_VALUES);CHKERRQ(ierr);} + } + } + ierr = MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); + ierr = MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); + + ierr = PetscOptionsGetBool(NULL,"-test_newMat",&testnewC,NULL);CHKERRQ(ierr); + if (testnewC) { + /* + User may use a new matrix C with same nonzero pattern, e.g. + ./ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_package mumps -test_newMat + */ + Mat Ctmp; + ierr = MatDuplicate(C,MAT_COPY_VALUES,&Ctmp);CHKERRQ(ierr); + ierr = MatDestroy(&C);CHKERRQ(ierr); + ierr = MatDuplicate(Ctmp,MAT_COPY_VALUES,&C);CHKERRQ(ierr); + ierr = MatDestroy(&Ctmp);CHKERRQ(ierr); + } + /* + Compute another right-hand-side vector + */ + ierr = MatMult(C,u,b);CHKERRQ(ierr); + + /* + Set operators. Here the matrix that defines the linear system + also serves as the preconditioning matrix. + */ + ierr = KSPSetOperators(ksp,C,C);CHKERRQ(ierr); + + /* + Solve linear system + */ + ierr = KSPSetUp(ksp);CHKERRQ(ierr); + ierr = KSPSolve(ksp,b,x);CHKERRQ(ierr); + + /* + Check the error + */ + ierr = VecAXPY(x,none,u);CHKERRQ(ierr); + ierr = VecNorm(x,NORM_2,&norm);CHKERRQ(ierr); + ierr = KSPGetIterationNumber(ksp,&its);CHKERRQ(ierr); + if (norm > 1.e-4) { + ierr = PetscPrintf(PETSC_COMM_WORLD,"Norm of error %g, Iterations %D\n",(double)norm,its);CHKERRQ(ierr); + } + + /* + Free work space. All PETSc objects should be destroyed when they + are no longer needed. + */ + ierr = KSPDestroy(&ksp);CHKERRQ(ierr); + ierr = VecDestroy(&u);CHKERRQ(ierr); + ierr = VecDestroy(&x);CHKERRQ(ierr); + ierr = VecDestroy(&b);CHKERRQ(ierr); + ierr = MatDestroy(&C);CHKERRQ(ierr); + + /* + Indicate to PETSc profiling that we're concluding the second stage + */ + ierr = PetscLogStagePop();CHKERRQ(ierr); + + ierr = PetscFinalize(); + return 0; +} + + diff --git a/krylov/run.sh b/krylov/run.sh new file mode 100755 index 0000000..27d49d8 --- /dev/null +++ b/krylov/run.sh @@ -0,0 +1,11 @@ +mpirun -n 1 ./2_sparse_krylov_simple -ksp_monitor_short -m 5 -n 5 -ksp_gmres_cgs_refinement_type refine_always +mpirun -n 2 ./2_sparse_krylov_simple -ksp_monitor_short -m 5 -n 5 -mat_view draw -ksp_gmres_cgs_refinement_type refine_always -nox +mpirun -n 4 ./2_sparse_krylov_simple -pc_type bjacobi -pc_bjacobi_blocks 1 -ksp_monitor_short -sub_pc_type jacobi -sub_ksp_type gmres +mpirun -n 3 ./2_sparse_krylov_simple -ksp_type fbcgsr -pc_type bjacobi + +mpirun -n 1 ./2_sparse_krylov_twosys -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always +mpirun -n 2 ./2_sparse_krylov_twosys -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001 +mpirun -n 2 ./2_sparse_krylov_twosys -ksp_gmres_cgs_refinement_type refine_always +mpirun -n 5 ./2_sparse_krylov_twosys -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi -ksp_rtol 1.e-10 + +mpirun -n 2 ./2_sparse_krylov_multirhs -ntimes 4 -ksp_gmres_cgs_refinement_type refine_always -- GitLab