Actual source code: ex7f.F90

      2: ! Block Jacobi preconditioner for solving a linear system in parallel with KSP
      3: ! The code indicates the procedures for setting the particular block sizes and
      4: ! for using different linear solvers on the individual blocks

      6: ! This example focuses on ways to customize the block Jacobi preconditioner.
      7: ! See ex1.c and ex2.c for more detailed comments on the basic usage of KSP
      8: ! (including working with matrices and vectors)

     10: ! Recall: The block Jacobi method is equivalent to the ASM preconditioner with zero overlap.

     12: !/*T
     13: !   Concepts: KSP^customizing the block Jacobi preconditioner
     14: !   Processors: n
     15: !T*/

     18: program main
     19: #include <petsc/finclude/petscksp.h>
     20:       use petscksp

     22:       implicit none
     23:       Vec             :: x,b,u      ! approx solution, RHS, exact solution
     24:       Mat             :: A            ! linear system matrix
     25:       KSP             :: ksp         ! KSP context
     26:       PC              :: myPc           ! PC context
     27:       PC              :: subpc        ! PC context for subdomain
     28:       PetscReal       :: norm         ! norm of solution error
     29:       PetscReal,parameter :: tol = 1.e-6
     30:       PetscErrorCode  :: ierr
     31:       PetscInt        :: i,j,Ii,JJ,n
     32:       PetscInt, parameter :: m = 4
     33:       PetscMPIInt     :: myRank,mySize
     34:       PetscInt        :: its,nlocal,first,Istart,Iend
     35:       PetscScalar     :: v
     36:       PetscScalar, parameter :: &
     37:         myNone = -1.0, &
     38:         sone   = 1.0
     39:       PetscBool       :: isbjacobi,flg
     40:       KSP,allocatable,dimension(:)      ::   subksp     ! array of local KSP contexts on this processor
     41:       PetscInt,allocatable,dimension(:) :: blks
     42:       character(len=PETSC_MAX_PATH_LEN) :: outputString
     43:       PetscInt,parameter :: one = 1, five = 5

     45:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
     46:       if (ierr /= 0) then
     47:         write(6,*)'Unable to initialize PETSc'
     48:         stop
     49:       endif

     51:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
     52:       CHKERRA(ierr)
     53:       call MPI_Comm_rank(PETSC_COMM_WORLD,myRank,ierr)
     54:       CHKERRA(ierr)
     55:       call MPI_Comm_size(PETSC_COMM_WORLD,mySize,ierr)
     56:       CHKERRA(ierr)
     57:       n=m+2

     59:       !-------------------------------------------------------------------
     60:       ! Compute the matrix and right-hand-side vector that define
     61:       ! the linear system, Ax = b.
     62:       !---------------------------------------------------------------

     64:       ! Create and assemble parallel matrix

     66:       call  MatCreate(PETSC_COMM_WORLD,A,ierr)
     67:       call  MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
     68:       call  MatSetFromOptions(A,ierr)
     69:       call  MatMPIAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,five,PETSC_NULL_INTEGER,ierr)
     70:       call  MatSeqAIJSetPreallocation(A,five,PETSC_NULL_INTEGER,ierr)
     71:       call  MatGetOwnershipRange(A,Istart,Iend,ierr)

     73:       do Ii=Istart,Iend-1
     74:           v =-1.0; i = Ii/n; j = Ii - i*n
     75:           if (i>0) then
     76:             JJ = Ii - n
     77:             call MatSetValues(A,one,Ii,one,JJ,v,ADD_VALUES,ierr);CHKERRA(ierr)
     78:           endif

     80:           if (i<m-1) then
     81:             JJ = Ii + n
     82:             call MatSetValues(A,one,Ii,one,JJ,v,ADD_VALUES,ierr);CHKERRA(ierr)
     83:           endif

     85:           if (j>0) then
     86:             JJ = Ii - 1
     87:             call MatSetValues(A,one,Ii,one,JJ,v,ADD_VALUES,ierr);CHKERRA(ierr)
     88:           endif

     90:           if (j<n-1) then
     91:             JJ = Ii + 1
     92:             call MatSetValues(A,one,Ii,one,JJ,v,ADD_VALUES,ierr);CHKERRA(ierr)
     93:           endif

     95:           v=4.0
     96:           call MatSetValues(A,one,Ii,one,Ii,v,ADD_VALUES,ierr);CHKERRA(ierr)

     98:         enddo

    100:       call MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)
    101:       call MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr);CHKERRA(ierr)

    103:       ! Create parallel vectors

    105:       call  VecCreate(PETSC_COMM_WORLD,u,ierr)
    106:       CHKERRA(ierr)
    107:       call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
    108:       CHKERRA(ierr)
    109:       call  VecSetFromOptions(u,ierr)
    110:       CHKERRA(ierr)
    111:       call  VecDuplicate(u,b,ierr)
    112:       call  VecDuplicate(b,x,ierr)

    114:       ! Set exact solution; then compute right-hand-side vector.

    116:       call Vecset(u,sone,ierr)
    117:       CHKERRA(ierr)
    118:       call MatMult(A,u,b,ierr)
    119:       CHKERRA(ierr)

    121:       ! Create linear solver context

    123:       call KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
    124:       CHKERRA(ierr)

    126:       ! Set operators. Here the matrix that defines the linear system
    127:       ! also serves as the preconditioning matrix.

    129:       call KSPSetOperators(ksp,A,A,ierr)
    130:       CHKERRA(ierr)

    132:       ! Set default preconditioner for this program to be block Jacobi.
    133:       ! This choice can be overridden at runtime with the option
    134:       ! -pc_type <type>

    136:       call KSPGetPC(ksp,myPc,ierr)
    137:       CHKERRA(ierr)
    138:       call PCSetType(myPc,PCBJACOBI,ierr)
    139:       CHKERRA(ierr)

    141:       ! -----------------------------------------------------------------
    142:       !            Define the problem decomposition
    143:       !-------------------------------------------------------------------

    145:       ! Call PCBJacobiSetTotalBlocks() to set individually the size of
    146:       ! each block in the preconditioner.  This could also be done with
    147:       ! the runtime option -pc_bjacobi_blocks <blocks>
    148:       ! Also, see the command PCBJacobiSetLocalBlocks() to set the
    149:       ! local blocks.

    151:       ! Note: The default decomposition is 1 block per processor.

    153:       allocate(blks(m),source = n)

    155:       call PCBJacobiSetTotalBlocks(myPc,m,blks,ierr)
    156:       CHKERRA(ierr)
    157:       deallocate(blks)

    159:       !-------------------------------------------------------------------
    160:       !       Set the linear solvers for the subblocks
    161:       !-------------------------------------------------------------------

    163:       !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    164:       ! Basic method, should be sufficient for the needs of most users.
    165:       !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    166:       ! By default, the block Jacobi method uses the same solver on each
    167:       ! block of the problem.  To set the same solver options on all blocks,
    168:       ! use the prefix -sub before the usual PC and KSP options, e.g.,
    169:       ! -sub_pc_type <pc> -sub_ksp_type <ksp> -sub_ksp_rtol 1.e-4

    171:       !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    172:       !  Advanced method, setting different solvers for various blocks.
    173:       !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

    175:       ! Note that each block's KSP context is completely independent of
    176:       ! the others, and the full range of uniprocessor KSP options is
    177:       ! available for each block. The following section of code is intended
    178:       ! to be a simple illustration of setting different linear solvers for
    179:       ! the individual blocks.  These choices are obviously not recommended
    180:       ! for solving this particular problem.

    182:       call PetscObjectTypeCompare(myPc,PCBJACOBI,isbjacobi,ierr)

    185:       if (isbjacobi) then

    187:         ! Call KSPSetUp() to set the block Jacobi data structures (including
    188:         ! creation of an internal KSP context for each block).
    189:         ! Note: KSPSetUp() MUST be called before PCBJacobiGetSubKSP()

    191:         call KSPSetUp(ksp,ierr)

    193:         ! Extract the array of KSP contexts for the local blocks
    194:         call PCBJacobiGetSubKSP(myPc,nlocal,first,PETSC_NULL_KSP,ierr)
    195:         allocate(subksp(nlocal))
    196:         call PCBJacobiGetSubKSP(myPc,nlocal,first,subksp,ierr)

    199:         ! Loop over the local blocks, setting various KSP options for each block

    201:         do i=0,nlocal-1

    203:           call KSPGetPC(subksp(i+1),subpc,ierr)

    205:           if (myRank>0) then

    207:             if (mod(i,2)==1) then
    208:               call PCSetType(subpc,PCILU,ierr); CHKERRA(ierr)

    210:             else
    211:               call PCSetType(subpc,PCNONE,ierr); CHKERRA(ierr)
    212:               call KSPSetType(subksp(i+1),KSPBCGS,ierr); CHKERRA(ierr)
    213:               call KSPSetTolerances(subksp(i+1),tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    214:               CHKERRA(ierr)
    215:             endif

    217:           else
    218:             call PCSetType(subpc,PCJACOBI,ierr); CHKERRA(ierr)
    219:             call KSPSetType(subksp(i+1),KSPGMRES,ierr); CHKERRA(ierr)
    220:             call KSPSetTolerances(subksp(i+1),tol,PETSC_DEFAULT_REAL,PETSC_DEFAULT_REAL,PETSC_DEFAULT_INTEGER,ierr)
    221:             CHKERRA(ierr)
    222:           endif

    224:         end do

    226:       endif

    228:       !----------------------------------------------------------------
    229:       !                Solve the linear system
    230:       !-----------------------------------------------------------------

    232:       ! Set runtime options

    234:       call KSPSetFromOptions(ksp,ierr); CHKERRA(ierr)

    236:       ! Solve the linear system

    238:       call KSPSolve(ksp,b,x,ierr); CHKERRA(ierr)

    240:       !  -----------------------------------------------------------------
    241:       !               Check solution and clean up
    242:       !-------------------------------------------------------------------

    245:       !  -----------------------------------------------------------------
    246:       ! Check the error
    247:       !  -----------------------------------------------------------------

    249:       !call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)

    251:       call VecAXPY(x,myNone,u,ierr)

    253:       !call VecView(x,PETSC_VIEWER_STDOUT_WORLD,ierr)

    256:       call VecNorm(x,NORM_2,norm,ierr); CHKERRA(ierr)
    257:       call KSPGetIterationNumber(ksp,its,ierr); CHKERRA(ierr)
    258:       write(outputString,*)'Norm of error',real(norm),'Iterations',its,'\n'         ! PETScScalar might be of complex type
    259:       call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr); CHKERRA(ierr)

    261:       ! Free work space.  All PETSc objects should be destroyed when they
    262:       ! are no longer needed.
    263:       deallocate(subksp)
    264:       call KSPDestroy(ksp,ierr);CHKERRA(ierr)
    265:       call VecDestroy(u,ierr); CHKERRA(ierr)
    266:       call VecDestroy(b,ierr); CHKERRA(ierr)
    267:       call MatDestroy(A,ierr); CHKERRA(ierr)
    268:       call VecDestroy(x,ierr); CHKERRA(ierr)
    269:       call PetscFinalize(ierr); CHKERRA(ierr)

    271: end program main

    273: !/*TEST
    274: !
    275: !   test:
    276: !      nsize: 2
    277: !      args: -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
    278: !
    279: !   test:
    280: !      suffix: 2
    281: !      nsize: 2
    282: !      args: -ksp_view ::ascii_info_detail
    283: !
    284: !TEST*/