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: program main
  9: #include <petsc/finclude/petscksp.h>
 10:       use petscksp

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

 34:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 35:       if (ierr /= 0) then
 36:         write(6,*)'Unable to initialize PETSc'
 37:         stop
 38:       endif

 40:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-orthog',orthog,flg,ierr)
 41:       CHKERRA(ierr)
 42:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
 43:       CHKERRA(ierr)
 44:       call MPI_Comm_rank(PETSC_COMM_WORLD,rank,ierr)
 45:       CHKERRA(ierr)
 46:       call MPI_Comm_size(PETSC_COMM_WORLD,size,ierr)
 47:       CHKERRA(ierr)
 48:       n=2*size

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

 52:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-mat_nonsym",mat_nonsymmetric,flg,ierr)
 53:       CHKERRA(ierr)
 54:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_scaledMat",testscaledMat,flg,ierr)
 55:       CHKERRA(ierr)

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

 61:       call PetscLogStageRegister("Original Solve",stages(0),ierr)
 62:       CHKERRA(ierr)
 63:       call PetscLogStageRegister("Second Solve",stages(1),ierr)
 64:       CHKERRA(ierr)

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

 69:       call PetscLogStagePush(stages(0),ierr)
 70:       CHKERRA(ierr)

 72:       ! Create parallel matrix, specifying only its global dimensions.
 73:       ! When using MatCreate(), the matrix format can be specified at
 74:       ! runtime. Also, the parallel partitioning of the matrix is
 75:       ! determined by PETSc at runtime.

 77:       call MatCreate(PETSC_COMM_WORLD,C,ierr)
 78:       CHKERRA(ierr)
 79:       call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
 80:       CHKERRA(ierr)
 81:       call MatSetFromOptions(C,ierr)
 82:       CHKERRA(ierr)
 83:       call MatSetUp(C,ierr)
 84:       CHKERRA(ierr)

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

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

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

 98:       intitializeC: do Ii=Istart,Iend-1
 99:           v =-1.0; i = Ii/n; j = Ii - i*n
100:           if (i>0) then
101:             JJ = Ii - n
102:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
103:             CHKERRA(ierr)
104:           endif

106:           if (i<m-1) then
107:             JJ = Ii + n
108:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
109:             CHKERRA(ierr)
110:           endif

112:           if (j>0) then
113:             JJ = Ii - 1
114:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
115:             CHKERRA(ierr)
116:           endif

118:           if (j<n-1) then
119:             JJ = Ii + 1
120:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
121:             CHKERRA(ierr)
122:           endif

124:           v=4.0
125:           call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
126:           CHKERRA(ierr)

128:       enddo intitializeC

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

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

152:       call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
153:       CHKERRA(ierr)
154:       call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
155:       CHKERRA(ierr)

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

162:       call  VecCreate(PETSC_COMM_WORLD,u,ierr)
163:       call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
164:       call  VecSetFromOptions(u,ierr)
165:       call  VecDuplicate(u,b,ierr)
166:       call  VecDuplicate(b,x,ierr)

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

172:       call  VecGetOwnershipRange(x,low,high,ierr)
173:       CHKERRA(ierr)

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

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

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

195:       call  VecAssemblyBegin(u,ierr)
196:       CHKERRA(ierr)
197:       call  VecAssemblyEnd(u,ierr)
198:       CHKERRA(ierr)

200:       ! Compute right-hand-side vector

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

204:       CHKERRA(ierr)

206:       ! Create linear solver context

208:       call  KSPCreate(PETSC_COMM_WORLD,ksp,ierr)
209:       CHKERRA(ierr)
210:       ! Set operators. Here the matrix that defines the linear system
211:       ! also serves as the preconditioning matrix.

213:       call  KSPSetOperators(ksp,C,C,ierr)
214:       CHKERRA(ierr)
215:       ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)

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

225:       call  KSPSetUp(ksp,ierr)
226:       CHKERRA(ierr)

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

