Actual source code: ex5f.F90

petsc-3.11.2 2019-05-18
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
 15: 
 16:       implicit none
 17:       KSP            :: myKsp            ! 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, parameter :: m = 3
 27:       PetscMPIInt    :: mySize,myRank
 28:       PetscBool, parameter :: &
 29:         testnewC         = PETSC_FALSE, &
 30:         testscaledMat    = PETSC_FALSE, &
 31:         mat_nonsymmetric = PETSC_FALSE
 32:       PetscBool      :: flg
 33:       PetscRandom    :: rctx
 34:       PetscLogStage,dimension(2) :: stages
 35:       character(len=80)          :: outputString
 36: 
 37:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
 38:       if (ierr /= 0) then
 39:         write(6,*)'Unable to initialize PETSc'
 40:         stop
 41:       endif
 42: 
 43:       call PetscOptionsGetInt(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,'-m',m,flg,ierr)
 44:       CHKERRA(ierr)
 45:       call MPI_Comm_rank(PETSC_COMM_WORLD,myRank,ierr)
 46:       CHKERRA(ierr)
 47:       call MPI_Comm_size(PETSC_COMM_WORLD,mySize,ierr)
 48:       CHKERRA(ierr)
 49:       n=2*mySize

 51: 
 52:       ! Set flag if we are doing a nonsymmetric problem; the default is symmetric.
 53: 
 54:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-mat_nonsym",mat_nonsymmetric,flg,ierr)
 55:       CHKERRA(ierr)
 56:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_scaledMat",testscaledMat,flg,ierr)
 57:       CHKERRA(ierr)
 58: 
 59:       ! Register two stages for separate profiling of the two linear solves.
 60:       ! Use the runtime option -log_view for a printout of performance
 61:       ! statistics at the program's conlusion.

 63:       call PetscLogStageRegister("Original Solve",stages(1),ierr)
 64:       CHKERRA(ierr)
 65:       call PetscLogStageRegister("Second Solve",stages(2),ierr)
 66:       CHKERRA(ierr)
 67: 
 68:       ! -------------- Stage 0: Solve Original System ----------------------
 69:       ! Indicate to PETSc profiling that we're beginning the first stage
 70: 
 71:       call PetscLogStagePush(stages(1),ierr)
 72:       CHKERRA(ierr)
 73: 
 74:       ! Create parallel matrix, specifying only its global dimensions.
 75:       ! When using MatCreate(), the matrix format can be specified at
 76:       ! runtime. Also, the parallel partitioning of the matrix is
 77:       ! determined by PETSc at runtime.
 78: 
 79:       call MatCreate(PETSC_COMM_WORLD,C,ierr)
 80:       CHKERRA(ierr)
 81:       call MatSetSizes(C,PETSC_DECIDE,PETSC_DECIDE,m*n,m*n,ierr)
 82:       CHKERRA(ierr)
 83:       call MatSetFromOptions(C,ierr)
 84:       CHKERRA(ierr)
 85:       call MatSetUp(C,ierr)
 86:       CHKERRA(ierr)

 88:       ! Currently, all PETSc parallel matrix formats are partitioned by
 89:       ! contiguous chunks of rows across the processors.  Determine which
 90:       ! rows of the matrix are locally owned.
 91: 
 92:       call MatGetOwnershipRange(C,Istart,Iend,ierr)
 93: 
 94:       ! Set matrix entries matrix in parallel.
 95:       ! - Each processor needs to insert only elements that it owns
 96:       ! locally (but any non-local elements will be sent to the
 97:       ! appropriate processor during matrix assembly).
 98:       !- Always specify global row and columns of matrix entries.

100:       intitializeC: do Ii=Istart,Iend-1
101:           v =-1.0; i = Ii/n; j = Ii - i*n
102:           if (i>0) then
103:             JJ = Ii - n
104:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
105:             CHKERRA(ierr)
106:           endif
107: 
108:           if (i<m-1) then
109:             JJ = Ii + n
110:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
111:             CHKERRA(ierr)
112:           endif
113: 
114:           if (j>0) then
115:             JJ = Ii - 1
116:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
117:             CHKERRA(ierr)
118:           endif
119: 
120:           if (j<n-1) then
121:             JJ = Ii + 1
122:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
123:             CHKERRA(ierr)
124:           endif
125: 
126:           v=4.0
127:           call MatSetValues(C,1,Ii,1,Ii,v,ADD_VALUES,ierr)
128:           CHKERRA(ierr)
129: 
130:       enddo intitializeC
131: 
132:       ! Make the matrix nonsymmetric if desired
133:       if (mat_nonsymmetric) then
134:         do Ii=Istart,Iend-1
135:           v=-1.5; i=Ii/n
136:           if (i>1) then
137:             JJ=Ii-n-1
138:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
139:             CHKERRA(ierr)
140:           endif
141:         enddo
142:       else
143:         call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
144:         CHKERRA(ierr)
145:         call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
146:         CHKERRA(ierr)
147:       endif
148: 
149:       ! Assemble matrix, using the 2-step process:
150:       ! MatAssemblyBegin(), MatAssemblyEnd()
151:       ! Computations can be done while messages are in transition
152:       ! by placing code between these two statements.
153: 
154:       call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
155:       CHKERRA(ierr)
156:       call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
157:       CHKERRA(ierr)

