Actual source code: ex8f.F90
1: !
2: ! Tests PCMGSetResidual
3: !
4: ! -----------------------------------------------------------------------
6: program main
7: #include <petsc/finclude/petscksp.h>
8: use petscksp
9: implicit none
11: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
12: ! Variable declarations
13: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
14: !
15: ! Variables:
16: ! ksp - linear solver context
17: ! x, b, u - approx solution, right-hand-side, exact solution vectors
18: ! A - matrix that defines linear system
19: ! its - iterations for convergence
20: ! norm - norm of error in solution
21: ! rctx - random number context
22: !
24: Mat A
25: Vec x,b,u
26: PC pc
27: PetscInt n,dim,istart,iend
28: PetscInt i,j,jj,ii,one,zero
29: PetscErrorCode ierr
30: PetscScalar v
31: external MyResidual
32: PetscScalar pfive
33: KSP ksp
35: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
36: ! Beginning of program
37: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
39: PetscCallA(PetscInitialize(ierr))
40: pfive = .5
41: n = 6
42: dim = n*n
43: one = 1
44: zero = 0
46: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
47: ! Compute the matrix and right-hand-side vector that define
48: ! the linear system, Ax = b.
49: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
51: ! Create parallel matrix, specifying only its global dimensions.
52: ! When using MatCreate(), the matrix format can be specified at
53: ! runtime. Also, the parallel partitioning of the matrix is
54: ! determined by PETSc at runtime.
56: PetscCallA(MatCreate(PETSC_COMM_WORLD,A,ierr))
57: PetscCallA(MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,dim,dim,ierr))
58: PetscCallA(MatSetFromOptions(A,ierr))
59: PetscCallA(MatSetUp(A,ierr))
61: ! Currently, all PETSc parallel matrix formats are partitioned by
62: ! contiguous chunks of rows across the processors. Determine which
63: ! rows of the matrix are locally owned.
65: PetscCallA(MatGetOwnershipRange(A,Istart,Iend,ierr))
67: ! Set matrix elements in parallel.
68: ! - Each processor needs to insert only elements that it owns
69: ! locally (but any non-local elements will be sent to the
70: ! appropriate processor during matrix assembly).
71: ! - Always specify global rows and columns of matrix entries.
73: do 10, II=Istart,Iend-1
74: v = -1.0
75: i = II/n
76: j = II - i*n
77: if (i.gt.0) then
78: JJ = II - n
79: PetscCallA(MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr))
80: endif
81: if (i.lt.n-1) then
82: JJ = II + n
83: PetscCallA(MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr))
84: endif
85: if (j.gt.0) then
86: JJ = II - 1
87: PetscCallA(MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr))
88: endif
89: if (j.lt.n-1) then
90: JJ = II + 1
91: PetscCallA(MatSetValues(A,one,II,one,JJ,v,ADD_VALUES,ierr))
92: endif
93: v = 4.0
94: PetscCallA( MatSetValues(A,one,II,one,II,v,ADD_VALUES,ierr))
95: 10 continue
97: ! Assemble matrix, using the 2-step process:
98: ! MatAssemblyBegin(), MatAssemblyEnd()
99: ! Computations can be done while messages are in transition
100: ! by placing code between these two statements.
102: PetscCallA(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY,ierr))
103: PetscCallA(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY,ierr))
105: ! Create parallel vectors.
106: ! - Here, the parallel partitioning of the vector is determined by
107: ! PETSc at runtime. We could also specify the local dimensions
108: ! if desired.
109: ! - Note: We form 1 vector from scratch and then duplicate as needed.
111: PetscCallA(VecCreate(PETSC_COMM_WORLD,u,ierr))
112: PetscCallA(VecSetSizes(u,PETSC_DECIDE,dim,ierr))
113: PetscCallA(VecSetFromOptions(u,ierr))
114: PetscCallA(VecDuplicate(u,b,ierr))
115: PetscCallA(VecDuplicate(b,x,ierr))
117: ! Set exact solution; then compute right-hand-side vector.
119: PetscCallA(VecSet(u,pfive,ierr))
120: PetscCallA(MatMult(A,u,b,ierr))
122: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: ! Create the linear solver and set various options
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: ! Create linear solver context
128: PetscCallA(KSPCreate(PETSC_COMM_WORLD,ksp,ierr))
129: PetscCallA(KSPGetPC(ksp,pc,ierr))
130: PetscCallA(PCSetType(pc,PCMG,ierr))
131: PetscCallA(PCMGSetLevels(pc,one,PETSC_NULL_MPI_COMM,ierr))
132: PetscCallA(PCMGSetResidual(pc,zero,MyResidual,A,ierr))
134: ! Set operators. Here the matrix that defines the linear system
135: ! also serves as the preconditioning matrix.
137: PetscCallA(KSPSetOperators(ksp,A,A,ierr))
139: PetscCallA(KSPDestroy(ksp,ierr))
140: PetscCallA(VecDestroy(u,ierr))
141: PetscCallA(VecDestroy(x,ierr))
142: PetscCallA(VecDestroy(b,ierr))
143: PetscCallA(MatDestroy(A,ierr))
145: PetscCallA(PetscFinalize(ierr))
146: end
148: subroutine MyResidual(A,b,x,r,ierr)
149: use petscksp
150: implicit none
151: Mat A
152: Vec b,x,r
153: integer ierr
154: return
155: end
157: !/*TEST
158: !
159: ! test:
160: ! nsize: 1
161: !TEST*/