! ! Description: Solves a complex linear system in parallel with KSP (Fortran code). ! !!/*T ! Concepts: KSP^solving a Helmholtz equation ! Concepts: complex numbers ! Processors: n !T*/ ! ! The model problem: ! Solve Helmholtz equation on the unit square: (0,1) x (0,1) ! -delta u - sigma1*u + i*sigma2*u = f, ! where delta = Laplace operator ! Dirichlet b.c.'s on all sides ! Use the 2-D, five-point finite difference stencil. ! ! Compiling the code: ! This code uses the complex numbers version of PETSc, so configure ! must be run to enable this ! ! ! ----------------------------------------------------------------------- program main #include use petscksp implicit none ! ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Variable declarations ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! ! Variables: ! ksp - linear solver context ! x, b, u - approx solution, right-hand-side, exact solution vectors ! A - matrix that defines linear system ! its - iterations for convergence ! norm - norm of error in solution ! rctx - random number context ! KSP ksp Mat A Vec x,b,u PetscRandom rctx PetscReal norm,h2,sigma1 PetscScalar none,sigma2,v,pfive,czero PetscScalar cone PetscInt dim,its,n,Istart PetscInt Iend,i,j,II,JJ,one PetscErrorCode ierr PetscMPIInt rank PetscBool flg integer irun,nruns integer comm_local, size_local, rank_local logical use_random, use_cuda ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Beginning of program ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - none = -1.0 n = 6 sigma1 = 100.0 czero = 0.0 nruns = 3 call MPI_INIT(ierr) comm_local = MPI_COMM_WORLD call MPI_COMM_RANK( comm_local, rank_local, ierr) call MPI_COMM_SIZE( comm_local, size_local, ierr) write(6,*) 'hello from', rank_local, ' of ', size_local PETSC_COMM_WORLD = comm_local do irun = 1,nruns call PetscInitialize(PETSC_NULL_CHARACTER,ierr) write(6,*) 'petsc initalized for ', irun, ' time(s)' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local if (ierr .ne. 0) then print*,'Unable to initialize PETSc' stop endif CHKERRA(ierr) cone = PETSC_i call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr) call PetscOptionsGetReal(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,& & '-sigma1',sigma1,flg,ierr) call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER, & & '-n',n,flg,ierr) dim = n*n ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! 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. call MatCreate(PETSC_COMM_WORLD,A,ierr) call MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim,ierr) call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,& & '-usecuda',flg,ierr) if (flg) then use_cuda = .true. call MatSetType(A,MATMPIAIJCUSPARSE,ierr) write(6,*) 'set mat for gpu' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) write(6,*) 'matsettype' else use_cuda = .false. ! just need this to fall back ! to CPUs call MatSetType(A,MATMPIAIJ,ierr) endif call MatSetFromOptions(A,ierr) call MatSetUp(A,ierr) write(6,*) 'mat setup' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(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. call MatGetOwnershipRange(A,Istart,Iend,ierr) write(6,*) 'mat ownership range' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) ! Set matrix elements 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. call PetscOptionsHasName(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,& & '-norandom',flg,ierr) if (flg) then use_random = .false. sigma2 = 10.0*PETSC_i else use_random = .true. call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr) call PetscRandomSetFromOptions(rctx,ierr) call PetscRandomSetInterval(rctx,czero,cone,ierr) endif h2 = 1.0/real((n+1)*(n+1)) one = 1 do 10, II=Istart,Iend-1 v = -1.0 i = II/n j = II - i*n if (i.gt.0) then JJ = II - n call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (i.lt.n-1) then JJ = II + n call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (j.gt.0) then JJ = II - 1 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (j.lt.n-1) then JJ = II + 1 call MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr) endif if (use_random) call PetscRandomGetValue(rctx,sigma2,ierr) v = 4.0 - sigma1*h2 + sigma2*h2 call MatSetValues(A,one,II,one,II,v,ADD_VALUES,ierr) 10 continue if (use_random) call PetscRandomDestroy(rctx,ierr) write(6,*) 'mat value set' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(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. call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr) call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr) write(6,*) 'mat assumbly' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) ! Create parallel vectors. ! - Here, the parallel partitioning of the vector is determined by ! PETSc at runtime. We could also specify the local dimensions ! if desired. ! - Note: We form 1 vector from scratch and then duplicate as needed. ! call VecCreate(PETSC_COMM_WORLD,u,ierr) call MatCreateVecs(A, u, PETSC_NULL_VEC, ierr) ! call VecSetSizes(u,PETSC_DECIDE,dim,ierr) write(6,*) 'setup vec from mat' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) call VecSetFromOptions(u,ierr) call VecDuplicate(u,b,ierr) call VecDuplicate(b,x,ierr) write(6,*) 'vec duplication' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) ! Set exact solution; then compute right-hand-side vector. if (use_random) then call PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr) call PetscRandomSetFromOptions(rctx,ierr) call VecSetRandom(u,rctx,ierr) else pfive = 0.5 call VecSet(u,pfive,ierr) endif call MatMult(A,u,b,ierr) write(6,*) 'MatMult for rhs' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create the linear solver and set various options ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Create linear solver context call KSPCreate(PETSC_COMM_WORLD,ksp,ierr) ! Set operators. Here the matrix that defines the linear system ! also serves as the preconditioning matrix. call KSPSetOperators(ksp,A,A,ierr) ! Set runtime options, e.g., ! -ksp_type -pc_type -ksp_monitor -ksp_rtol call KSPSetFromOptions(ksp,ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Solve the linear system ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - call KSPSolve(ksp,b,x,ierr) write(6,*) 'KSP solve the system' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Check solution and clean up ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - ! Check the error call VecAXPY(x,none,u,ierr) write(6,*) 'AxPy for error estimation' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) call VecNorm(x,NORM_2,norm,ierr) call KSPGetIterationNumber(ksp,its,ierr) if (rank .eq. 0) then write(6,*) 'run ', irun if (norm .gt. 1.e-12) then write(6,100) norm,its else write(6,110) its endif endif 100 format('Norm of error ',e11.4,',iterations ',i5) 110 format('Norm of error < 1.e-12,iterations ',i5) ! Free work space. All PETSc objects should be destroyed when they ! are no longer needed. if (use_random) then call PetscRandomDestroy(rctx,ierr) end if call KSPDestroy(ksp,ierr) write(6,*) 'KSP distroyed' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) call VecDestroy(u,ierr) call VecDestroy(x,ierr) call VecDestroy(b,ierr) write(6,*) 'vecs distroyed' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) call MatDestroy(A,ierr) write(6,*) 'mat distroyed' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) call PetscFinalize(ierr) write(6,*) 'PETSc finalized' write(6,*) 'with errcode= ', ierr, ' @ ', rank_local CHKERRA(ierr) enddo write(6,*) 'goodbye from', rank_local, ' of ', size_local call MPI_FINALIZE(ierr) end ! !/*TEST ! ! build: ! requires: complex ! ! test: ! args: -n 6 -norandom -pc_type none -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always ! output_file: output/ex11f_1.out ! !TEST*/