159:       ! Create parallel vectors.
160:       ! - When using VecSetSizes(), we specify only the vector's global
161:       !   dimension; the parallel partitioning is determined at runtime.
162:       ! - Note: We form 1 vector from scratch and then duplicate as needed.
163: 
164:       call  VecCreate(PETSC_COMM_WORLD,u,ierr)
165:       call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
166:       call  VecSetFromOptions(u,ierr)
167:       call  VecDuplicate(u,b,ierr)
168:       call  VecDuplicate(b,x,ierr)

170:       ! Currently, all parallel PETSc vectors are partitioned by
171:       ! contiguous chunks across the processors.  Determine which
172:       ! range of entries are locally owned.
173: 
174:       call  VecGetOwnershipRange(x,low,high,ierr)
175:       CHKERRA(ierr)

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

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

192:       ! Assemble vector, using the 2-step process:
193:       ! VecAssemblyBegin(), VecAssemblyEnd()
194:       ! Computations can be done while messages are in transition,
195:       ! by placing code between these two statements.
196: 
197:       call  VecAssemblyBegin(u,ierr)
198:       CHKERRA(ierr)
199:       call  VecAssemblyEnd(u,ierr)
200:       CHKERRA(ierr)
201: 
202:       ! Compute right-hand-side vector
203: 
204:       call  MatMult(C,u,b,ierr)

206:       CHKERRA(ierr)
207: 
208:       ! Create linear solver context

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

215:       call  KSPSetOperators(myKsp,C,C,ierr)
216:       CHKERRA(ierr)
217:       ! Set runtime options (e.g., -ksp_type <type> -pc_type <type>)

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

227:       call  KSPSetUp(myKsp,ierr)
228:       CHKERRA(ierr)

230:       call  KSPSolve(myKsp,b,x,ierr)
231:       CHKERRA(ierr)

233:       ! Check the error
234: 
235:       call VecAXPY(x,myNone,u,ierr)
236:       call VecNorm(x,NORM_2,norm,ierr)

238:       call KSPGetIterationNumber(myKsp,its,ierr)
239:       if (.not. testscaledMat .or. norm > 1.e-7) then
240:         write(outputString,*)'Norm of error',norm,'Iterations',its,'\n'
241:         call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
242:       endif
243: 
244:       ! -------------- Stage 1: Solve Second System ----------------------
245: 
246:       ! Solve another linear system with the same method.  We reuse the KSP
247:       ! context, matrix and vector data structures, and hence save the
248:       ! overhead of creating new ones.

250:       ! Indicate to PETSc profiling that we're concluding the first
251:       ! stage with PetscLogStagePop(), and beginning the second stage with
252:       ! PetscLogStagePush().
253: 
254:       call PetscLogStagePop(ierr)
255:       CHKERRA(ierr)
256:       call PetscLogStagePush(stages(2),ierr)
257:       CHKERRA(ierr)
258: 
259:       ! Initialize all matrix entries to zero.  MatZeroEntries() retains the
260:       ! nonzero structure of the matrix for sparse formats.
261: 
262:       call MatZeroEntries(C,ierr)
263:       CHKERRA(ierr)
264: 
265:       ! Assemble matrix again.  Note that we retain the same matrix data
266:       ! structure and the same nonzero pattern; we just change the values
267:       ! of the matrix entries.
268: 
269:       do i=0,m-1
270:         do j=2*myRank,2*myRank+1
271:           v =-1.0; Ii=j + n*i
272:           if (i>0) then
273:             JJ = Ii - n
274:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
275:             CHKERRA(ierr)
276:           endif
277: 
278:           if (i<m-1) then
279:             JJ = Ii + n
280:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
281:             CHKERRA(ierr)
282:           endif
283: 
284:           if (j>0) then
285:             JJ = Ii - 1
286:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
287:             CHKERRA(ierr)
288:           endif
289: 
290:           if (j<n-1) then
291:             JJ = Ii + 1
292:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
293:             CHKERRA(ierr)
294:           endif
295: 
296:           v=6.0
297:           call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
298:           CHKERRA(ierr)
299: 
300:         enddo
301:       enddo
302: 
303:       ! Make the matrix nonsymmetric if desired
304: 
305:       if (mat_nonsymmetric) then
306:         do Ii=Istart,Iend-1
307:           v=-1.5;  i=Ii/n
308:           if (i>1) then
309:             JJ=Ii-n-1
310:             call MatSetValues(C,1,Ii,1,JJ,v,ADD_VALUES,ierr)
311:             CHKERRA(ierr)
312:           endif
313:         enddo
314:       else
315:           call MatSetOption(C,MAT_SYMMETRIC,PETSC_TRUE,ierr)
316:           CHKERRA(ierr)
317:           call MatSetOption(C,MAT_SYMMETRY_ETERNAL,PETSC_TRUE,ierr)
318:           CHKERRA(ierr)
319:       endif
320: 
321:       ! Assemble matrix, using the 2-step process:
322:       ! MatAssemblyBegin(), MatAssemblyEnd()
323:       ! Computations can be done while messages are in transition
324:       ! by placing code between these two statements.
325: 
326:       call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
327:       CHKERRA(ierr)
328:       call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
329:       CHKERRA(ierr)
330: 
331:       if (testscaledMat) then
332:         ! Scale a(0,0) and a(M-1,M-1)

