Actual source code: ex5f.F90

petsc-3.14.2 2020-12-03
Report Typos and Errors
  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


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

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

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

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

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

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

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

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

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

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

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

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

111:           if (i<m-1) then
112:             JJ = Ii + n
113:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
114:             CHKERRA(ierr)
115:           endif

117:           if (j>0) then
118:             JJ = Ii - 1
119:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
120:             CHKERRA(ierr)
121:           endif

123:           if (j<n-1) then
124:             JJ = Ii + 1
125:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
126:             CHKERRA(ierr)
127:           endif

129:           v=4.0
130:           call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
131:           CHKERRA(ierr)

133:       enddo intitializeC

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

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

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

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

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

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

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

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

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

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

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

205:       ! Compute right-hand-side vector

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

209:       CHKERRA(ierr)

211:       ! Create linear solver context

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

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

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

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

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

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

244:       ! Check the error

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

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

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

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

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

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

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

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

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

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

289:           if (i<m-1) then
290:             JJ = Ii + n
291:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
292:             CHKERRA(ierr)
293:           endif

295:           if (j>0) then
296:             JJ = Ii - 1
297:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
298:             CHKERRA(ierr)
299:           endif

301:           if (j<n-1) then
302:             JJ = Ii + 1
303:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
304:             CHKERRA(ierr)
305:           endif

307:           v=6.0
308:           call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
309:           CHKERRA(ierr)

311:         enddo
312:       enddo

314:       ! Make the matrix nonsymmetric if desired

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

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

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

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

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

348:         endif

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

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

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

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

367:       endif

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

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

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

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

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

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

389:       ! Solve linear system

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

394:       ! Check the error

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

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

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

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

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

418:       call PetscFinalize(ierr)

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

509: end program main