236:       call  KSPSolve(ksp,b,x,ierr)
237:       CHKERRA(ierr)

239:       ! Check the error

241:       call VecAXPY(x,myNone,u,ierr)
242:       call VecNorm(x,NORM_2,norm,ierr)

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

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

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

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

260:       call PetscLogStagePop(ierr)
261:       CHKERRA(ierr)
262:       call PetscLogStagePush(stages(1),ierr)
263:       CHKERRA(ierr)

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

268:       call MatZeroEntries(C,ierr)
269:       CHKERRA(ierr)

271:       ! Assemble matrix again.  Note that we retain the same matrix data
272:       ! structure and the same nonzero pattern; we just change the values
273:       ! of the matrix entries.

275:       do i=0,m-1
276:         do j=2*rank,2*rank+1
277:           v =-1.0; Ii=j + n*i
278:           if (i>0) then
279:             JJ = Ii - n
280:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
281:             CHKERRA(ierr)
282:           endif

284:           if (i<m-1) then
285:             JJ = Ii + n
286:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
287:             CHKERRA(ierr)
288:           endif

290:           if (j>0) then
291:             JJ = Ii - 1
292:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
293:             CHKERRA(ierr)
294:           endif

296:           if (j<n-1) then
297:             JJ = Ii + 1
298:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
299:             CHKERRA(ierr)
300:           endif

302:           v=6.0
303:           call MatSetValues(C,one,Ii,one,Ii,v,ADD_VALUES,ierr)
304:           CHKERRA(ierr)

306:         enddo
307:       enddo

309:       ! Make the matrix nonsymmetric if desired

311:       if (mat_nonsymmetric) then
312:         do Ii=Istart,Iend-1
313:           v=-1.5;  i=Ii/n
314:           if (i>1) then
315:             JJ=Ii-n-1
316:             call MatSetValues(C,one,Ii,one,JJ,v,ADD_VALUES,ierr)
317:             CHKERRA(ierr)
318:           endif
319:         enddo
320:       endif

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

327:       call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
328:       CHKERRA(ierr)
329:       call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
330:       CHKERRA(ierr)

332:       if (testscaledMat) then
333:         ! Scale a(0,0) and a(M-1,M-1)

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

343:         endif

345:         call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
346:         call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)

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

350:         call  VecDestroy(u,ierr)
351:         call  VecCreate(PETSC_COMM_WORLD,u,ierr)
352:         call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
353:         call  VecSetFromOptions(u,ierr)

355:         call  PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
356:         call  PetscRandomSetFromOptions(rctx,ierr)
357:         call  VecSetRandom(u,rctx,ierr)
358:         call  PetscRandomDestroy(rctx,ierr)
359:         call  VecAssemblyBegin(u,ierr)
360:         call  VecAssemblyEnd(u,ierr)

362:       endif

364:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_newMat",testnewC,flg,ierr)
365:       CHKERRA(ierr)

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

371:         call  MatDuplicate(C,MAT_COPY_VALUES,Ctmp,ierr)
372:         call  MatDestroy(C,ierr)
373:         call  MatDuplicate(Ctmp,MAT_COPY_VALUES,C,ierr)
374:         call  MatDestroy(Ctmp,ierr)
375:       endif

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

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

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

384:       ! Solve linear system

386:       call  KSPSetUp(ksp,ierr); CHKERRA(ierr)
387:       call  KSPSolve(ksp,b,x,ierr); CHKERRA(ierr)

389:       ! Check the error

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

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

402:       call  KSPDestroy(ksp,ierr); CHKERRA(ierr)
403:       call  VecDestroy(u,ierr); CHKERRA(ierr)
404:       call  VecDestroy(x,ierr); CHKERRA(ierr)
405:       call  VecDestroy(b,ierr); CHKERRA(ierr)
406:       call  MatDestroy(C,ierr); CHKERRA(ierr)

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

410:       call PetscLogStagePop(ierr)
411:       CHKERRA(ierr)

413:       call PetscFinalize(ierr)

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

504: end program main