334:         if (myRank /= 0) then
335:           v = 6.0*0.00001; Ii = 0; JJ = 0
336:           call MatSetValues(C,1,Ii,1,JJ,v,INSERT_VALUES,ierr)
337:           CHKERRA(ierr)
338:         elseif (myRank == mySize -1) then
339:           v = 6.0*0.00001; Ii = m*n-1; JJ = m*n-1
340:           call MatSetValues(C,1,Ii,1,JJ,v,INSERT_VALUES,ierr)
341: 
342:         endif
343: 
344:         call MatAssemblyBegin(C,MAT_FINAL_ASSEMBLY,ierr)
345:         call MatAssemblyEnd(C,MAT_FINAL_ASSEMBLY,ierr)
346: 
347:         ! Compute a new right-hand-side vector
348: 
349:         call  VecDestroy(u,ierr)
350:         call  VecCreate(PETSC_COMM_WORLD,u,ierr)
351:         call  VecSetSizes(u,PETSC_DECIDE,m*n,ierr)
352:         call  VecSetFromOptions(u,ierr)

354:         call  PetscRandomCreate(PETSC_COMM_WORLD,rctx,ierr)
355:         call  PetscRandomSetFromOptions(rctx,ierr)
356:         call  VecSetRandom(u,rctx,ierr)
357:         call  PetscRandomDestroy(rctx,ierr)
358:         call  VecAssemblyBegin(u,ierr)
359:         call  VecAssemblyEnd(u,ierr)
360: 
361:       endif
362: 
363:       call PetscOptionsGetBool(PETSC_NULL_OPTIONS,PETSC_NULL_CHARACTER,"-test_newMat",testnewC,flg,ierr)
364:       CHKERRA(ierr)
365: 
366:       if (testnewC) then
367:       ! User may use a new matrix C with same nonzero pattern, e.g.
368:       ! ex5 -ksp_monitor -mat_type sbaij -pc_type cholesky -pc_factor_mat_solver_type mumps -test_newMat

370:         call  MatDuplicate(C,MAT_COPY_VALUES,Ctmp,ierr)
371:         call  MatDestroy(C,ierr)
372:         call  MatDuplicate(Ctmp,MAT_COPY_VALUES,C,ierr)
373:         call  MatDestroy(Ctmp,ierr)
374:       endif
375: 
376:       call MatMult(C,u,b,ierr);CHKERRA(ierr)

378:       ! Set operators. Here the matrix that defines the linear system
379:       ! also serves as the preconditioning matrix.
380: 
381:       call KSPSetOperators(myKsp,C,C,ierr);CHKERRA(ierr)
382: 
383:       ! Solve linear system
384: 
385:       call  KSPSetUp(myKsp,ierr); CHKERRA(ierr)
386:       call  KSPSolve(myKsp,b,x,ierr); CHKERRA(ierr)
387: 
388:       ! Check the error
389: 
390:       call VecAXPY(x,myNone,u,ierr); CHKERRA(ierr)
391:       call VecNorm(x,NORM_2,norm,ierr); CHKERRA(ierr)
392:       call KSPGetIterationNumber(myKsp,its,ierr); CHKERRA(ierr)
393:       if (.not. testscaledMat .or. norm > 1.e-7) then
394:         write(outputString,*)'Norm of error',norm,'Iterations',its,'\n'
395:         call PetscPrintf(PETSC_COMM_WORLD,outputString,ierr)
396:       endif
397: 
398:       ! Free work space.  All PETSc objects should be destroyed when they
399:       ! are no longer needed.
400: 
401:       call  KSPDestroy(myKsp,ierr); CHKERRA(ierr)
402:       call  VecDestroy(u,ierr); CHKERRA(ierr)
403:       call  VecDestroy(x,ierr); CHKERRA(ierr)
404:       call  VecDestroy(b,ierr); CHKERRA(ierr)
405:       call  MatDestroy(C,ierr); CHKERRA(ierr)
406: 
407:       ! Indicate to PETSc profiling that we're concluding the second stage
408: 
409:       call PetscLogStagePop(ierr)
410:       CHKERRA(ierr)
411: 
412:       call PetscFinalize(ierr)
413: 

415: end program main