Actual source code: ex9f.F90

  1: !
  2: !   Description: Tests PCFieldSplitGetIS and PCFieldSplitSetIS from Fortran.
  3: !
  4:       program main
  5: #include <petsc/finclude/petscksp.h>
  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,i1,i2,i3
 16:       PetscInt ione,izero
 17:       PetscBool  flg
 18:       PetscMPIInt size
 19:       PetscScalar      none,one,value(3)
 20:       IS isin,isout

 22: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 23: !                 Beginning of program
 24: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 26:       PetscCallA(PetscInitialize(ierr))
 27:       PetscCallMPIA(MPI_Comm_size(PETSC_COMM_WORLD,size,ierr))
 28:       if (size .ne. 1) then; SETERRA(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,'This is a uniprocessor example only'); endif
 29:       none = -1.0
 30:       one  = 1.0
 31:       n    = 10
 32:       i1 = 1
 33:       i2 = 2
 34:       i3 = 3
 35:       PetscCallA(PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-n',n,flg,ierr))

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

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

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

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

 54:       value(1) = -1.0
 55:       value(2) = 2.0
 56:       value(3) = -1.0
 57:       do 50 i=1,n-2
 58:          col(1) = i-1
 59:          col(2) = i
 60:          col(3) = i+1
 61:          PetscCallA(MatSetValues(A,i1,i,i3,col,value,INSERT_VALUES,ierr))
 62:   50  continue
 63:       i = n - 1
 64:       col(1) = n - 2
 65:       col(2) = n - 1
 66:       PetscCallA(MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr))
 67:       i = 0
 68:       col(1) = 0
 69:       col(2) = 1
 70:       value(1) = 2.0
 71:       value(2) = -1.0
 72:       PetscCallA(MatSetValues(A,i1,i,i2,col,value,INSERT_VALUES,ierr))
 73:       PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
 74:       PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))

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

 79:       PetscCallA(VecCreate(PETSC_COMM_WORLD,x,ierr))
 80:       PetscCallA(VecSetSizes(x,PETSC_DECIDE,n,ierr))
 81:       PetscCallA(VecSetFromOptions(x,ierr))
 82:       PetscCallA(VecDuplicate(x,b,ierr))
 83:       PetscCallA(VecDuplicate(x,u,ierr))

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

 87:       PetscCallA(VecSet(u,one,ierr))
 88:       PetscCallA(MatMult(A,u,b,ierr))

 90: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 91: !          Create the linear solver and set various options
 92: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 94: !  Create linear solver context

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

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

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

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

111:       PetscCallA(KSPGetPC(ksp,pc,ierr))
112:       PetscCallA(PCSetType(pc,PCFIELDSPLIT,ierr))
113:       izero = 0
114:       ione  = 1
115:       PetscCallA(ISCreateStride(PETSC_COMM_SELF,n,izero,ione,isin,ierr))
116:       PetscCallA(PCFieldSplitSetIS(pc,"splitname",isin,ierr))
117:       PetscCallA(PCFieldSplitGetIS(pc,"splitname",isout,ierr))
118:       if (isin .ne. isout) then ; SETERRA(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCFieldSplitGetIS() failed"); endif

120: !  Set runtime options, e.g.,
121: !      -ksp_type <type> -pc_type <type> -ksp_monitor -ksp_rtol <rtol>
122: !  These options will override those specified above as long as
123: !  KSPSetFromOptions() is called _after_ any other customization
124: !  routines.

126:       PetscCallA(KSPSetFromOptions(ksp,ierr))

128: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: !                      Solve the linear system
130: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

132:       PetscCallA(KSPSolve(ksp,b,x,ierr))
133:       PetscCallA(PetscLogStagePop(ierr))

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

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

139: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: !                      Check solution and clean up
141: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

143: !  Check the error

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

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

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

167:       end

169: !/*TEST
170: !
171: !     test:
172: !       args: -ksp_monitor
173: !
174: !TEST*/