Actual source code: ex5f.F90

```  1:       !Solves two linear systems in parallel with KSP.  The code
2:       !illustrates repeated solution of linear systems with the same preconditioner
3:       !method but different matrices (having the same nonzero structure).  The code
4:       !also uses multiple profiling stages.  Input arguments are
5:       !  -m <size> : problem size
6:       !  -mat_nonsym : use nonsymmetric matrix (default is symmetric)

8:       !Concepts: KSP^repeatedly solving linear systems;
9:       !Concepts: PetscLog^profiling multiple stages of code;
10:       !Processors: n

12: program main
13: #include <petsc/finclude/petscksp.h>
14:       use petscksp

16:       implicit none
17:       KSP            :: ksp            ! linear solver context
18:       Mat            :: C,Ctmp           ! matrix
19:       Vec            :: x,u,b            ! approx solution, RHS, exact solution
20:       PetscReal      :: norm             ! norm of solution error
21:       PetscScalar    :: v
22:       PetscScalar, parameter :: myNone = -1.0
23:       PetscInt       :: Ii,JJ,ldim,low,high,iglobal,Istart,Iend
24:       PetscErrorCode :: ierr
25:       PetscInt       :: i,j,its,n
26:       PetscInt       :: m = 3, orthog = 0
27:       PetscMPIInt    :: size,rank
28:       PetscBool :: &
29:         testnewC         = PETSC_FALSE, &
30:         testscaledMat    = PETSC_FALSE, &
31:         mat_nonsymmetric = PETSC_FALSE
32:       PetscBool      :: flg
33:       PetscRandom    :: rctx
34:       PetscLogStage,dimension(0:1) :: stages
35:       character(len=PETSC_MAX_PATH_LEN) :: outputString
36:       PetscInt,parameter :: one = 1

38:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
39:       if (ierr /= 0) then
40:         write(6,*)'Unable to initialize PETSc'
41:         stop
42:       endif

44:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-orthog',orthog,flg,ierr)
45:       CHKERRA(ierr)
46:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
47:       CHKERRA(ierr)
48:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
49:       CHKERRA(ierr)
50:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
51:       CHKERRA(ierr)
52:       n=2*size

54:       ! Set flag if we are doing a nonsymmetric problem; the default is symmetric.

56:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-mat_nonsym",mat_nonsymmetric,flg,ierr)
57:       CHKERRA(ierr)
58:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_scaledMat",testscaledMat,flg,ierr)
59:       CHKERRA(ierr)

61:       ! Register two stages for separate profiling of the two linear solves.
62:       ! Use the runtime option -log_view for a printout of performance
63:       ! statistics at the program's conlusion.

65:       call PetscLogStageRegister("Original Solve",stages(0),ierr)
66:       CHKERRA(ierr)
67:       call PetscLogStageRegister("Second Solve",stages(1),ierr)
68:       CHKERRA(ierr)

70:       ! -------------- Stage 0: Solve Original System ----------------------
71:       ! Indicate to PETSc profiling that we're beginning the first stage

73:       call PetscLogStagePush(stages(0),ierr)
74:       CHKERRA(ierr)

76:       ! Create parallel matrix, specifying only its global dimensions.
77:       ! When using MatCreate(), the matrix format can be specified at
78:       ! runtime. Also, the parallel partitioning of the matrix is
79:       ! determined by PETSc at runtime.

81:       call MatCreate(PETSC_COMM_WORLD,C,ierr)
82:       CHKERRA(ierr)
83:       call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
84:       CHKERRA(ierr)
85:       call MatSetFromOptions(C,ierr)
86:       CHKERRA(ierr)
87:       call MatSetUp(C,ierr)
88:       CHKERRA(ierr)

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

94:       call MatGetOwnershipRange(C,Istart,Iend,ierr)

96:       ! Set matrix entries matrix in parallel.
97:       ! - Each processor needs to insert only elements that it owns
98:       ! locally (but any non-local elements will be sent to the
99:       ! appropriate processor during matrix assembly).
100:       !- Always specify global row and columns of matrix entries.

102:       intitializeC: do Ii=Istart,Iend-1
103:           v =-1.0; i = Ii/n; j = Ii - i*n
104:           if (i>0) then
105:             JJ = Ii - n
107:             CHKERRA(ierr)
108:           endif

110:           if (i<m-1) then
111:             JJ = Ii + n
113:             CHKERRA(ierr)
114:           endif

116:           if (j>0) then
117:             JJ = Ii - 1
119:             CHKERRA(ierr)
120:           endif

122:           if (j<n-1) then
123:             JJ = Ii + 1
125:             CHKERRA(ierr)
126:           endif

128:           v=4.0
130:           CHKERRA(ierr)

132:       enddo intitializeC

134:       ! Make the matrix nonsymmetric if desired
135:       if (mat_nonsymmetric) then
136:         do Ii=Istart,Iend-1
137:           v=-1.5; i=Ii/n
138:           if (i>1) then
139:             JJ=Ii-n-1
141:             CHKERRA(ierr)
142:           endif
143:         enddo
144:       else
145:         call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
146:         CHKERRA(ierr)
147:         call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
148:         CHKERRA(ierr)
149:       endif

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

156:       call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
157:       CHKERRA(ierr)
158:       call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
159:       CHKERRA(ierr)

161:       ! Create parallel vectors.
162:       ! - When using VecSetSizes(), we specify only the vector's global
163:       !   dimension; the parallel partitioning is determined at runtime.
164:       ! - Note: We form 1 vector from scratch and then duplicate as needed.

166:       call  VecCreate(PETSC_COMM_WORLD,u,ierr)
167:       call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
168:       call  VecSetFromOptions(u,ierr)
169:       call  VecDuplicate(u,b,ierr)
170:       call  VecDuplicate(b,x,ierr)

172:       ! Currently, all parallel PETSc vectors are partitioned by
173:       ! contiguous chunks across the processors.  Determine which
174:       ! range of entries are locally owned.

176:       call  VecGetOwnershipRange(x,low,high,ierr)
177:       CHKERRA(ierr)

179:       !Set elements within the exact solution vector in parallel.
180:       ! - Each processor needs to insert only elements that it owns
181:       ! locally (but any non-local entries will be sent to the
182:       ! appropriate processor during vector assembly).
183:       ! - Always specify global locations of vector entries.

185:       call VecGetLocalSize(x,ldim,ierr)
186:       CHKERRA(ierr)
187:       do i=0,ldim-1
188:         iglobal = i + low
189:         v = real(i + 100*rank)
190:         call VecSetValues(u,one,iglobal,v,INSERT_VALUES,ierr)
191:         CHKERRA(ierr)
192:       enddo

194:       ! Assemble vector, using the 2-step process:
195:       ! VecAssemblyBegin(), VecAssemblyEnd()
196:       ! Computations can be done while messages are in transition,
197:       ! by placing code between these two statements.

199:       call  VecAssemblyBegin(u,ierr)
200:       CHKERRA(ierr)
201:       call  VecAssemblyEnd(u,ierr)
202:       CHKERRA(ierr)

204:       ! Compute right-hand-side vector

206:       call  MatMult(C,u,b,ierr)

208:       CHKERRA(ierr)

210:       ! Create linear solver context

212:       call  KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
213:       CHKERRA(ierr)
214:       ! Set operators. Here the matrix that defines the linear system
215:       ! also serves as the preconditioning matrix.

217:       call  KSPSetOperators(ksp,C,C,ierr)
218:       CHKERRA(ierr)
219:       ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)

221:       call  KSPSetFromOptions(ksp,ierr)
222:       CHKERRA(ierr)
223:       ! Solve linear system.  Here we explicitly call KSPSetUp() for more
224:       ! detailed performance monitoring of certain preconditioners, such
225:       ! as ICC and ILU.  This call is optional, as KSPSetUp() will
226:       ! automatically be called within KSPSolve() if it hasn't been

229:       call  KSPSetUp(ksp,ierr)
230:       CHKERRA(ierr)

232:       ! Do not do this in application code, use -ksp_gmres_modifiedgramschmidt or -ksp_gmres_modifiedgramschmidt
233:       if (orthog .eq. 1) then
234:          call KSPGMRESSetOrthogonalization(ksp,KSPGMRESModifiedGramSchmidtOrthogonalization,ierr)
235:       else if (orthog .eq. 2) then
236:          call KSPGMRESSetOrthogonalization(ksp,KSPGMRESClassicalGramSchmidtOrthogonalization,ierr)
237:       endif
238:       CHKERRA(ierr)

240:       call  KSPSolve(ksp,b,x,ierr)
241:       CHKERRA(ierr)

243:       ! Check the error

245:       call VecAXPY(x,myNone,u,ierr)
246:       call VecNorm(x,NORM_2,norm,ierr)

248:       call KSPGetIterationNumber(ksp,its,ierr)
249:       if (.not. testscaledMat .or. norm > 1.e-7) then
250:         write(outputString,'(a,f11.9,a,i2.2,a)') 'Norm of error ',norm,', Iterations ',its,'\n'
251:         call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
252:       endif

254:       ! -------------- Stage 1: Solve Second System ----------------------

256:       ! Solve another linear system with the same method.  We reuse the KSP
257:       ! context, matrix and vector data structures, and hence save the
258:       ! overhead of creating new ones.

260:       ! Indicate to PETSc profiling that we're concluding the first
261:       ! stage with PetscLogStagePop(), and beginning the second stage with
262:       ! PetscLogStagePush().

264:       call PetscLogStagePop(ierr)
265:       CHKERRA(ierr)
266:       call PetscLogStagePush(stages(1),ierr)
267:       CHKERRA(ierr)

269:       ! Initialize all matrix entries to zero.  MatZeroEntries() retains the
270:       ! nonzero structure of the matrix for sparse formats.

272:       call MatZeroEntries(C,ierr)
273:       CHKERRA(ierr)

275:       ! Assemble matrix again.  Note that we retain the same matrix data
276:       ! structure and the same nonzero pattern; we just change the values
277:       ! of the matrix entries.

279:       do i=0,m-1
280:         do j=2*rank,2*rank+1
281:           v =-1.0; Ii=j + n*i
282:           if (i>0) then
283:             JJ = Ii - n
285:             CHKERRA(ierr)
286:           endif

288:           if (i<m-1) then
289:             JJ = Ii + n
291:             CHKERRA(ierr)
292:           endif

294:           if (j>0) then
295:             JJ = Ii - 1
297:             CHKERRA(ierr)
298:           endif

300:           if (j<n-1) then
301:             JJ = Ii + 1
303:             CHKERRA(ierr)
304:           endif

306:           v=6.0
308:           CHKERRA(ierr)

310:         enddo
311:       enddo

313:       ! Make the matrix nonsymmetric if desired

315:       if (mat_nonsymmetric) then
316:         do Ii=Istart,Iend-1
317:           v=-1.5;  i=Ii/n
318:           if (i>1) then
319:             JJ=Ii-n-1
321:             CHKERRA(ierr)
322:           endif
323:         enddo
324:       endif

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

331:       call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
332:       CHKERRA(ierr)
333:       call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
334:       CHKERRA(ierr)

336:       if (testscaledMat) then
337:         ! Scale a(0,0) and a(M-1,M-1)

339:         if (rank /= 0) then
340:           v = 6.0*0.00001; Ii = 0; JJ = 0
341:           call MatSetValues(C,one,Ii,one,JJ,v,INSERT_VALUES,ierr)
342:           CHKERRA(ierr)
343:         elseif (rank == size -1) then
344:           v = 6.0*0.00001; Ii = m*n-1; JJ = m*n-1
345:           call MatSetValues(C,one,Ii,one,JJ,v,INSERT_VALUES,ierr)

347:         endif

349:         call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
350:         call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)

352:         ! Compute a new right-hand-side vector

354:         call  VecDestroy(u,ierr)
355:         call  VecCreate(PETSC_COMM_WORLD,u,ierr)
356:         call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
357:         call  VecSetFromOptions(u,ierr)

359:         call  PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
360:         call  PetscRandomSetFromOptions(rctx,ierr)
361:         call  VecSetRandom(u,rctx,ierr)
362:         call  PetscRandomDestroy(rctx,ierr)
363:         call  VecAssemblyBegin(u,ierr)
364:         call  VecAssemblyEnd(u,ierr)

366:       endif

368:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_newMat",testnewC,flg,ierr)
369:       CHKERRA(ierr)

371:       if (testnewC) then
372:       ! User may use a new matrix C with same nonzero pattern, e.g.
373:       ! ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat

375:         call  MatDuplicate(C,MAT_COPY_VALUES,Ctmp,ierr)
376:         call  MatDestroy(C,ierr)
377:         call  MatDuplicate(Ctmp,MAT_COPY_VALUES,C,ierr)
378:         call  MatDestroy(Ctmp,ierr)
379:       endif

381:       call MatMult(C,u,b,ierr);CHKERRA(ierr)

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

386:       call KSPSetOperators(ksp,C,C,ierr);CHKERRA(ierr)

388:       ! Solve linear system

390:       call  KSPSetUp(ksp,ierr); CHKERRA(ierr)
391:       call  KSPSolve(ksp,b,x,ierr); CHKERRA(ierr)

393:       ! Check the error

395:       call VecAXPY(x,myNone,u,ierr); CHKERRA(ierr)
396:       call VecNorm(x,NORM_2,norm,ierr); CHKERRA(ierr)
397:       call KSPGetIterationNumber(ksp,its,ierr); CHKERRA(ierr)
398:       if (.not. testscaledMat .or. norm > 1.e-7) then
399:         write(outputString,'(a,f11.9,a,i2.2,a)') 'Norm of error ',norm,', Iterations ',its,'\n'
400:         call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
401:       endif

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

406:       call  KSPDestroy(ksp,ierr); CHKERRA(ierr)
407:       call  VecDestroy(u,ierr); CHKERRA(ierr)
408:       call  VecDestroy(x,ierr); CHKERRA(ierr)
409:       call  VecDestroy(b,ierr); CHKERRA(ierr)
410:       call  MatDestroy(C,ierr); CHKERRA(ierr)

412:       ! Indicate to PETSc profiling that we're concluding the second stage

414:       call PetscLogStagePop(ierr)
415:       CHKERRA(ierr)

417:       call PetscFinalize(ierr)

419: !/*TEST
420: !
421: !   test:
422: !      args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always
423: !
424: !   test:
425: !      suffix: 2
426: !      nsize: 2
427: !      args: -pc_type jacobi -ksp_monitor_short -ksp_gmres_cgs_refinement_type refine_always -ksp_rtol .000001
428: !
429: !   test:
430: !      suffix: 5
431: !      nsize: 2
432: !      args: -ksp_gmres_cgs_refinement_type refine_always -ksp_monitor draw::draw_lg -ksp_monitor_true_residual draw::draw_lg
433: !      output_file: output/ex5f_5.out
434: !
435: !   test:
436: !      suffix: asm
437: !      nsize: 4
438: !      args: -pc_type asm
439: !      output_file: output/ex5f_asm.out
440: !
441: !   test:
442: !      suffix: asm_baij
443: !      nsize: 4
444: !      args: -pc_type asm -mat_type baij
445: !      output_file: output/ex5f_asm.out
446: !
447: !   test:
448: !      suffix: redundant_0
449: !      args: -m 1000 -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
450: !
451: !   test:
452: !      suffix: redundant_1
453: !      nsize: 5
454: !      args: -pc_type redundant -pc_redundant_number 1 -redundant_ksp_type gmres -redundant_pc_type jacobi
455: !
456: !   test:
457: !      suffix: redundant_2
458: !      nsize: 5
459: !      args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi
460: !
461: !   test:
462: !      suffix: redundant_3
463: !      nsize: 5
464: !      args: -pc_type redundant -pc_redundant_number 5 -redundant_ksp_type gmres -redundant_pc_type jacobi
465: !
466: !   test:
467: !      suffix: redundant_4
468: !      nsize: 5
469: !      args: -pc_type redundant -pc_redundant_number 3 -redundant_ksp_type gmres -redundant_pc_type jacobi -psubcomm_type interlaced
470: !
471: !   test:
472: !      suffix: superlu_dist
473: !      nsize: 15
474: !      requires: superlu_dist
475: !      args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat
476: !
477: !   test:
478: !      suffix: superlu_dist_2
479: !      nsize: 15
480: !      requires: superlu_dist
481: !      args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 150 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact SamePattern_SameRowPerm
482: !      output_file: output/ex5f_superlu_dist.out
483: !
484: !   test:
485: !      suffix: superlu_dist_3
486: !      nsize: 15
487: !      requires: superlu_dist
488: !      args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_equil false -m 500 -mat_superlu_dist_r 3 -mat_superlu_dist_c 5 -test_scaledMat -mat_superlu_dist_fact DOFACT
489: !      output_file: output/ex5f_superlu_dist.out
490: !
491: !   test:
492: !      suffix: superlu_dist_0
493: !      nsize: 1
494: !      requires: superlu_dist
495: !      args: -pc_type lu -pc_factor_mat_solver_type superlu_dist -test_scaledMat
496: !      output_file: output/ex5f_superlu_dist.out
497: !
498: !   test:
499: !      suffix: orthog1
500: !      args: -orthog 1 -ksp_view
501: !
502: !   test:
503: !      suffix: orthog2
504: !      args: -orthog 2 -ksp_view
505: !
506: !TEST*/

508: end program main
```