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*/