1! 2! Description: Solves a complex linear system in parallel with KSP (Fortran code). 3! 4!!/*T 5! Concepts: KSP^solving a Helmholtz equation 6! Concepts: complex numbers 7! Processors: n 8!T*/ 9 10 11! 12! The model problem: 13! Solve Helmholtz equation on the unit square: (0,1) x (0,1) 14! -delta u - sigma1*u + i*sigma2*u = f, 15! where delta = Laplace operator 16! Dirichlet b.c.'s on all sides 17! Use the 2-D, five-point finite difference stencil. 18! 19! Compiling the code: 20! This code uses the complex numbers version of PETSc, so configure 21! must be run to enable this 22! 23! 24! ----------------------------------------------------------------------- 25 26 program main 27#include <petsc/finclude/petscksp.h> 28 use petscksp 29 implicit none 30 31! 32! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 33! Variable declarations 34! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 35! 36! Variables: 37! ksp - linear solver context 38! x, b, u - approx solution, right-hand-side, exact solution vectors 39! A - matrix that defines linear system 40! its - iterations for convergence 41! norm - norm of error in solution 42! rctx - random number context 43! 44 45 KSP ksp 46 Mat A 47 Vec x,b,u 48 PetscRandom rctx 49 PetscReal norm,h2,sigma1 50 PetscScalar none,sigma2,v,pfive,czero 51 PetscScalar cone 52 PetscInt dim,its,n,Istart 53 PetscInt Iend,i,j,II,JJ,one 54 PetscErrorCode ierr 55 PetscMPIInt rank 56 PetscBool flg 57 logical use_random 58 59! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 60! Beginning of program 61! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 62 63 call PetscInitialize(PETSC_NULL_CHARACTER,ierr) 64 if (ierr .ne. 0) then 65 print*,'Unable to initialize PETSc' 66 stop 67 endif 68 69 none = -1.0 70 n = 6 71 sigma1 = 100.0 72 czero = 0.0 73 cone = PETSC_i 74 call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) 75 call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-sigma1',sigma1,flg,ierr) 76 call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr) 77 dim = n*n 78 79! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 80! Compute the matrix and right-hand-side vector that define 81! the linear system, Ax = b. 82! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 83 84! Create parallel matrix, specifying only its global dimensions. 85! When using MatCreate(), the matrix format can be specified at 86! runtime. Also, the parallel partitioning of the matrix is 87! determined by PETSc at runtime. 88 89 call MatCreate(PETSC_COMM_WORLD,A,ierr) 90 call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim,ierr) 91 call MatSetFromOptions(A,ierr) 92 call MatSetUp(A,ierr) 93 94! Currently, all PETSc parallel matrix formats are partitioned by 95! contiguous chunks of rows across the processors. Determine which 96! rows of the matrix are locally owned. 97 98 call MatGetOwnershipRange(A,Istart,Iend,ierr) 99 100! Set matrix elements in parallel. 101! - Each processor needs to insert only elements that it owns 102! locally (but any non-local elements will be sent to the 103! appropriate processor during matrix assembly). 104! - Always specify global rows and columns of matrix entries. 105 106 call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-norandom',flg,ierr) 107 if (flg) then 108 use_random = .false. 109 sigma2 = 10.0*PETSC_i 110 else 111 use_random = .true. 112 call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr) 113 call PetscRandomSetFromOptions(rctx,ierr) 114 call PetscRandomSetInterval(rctx,czero,cone,ierr) 115 endif 116 h2 = 1.0/real((n+1)*(n+1)) 117 118 one = 1 119 do 10, II=Istart,Iend-1 120 v = -1.0 121 i = II/n 122 j = II - i*n 123 if (i.gt.0) then 124 JJ = II - n 125 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) 126 endif 127 if (i.lt.n-1) then 128 JJ = II + n 129 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) 130 endif 131 if (j.gt.0) then 132 JJ = II - 1 133 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) 134 endif 135 if (j.lt.n-1) then 136 JJ = II + 1 137 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) 138 endif 139 if (use_random) call PetscRandomGetValue(rctx,sigma2,ierr) 140 v = 4.0 - sigma1*h2 + sigma2*h2 141 call MatSetValues(A,one,II,one,II,v,ADD_VALUES,ierr) 142 10 continue 143 if (use_random) call PetscRandomDestroy(rctx,ierr) 144 145! Assemble matrix, using the 2-step process: 146! MatAssemblyBegin(), MatAssemblyEnd() 147! Computations can be done while messages are in transition 148! by placing code between these two statements. 149 150 call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) 151 call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) 152 153! Create parallel vectors. 154! - Here, the parallel partitioning of the vector is determined by 155! PETSc at runtime. We could also specify the local dimensions 156! if desired. 157! - Note: We form 1 vector from scratch and then duplicate as needed. 158 159 call VecCreate(PETSC_COMM_WORLD,u,ierr) 160 call VecSetSizes(u,PETSC_DECIDE,dim,ierr) 161 call VecSetFromOptions(u,ierr) 162 call VecDuplicate(u,b,ierr) 163 call VecDuplicate(b,x,ierr) 164 165! Set exact solution; then compute right-hand-side vector. 166 167 if (use_random) then 168 call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr) 169 call PetscRandomSetFromOptions(rctx,ierr) 170 call VecSetRandom(u,rctx,ierr) 171 else 172 pfive = 0.5 173 call VecSet(u,pfive,ierr) 174 endif 175 call MatMult(A,u,b,ierr) 176 177! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178! Create the linear solver and set various options 179! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 180 181! Create linear solver context 182 183 call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) 184 185! Set operators. Here the matrix that defines the linear system 186! also serves as the preconditioning matrix. 187 188 call KSPSetOperators(ksp,A,A,ierr) 189 190! Set runtime options, e.g., 191! -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol> 192 193 call KSPSetFromOptions(ksp,ierr) 194 195! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 196! Solve the linear system 197! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 198 199 call KSPSolve(ksp,b,x,ierr) 200 201! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202! Check solution and clean up 203! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 204 205! Check the error 206 207 call VecAXPY(x,none,u,ierr) 208 call VecNorm(x,NORM_2,norm,ierr) 209 call KSPGetIterationNumber(ksp,its,ierr) 210 if (rank .eq. 0) then 211 if (norm .gt. 1.e-12) then 212 write(6,100) norm,its 213 else 214 write(6,110) its 215 endif 216 endif 217 100 format('Norm of error ',e11.4,',iterations ',i5) 218 110 format('Norm of error < 1.e-12,iterations ',i5) 219 220! Free work space. All PETSc objects should be destroyed when they 221! are no longer needed. 222 223 if (use_random) call PetscRandomDestroy(rctx,ierr) 224 call KSPDestroy(ksp,ierr) 225 call VecDestroy(u,ierr) 226 call VecDestroy(x,ierr) 227 call VecDestroy(b,ierr) 228 call MatDestroy(A,ierr) 229 230 call PetscFinalize(ierr) 231 end 232 233! 234!/*TEST 235! 236! build: 237! requires: complex 238! 239! test: 240! args: -n 6 -norandom -pc_type none -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always 241! output_file: output/ex11f_1.out 242! 243!TEST*/ 244