1: !
  2: !   Description: Tests PCFieldSplitGetIS and PCFieldSplitSetIS from Fortran.
  3: !
  4: #include <petsc/finclude/petscksp.h>
  5: program main
  6:   use petscksp
  7:   implicit none

  9:   Vec x, b, u
 10:   Mat A
 11:   KSP ksp
 12:   PC pc
 13:   PetscReal norm
 14:   PetscErrorCode ierr
 15:   PetscInt i, n, col(3), its
 16:   PetscInt :: nksp
 17:   PetscBool flg
 18:   PetscMPIInt size
 19:   PetscScalar, parameter :: one = 1.0, none = -1.0
 20:   PetscScalar value(3)
 21:   KSP, pointer :: subksp(:)

 23:   IS isin, isout

 25: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 26: !                 Beginning of program
 27: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 29:   PetscCallA(PetscInitialize(ierr))
 30:   PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD, size, ierr))
 31:   PetscCheckA(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, 'This is a uniprocessor example only')
 32:   n = 10
 33:   PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS, PETSC_NULL_CHARACTER, '-n', n, flg, ierr))

 35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 36: !         Compute the matrix and right-hand-side vector that define
 37: !         the linear system, Ax = b.
 38: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 40: !  Create matrix.  When using MatCreate(), the matrix format can
 41: !  be specified at runtime.

 43:   PetscCallA(MatCreate(PETSC_COMM_WORLD, A, ierr))
 44:   PetscCallA(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, n, n, ierr))
 45:   PetscCallA(MatSetFromOptions(A, ierr))
 46:   PetscCallA(MatSetUp(A, ierr))

 48: !  Assemble matrix.
 49: !   - Note that MatSetValues() uses 0-based row and column numbers
 50: !     in Fortran as well as in C (as set here in the array "col").

 52:   value = [-1.0, 2.0, -1.0]
 53:   do i = 1, n - 2
 54:     col(1) = i - 1
 55:     col(2) = i
 56:     col(3) = i + 1
 57:     PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [i], 3_PETSC_INT_KIND, col, value, INSERT_VALUES, ierr))
 58:   end do
 59:   i = n - 1
 60:   col(1) = n - 2
 61:   col(2) = n - 1
 62:   PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, value, INSERT_VALUES, ierr))
 63:   i = 0
 64:   col(1) = 0
 65:   col(2) = 1
 66:   value(1) = 2.0
 67:   value(2) = -1.0
 68:   PetscCallA(MatSetValues(A, 1_PETSC_INT_KIND, [i], 2_PETSC_INT_KIND, col, value, INSERT_VALUES, ierr))
 69:   PetscCallA(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY, ierr))
 70:   PetscCallA(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY, ierr))

 72: !  Create vectors.  Note that we form 1 vector from scratch and
 73: !  then duplicate as needed.

 75:   PetscCallA(VecCreate(PETSC_COMM_WORLD, x, ierr))
 76:   PetscCallA(VecSetSizes(x, PETSC_DECIDE, n, ierr))
 77:   PetscCallA(VecSetFromOptions(x, ierr))
 78:   PetscCallA(VecDuplicate(x, b, ierr))
 79:   PetscCallA(VecDuplicate(x, u, ierr))

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

 83:   PetscCallA(VecSet(u, one, ierr))
 84:   PetscCallA(MatMult(A, u, b, ierr))

 86: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 87: !          Create the linear solver and set various options
 88: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 90: !  Create linear solver context

 92:   PetscCallA(KSPCreate(PETSC_COMM_WORLD, ksp, ierr))

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

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

 99: !  Set linear solver defaults for this problem (optional).
100: !   - By extracting the KSP and PC contexts from the KSP context,
101: !     we can then directly call any KSP and PC routines
102: !     to set various options.
103: !   - The following four statements are optional; all of these
104: !     parameters could alternatively be specified at runtime via
105: !     KSPSetFromOptions()

107:   PetscCallA(KSPGetPC(ksp, pc, ierr))
108:   PetscCallA(PCSetType(pc, PCFIELDSPLIT, ierr))
109:   PetscCallA(ISCreateStride(PETSC_COMM_SELF, n, 0_PETSC_INT_KIND, 1_PETSC_INT_KIND, isin, ierr))
110:   PetscCallA(PCFieldSplitSetIS(pc, 'splitname', isin, ierr))
111:   PetscCallA(PCFieldSplitGetIS(pc, 'splitname', isout, ierr))
112:   PetscCheckA(isin == isout, PETSC_COMM_SELF, PETSC_ERR_PLIB, 'PCFieldSplitGetIS() failed')

114: !  Set runtime options, e.g.,
115: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
116: !  These options will override those specified above as long as
117: !  KSPSetFromOptions() is called _after_ any other customization
118: !  routines.

120:   PetscCallA(KSPSetFromOptions(ksp, ierr))

122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: !                      Solve the linear system
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:   PetscCallA(PCSetUp(pc, ierr))
126:   PetscCallA(PCFieldSplitGetSubKSP(pc, nksp, subksp, ierr))
127:   PetscCheckA(nksp == 2, PETSC_COMM_WORLD, PETSC_ERR_PLIB, 'Number of KSP should be two')
128:   PetscCallA(KSPView(subksp(1), PETSC_VIEWER_STDOUT_WORLD, ierr))
129:   PetscCallA(PCFieldSplitRestoreSubKSP(pc, nksp, subksp, ierr))

131:   PetscCallA(KSPSolve(ksp, b, x, ierr))

133: !  View solver info; we could instead use the option -ksp_view

135:   PetscCallA(KSPView(ksp, PETSC_VIEWER_STDOUT_WORLD, ierr))

137: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: !                      Check solution and clean up
139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

141: !  Check the error

143:   PetscCallA(VecAXPY(x, none, u, ierr))
144:   PetscCallA(VecNorm(x, NORM_2, norm, ierr))
145:   PetscCallA(KSPGetIterationNumber(ksp, its, ierr))
146:   if (norm > 1.e-12) then
147:     write (6, 100) norm, its
148:   else
149:     write (6, 200) its
150:   end if
151: 100 format('Norm of error ', e11.4, ',  Iterations = ', i5)
152: 200 format('Norm of error < 1.e-12, Iterations = ', i5)

154: !  Free work space.  All PETSc objects should be destroyed when they
155: !  are no longer needed.

157:   PetscCallA(ISDestroy(isin, ierr))
158:   PetscCallA(VecDestroy(x, ierr))
159:   PetscCallA(VecDestroy(u, ierr))
160:   PetscCallA(VecDestroy(b, ierr))
161:   PetscCallA(MatDestroy(A, ierr))
162:   PetscCallA(KSPDestroy(ksp, ierr))
163:   PetscCallA(PetscFinalize(ierr))

165: end

167: !/*TEST
168: !
169: !     test:
170: !       args: -ksp_monitor
171: !
172: !TEST*/