1! 2! Description: Solves a tridiagonal linear system with KSP. 3! 4! ----------------------------------------------------------------------- 5! 6! Demonstrate a custom KSP convergence test that calls the default convergence test 7! 8#include <petsc/finclude/petscksp.h> 9module ex1fmodule 10 use petscksp 11 implicit none 12 13contains 14 subroutine MyKSPConverged(ksp, n, rnorm, flag, defaultctx, ierr) 15 16 KSP ksp 17 PetscErrorCode ierr 18 PetscInt n 19 integer(8) defaultctx 20 KSPConvergedReason flag 21 PetscReal rnorm 22 23 ! Must call default convergence test on the 0th iteration 24 PetscCall(KSPConvergedDefault(ksp, n, rnorm, flag, defaultctx, ierr)) 25 end subroutine MyKSPConverged 26end module ex1fmodule 27 28program main 29 use petscksp 30 use ex1fmodule 31 implicit none 32 33! 34! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 35! Variable declarations 36! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 37! 38! Variables: 39! ksp - linear solver context 40! ksp - Krylov subspace method context 41! pc - preconditioner context 42! x, b, u - approx solution, right-hand side, exact solution vectors 43! A - matrix that defines linear system 44! its - iterations for convergence 45! norm - norm of error in solution 46! 47 Vec x, b, u 48 Mat A 49 KSP ksp 50 PC pc 51 PetscReal norm, tol 52 PetscErrorCode ierr 53 PetscInt i, n, col(3), its, i1, i2, i3 54 PetscBool flg 55 PetscMPIInt size 56 PetscScalar none, one, value(3) 57 PetscLogStage stages(2) 58 PetscFortranAddr defaultctx 59 60! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 61! Beginning of program 62! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 63 64 PetscCallA(PetscInitialize(ierr)) 65 PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr)) 66 PetscCheckA(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only') 67 none = -1.0 68 one = 1.0 69 n = 10 70 i1 = 1 71 i2 = 2 72 i3 = 3 73 PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr)) 74 75 PetscCallA(PetscLogStageRegister('MatVec Assembly', stages(1), ierr)) 76 PetscCallA(PetscLogStageRegister('KSP Solve', stages(2), ierr)) 77 PetscCallA(PetscLogStagePush(stages(1), ierr)) 78! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 79! Compute the matrix and right-hand-side vector that define 80! the linear system, Ax = b. 81! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 82 83! Create matrix. When using MatCreate(), the matrix format can 84! be specified at runtime. 85 86 PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr)) 87 PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr)) 88 PetscCallA(MatSetFromOptions(A, ierr)) 89 PetscCallA(MatSetUp(A, ierr)) 90 91! Assemble matrix. 92! - Note that MatSetValues() uses 0-based row and column numbers 93! in Fortran as well as in C (as set here in the array "col"). 94 95 value(1) = -1.0 96 value(2) = 2.0 97 value(3) = -1.0 98 do i = 1, n - 2 99 col(1) = i - 1 100 col(2) = i 101 col(3) = i + 1 102 PetscCallA(MatSetValues(A, i1, [i], i3, col, value, INSERT_VALUES, ierr)) 103 end do 104 i = n - 1 105 col(1) = n - 2 106 col(2) = n - 1 107 PetscCallA(MatSetValues(A, i1, [i], i2, col, value, INSERT_VALUES, ierr)) 108 i = 0 109 col(1) = 0 110 col(2) = 1 111 value(1) = 2.0 112 value(2) = -1.0 113 PetscCallA(MatSetValues(A, i1, [i], i2, col, value, INSERT_VALUES, ierr)) 114 PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr)) 115 PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr)) 116 117! Create vectors. Note that we form 1 vector from scratch and 118! then duplicate as needed. 119 120 PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr)) 121 PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr)) 122 PetscCallA(VecSetFromOptions(x, ierr)) 123 PetscCallA(VecDuplicate(x, b, ierr)) 124 PetscCallA(VecDuplicate(x, u, ierr)) 125 126! Set exact solution; then compute right-hand-side vector. 127 128 PetscCallA(VecSet(u, one, ierr)) 129 PetscCallA(MatMult(A, u, b, ierr)) 130 PetscCallA(PetscLogStagePop(ierr)) 131 PetscCallA(PetscLogStagePush(stages(2), ierr)) 132 133! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 134! Create the linear solver and set various options 135! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 136 137! Create linear solver context 138 139 PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr)) 140 141! Set operators. Here the matrix that defines the linear system 142! also serves as the matrix from which the preconditioner is constructed. 143 144 PetscCallA(KSPConvergedDefaultCreate(defaultctx, ierr)) 145 ! use KSPConvergedDefaultDestroyCptr() since it has that is the Fortran interface definition 146 PetscCallA(KSPSetConvergenceTest(ksp, MyKSPConverged, defaultctx, KSPConvergedDefaultDestroyCptr, ierr)) 147 PetscCallA(KSPSetOperators(ksp, A, A, ierr)) 148 149! Set linear solver defaults for this problem (optional). 150! - By extracting the KSP and PC contexts from the KSP context, 151! we can then directly call any KSP and PC routines 152! to set various options. 153! - The following four statements are optional; all of these 154! parameters could alternatively be specified at runtime via 155! KSPSetFromOptions() 156 157 PetscCallA(KSPGetPC(ksp, pc, ierr)) 158 PetscCallA(PCSetType(pc, PCJACOBI, ierr)) 159 tol = .0000001 160 PetscCallA(KSPSetTolerances(ksp, tol, PETSC_CURRENT_REAL, PETSC_CURRENT_REAL, PETSC_CURRENT_INTEGER, ierr)) 161 PetscCallA(KSPGetTolerances(ksp, PETSC_NULL_REAL, tol, PETSC_NULL_REAL, PETSC_NULL_INTEGER, ierr)) 162 163! Set runtime options, e.g., 164! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> 165! These options will override those specified above as long as 166! KSPSetFromOptions() is called _after_ any other customization 167! routines. 168 169 PetscCallA(KSPSetFromOptions(ksp, ierr)) 170 171! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 172! Solve the linear system 173! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 174 175 PetscCallA(KSPSolve(ksp, b, x, ierr)) 176 PetscCallA(PetscLogStagePop(ierr)) 177 178! View solver converged reason; we could instead use the option -ksp_converged_reason 179 PetscCallA(KSPConvergedReasonView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr)) 180 181! View solver info; we could instead use the option -ksp_view 182 183 PetscCallA(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr)) 184 185! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 186! Check solution and clean up 187! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 188 189! Check the error 190 191 PetscCallA(VecAXPY(x, none, u, ierr)) 192 PetscCallA(VecNorm(x, NORM_2, norm, ierr)) 193 PetscCallA(KSPGetIterationNumber(ksp, its, ierr)) 194 if (norm > 1.e-12) then 195 write (6, 100) norm, its 196 else 197 write (6, 200) its 198 end if 199100 format('Norm of error ', e11.4, ', Iterations = ', i5) 200200 format('Norm of error < 1.e-12, Iterations = ', i5) 201 202! Free work space. All PETSc objects should be destroyed when they 203! are no longer needed. 204 205 PetscCallA(VecDestroy(x, ierr)) 206 PetscCallA(VecDestroy(u, ierr)) 207 PetscCallA(VecDestroy(b, ierr)) 208 PetscCallA(MatDestroy(A, ierr)) 209 PetscCallA(KSPDestroy(ksp, ierr)) 210 PetscCallA(PetscFinalize(ierr)) 211 212end 213 214!/*TEST 215! 216! test: 217! requires: !single 218! args: -ksp_monitor_short 219! 220!TEST*/ 221