1: !
  2: !   Tests PCMGSetResidual
  3: !
  4: ! -----------------------------------------------------------------------
  5: #include <petsc/finclude/petscksp.h>
  6: module ex8fmodule
  7:   use petscksp
  8:   implicit none

 10: contains
 11:   subroutine MyResidual(A, b, x, r, ierr)
 12:     Mat A
 13:     Vec b, x, r
 14:     integer ierr
 15:   end

 17: end module ex8fmodule

 19: program main
 20:   use petscksp
 21:   use ex8fmodule
 22:   implicit none

 24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 25: !                   Variable declarations
 26: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 27: !
 28: !  Variables:
 29: !     ksp     - linear solver context
 30: !     x, b, u  - approx solution, right-hand side, exact solution vectors
 31: !     A        - matrix that defines linear system
 32: !     its      - iterations for convergence
 33: !     norm     - norm of error in solution
 34: !     rctx     - random number context
 35: !

 37:   Mat A
 38:   Vec x, b, u
 39:   PC pc
 40:   PetscInt n, dim, istart, iend
 41:   PetscInt i, j, jj, ii, one, zero
 42:   PetscErrorCode ierr
 43:   PetscScalar v
 44:   PetscScalar pfive
 45:   KSP ksp

 47: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 48: !                 Beginning of program
 49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 51:   PetscCallA(PetscInitialize(ierr))
 52:   pfive = .5
 53:   n = 6
 54:   dim = n*n
 55:   one = 1
 56:   zero = 0

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

 63: !  Create parallel matrix, specifying only its global dimensions.
 64: !  When using MatCreate(), the matrix format can be specified at
 65: !  runtime. Also, the parallel partitioning of the matrix is
 66: !  determined by PETSc at runtime.

 68:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 69:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, dim, dim, ierr))
 70:   PetscCallA(MatSetFromOptions(A, ierr))
 71:   PetscCallA(MatSetUp(A, ierr))

 73: !  Currently, all PETSc parallel matrix formats are partitioned by
 74: !  contiguous chunks of rows across the processors.  Determine which
 75: !  rows of the matrix are locally owned.

 77:   PetscCallA(MatGetOwnershipRange(A, Istart, Iend, ierr))

 79: !  Set matrix elements in parallel.
 80: !   - Each processor needs to insert only elements that it owns
 81: !     locally (but any non-local elements will be sent to the
 82: !     appropriate processor during matrix assembly).
 83: !   - Always specify global rows and columns of matrix entries.

 85:   do II = Istart, Iend - 1
 86:     v = -1.0
 87:     i = II/n
 88:     j = II - i*n
 89:     if (i > 0) then
 90:       JJ = II - n
 91:       PetscCallA(MatSetValues(A, one, [II], one, [JJ], [v], ADD_VALUES, ierr))
 92:     end if
 93:     if (i < n - 1) then
 94:       JJ = II + n
 95:       PetscCallA(MatSetValues(A, one, [II], one, [JJ], [v], ADD_VALUES, ierr))
 96:     end if
 97:     if (j > 0) then
 98:       JJ = II - 1
 99:       PetscCallA(MatSetValues(A, one, [II], one, [JJ], [v], ADD_VALUES, ierr))
100:     end if
101:     if (j < n - 1) then
102:       JJ = II + 1
103:       PetscCallA(MatSetValues(A, one, [II], one, [JJ], [v], ADD_VALUES, ierr))
104:     end if
105:     v = 4.0
106:     PetscCallA(MatSetValues(A, one, [II], one, [II], [v], ADD_VALUES, ierr))
107:   end do

109: !  Assemble matrix, using the 2-step process:
110: !       MatAssemblyBegin(), MatAssemblyEnd()
111: !  Computations can be done while messages are in transition
112: !  by placing code between these two statements.

114:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
115:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

117: !  Create parallel vectors.
118: !   - Here, the parallel partitioning of the vector is determined by
119: !     PETSc at runtime.  We could also specify the local dimensions
120: !     if desired.
121: !   - Note: We form 1 vector from scratch and then duplicate as needed.

123:   PetscCallA(VecCreate(PETSC_COMM_WORLD, u, ierr))
124:   PetscCallA(VecSetSizes(u, PETSC_DECIDE, dim, ierr))
125:   PetscCallA(VecSetFromOptions(u, ierr))
126:   PetscCallA(VecDuplicate(u, b, ierr))
127:   PetscCallA(VecDuplicate(b, x, ierr))

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

131:   PetscCallA(VecSet(u, pfive, ierr))
132:   PetscCallA(MatMult(A, u, b, ierr))

134: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
135: !         Create the linear solver and set various options
136: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

138: !  Create linear solver context

140:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))
141:   PetscCallA(KSPGetPC(ksp, pc, ierr))
142:   PetscCallA(PCSetType(pc, PCMG, ierr))
143:   PetscCallA(PCMGSetLevels(pc, one, PETSC_NULL_MPI_COMM, ierr))
144:   PetscCallA(PCMGSetResidual(pc, zero, MyResidual, A, ierr))

146: !  Set operators. Here the matrix that defines the linear system
147: !  also serves as the matrix used to construct the preconditioner.

149:   PetscCallA(KSPSetOperators(ksp, A, A, ierr))

151:   PetscCallA(KSPDestroy(ksp, ierr))
152:   PetscCallA(VecDestroy(u, ierr))
153:   PetscCallA(VecDestroy(x, ierr))
154:   PetscCallA(VecDestroy(b, ierr))
155:   PetscCallA(MatDestroy(A, ierr))

157:   PetscCallA(PetscFinalize(ierr))
158: end

160: !/*TEST
161: !
162: !   test:
163: !      nsize: 1
164: !      output_file: output/empty.out
165: !TEST*/