Actual source code: ex14f.F

petsc-3.5.4 2015-05-23
Report Typos and Errors
```  1: !
2: !
3: !  Solves a nonlinear system in parallel with a user-defined
4: !  Newton method that uses KSP to solve the linearized Newton sytems.  This solver
5: !  is a very simplistic inexact Newton method.  The intent of this code is to
6: !  demonstrate the repeated solution of linear sytems with the same nonzero pattern.
7: !
8: !  This is NOT the recommended approach for solving nonlinear problems with PETSc!
9: !  We urge users to employ the SNES component for solving nonlinear problems whenever
10: !  possible, as it offers many advantages over coding nonlinear solvers independently.
11: !
12: !  We solve the  Bratu (SFI - solid fuel ignition) problem in a 2D rectangular
13: !  domain, using distributed arrays (DMDAs) to partition the parallel grid.
14: !
15: !  The command line options include:
16: !  -par <parameter>, where <parameter> indicates the problem's nonlinearity
17: !     problem SFI:  <parameter> = Bratu parameter (0 <= par <= 6.81)
18: !  -mx <xg>, where <xg> = number of grid points in the x-direction
19: !  -my <yg>, where <yg> = number of grid points in the y-direction
20: !  -Nx <npx>, where <npx> = number of processors in the x-direction
21: !  -Ny <npy>, where <npy> = number of processors in the y-direction
22: !  -mf use matrix free for matrix vector product
23: !
24: !/*T
25: !   Concepts: KSP^writing a user-defined nonlinear solver
26: !   Concepts: DMDA^using distributed arrays
27: !   Processors: n
28: !T*/
29: !  ------------------------------------------------------------------------
30: !
31: !    Solid Fuel Ignition (SFI) problem.  This problem is modeled by
32: !    the partial differential equation
33: !
34: !            -Laplacian u - lambda*exp(u) = 0,  0 < x,y < 1,
35: !
36: !    with boundary conditions
37: !
38: !             u = 0  for  x = 0, x = 1, y = 0, y = 1.
39: !
40: !    A finite difference approximation with the usual 5-point stencil
41: !    is used to discretize the boundary value problem to obtain a nonlinear
42: !    system of equations.
43: !
44: !    The SNES version of this problem is:  snes/examples/tutorials/ex5f.F
45: !
46: !  -------------------------------------------------------------------------

48:       program main
49:       implicit none

51: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
52: !                    Include files
53: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
54: !
55: !     petscsys.h       - base PETSc routines   petscvec.h - vectors
56: !     petscmat.h - matrices
57: !     petscis.h     - index sets            petscksp.h - Krylov subspace methods
58: !     petscviewer.h - viewers               petscpc.h  - preconditioners

60: #include <finclude/petscsys.h>
61: #include <finclude/petscis.h>
62: #include <finclude/petscvec.h>
63: #include <finclude/petscmat.h>
64: #include <finclude/petscpc.h>
65: #include <finclude/petscksp.h>
66: #include <finclude/petscdm.h>
67: #include <finclude/petscdmda.h>

69:       MPI_Comm comm
70:       Vec      X,Y,F,localX,localF
71:       Mat      J,B
72:       DM       da
73:       KSP      ksp

75:       PetscInt  Nx,Ny,N,mx,my,ifive,ithree
76:       PetscBool  flg,nooutput,usemf
77:       common   /mycommon/ mx,my,B,localX,localF,da
78: !
79: !
80: !      This is the routine to use for matrix-free approach
81: !
82:       external mymult

84: !     --------------- Data to define nonlinear solver --------------
85:       PetscReal   rtol,ttol
86:       PetscReal   fnorm,ynorm,xnorm
87:       PetscInt            max_nonlin_its,one
88:       PetscInt            lin_its
89:       PetscInt           i,m
90:       PetscScalar        mone
91:       PetscErrorCode ierr

93:       mone           = -1.d0
94:       rtol           = 1.d-8
95:       max_nonlin_its = 10
96:       one            = 1
97:       ifive          = 5
98:       ithree         = 3

100:       call PetscInitialize(PETSC_NULL_CHARACTER,ierr)
101:       comm = PETSC_COMM_WORLD

103: !  Initialize problem parameters

105: !
106:       mx = 4
107:       my = 4
108:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-mx',mx,flg,ierr)
109:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-my',my,flg,ierr)
110:       N = mx*my

112:       nooutput = .false.
113:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-no_output',       &
114:      &     nooutput,ierr)

116: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
117: !     Create linear solver context
118: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

120:       call KSPCreate(comm,ksp,ierr)

122: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
123: !     Create vector data structures
124: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

126: !
127: !  Create distributed array (DMDA) to manage parallel grid and vectors
128: !
129:       Nx = PETSC_DECIDE
130:       Ny = PETSC_DECIDE
131:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-Nx',Nx,flg,ierr)
132:       call PetscOptionsGetInt(PETSC_NULL_CHARACTER,'-Ny',Ny,flg,ierr)
133:       call DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,         &
134:      &     DMDA_STENCIL_STAR,mx,my,Nx,Ny,one,one,                        &
135:      &     PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)

137: !
138: !  Extract global and local vectors from DMDA then duplicate for remaining
139: !  vectors that are the same types
140: !
141:        call DMCreateGlobalVector(da,X,ierr)
142:        call DMCreateLocalVector(da,localX,ierr)
143:        call VecDuplicate(X,F,ierr)
144:        call VecDuplicate(X,Y,ierr)
145:        call VecDuplicate(localX,localF,ierr)

148: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
149: !     Create matrix data structure for Jacobian
150: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151: !
152: !     Note:  For the parallel case, vectors and matrices MUST be partitioned
153: !     accordingly.  When using distributed arrays (DMDAs) to create vectors,
154: !     the DMDAs determine the problem partitioning.  We must explicitly
155: !     specify the local matrix dimensions upon its creation for compatibility
156: !     with the vector distribution.
157: !
158: !     Note: Here we only approximately preallocate storage space for the
159: !     Jacobian.  See the users manual for a discussion of better techniques
160: !     for preallocating matrix memory.
161: !
162:       call VecGetLocalSize(X,m,ierr)
163:       call MatCreateAIJ(comm,m,m,N,N,ifive,PETSC_NULL_INTEGER,ithree,         &
164:      &     PETSC_NULL_INTEGER,B,ierr)

166: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167: !     if usemf is on then matrix vector product is done via matrix free
168: !     approach. Note this is just an example, and not realistic because
169: !     we still use the actual formed matrix, but in reality one would
170: !     provide their own subroutine that would directly do the matrix
171: !     vector product and not call MatMult()
172: !     Note: we put B into a common block so it will be visible to the
173: !     mymult() routine
174:       usemf = .false.
175:       call PetscOptionsHasName(PETSC_NULL_CHARACTER,'-mf',usemf,ierr)
176:       if (usemf) then
177:          call MatCreateShell(comm,m,m,N,N,PETSC_NULL_INTEGER,J,ierr)
178:          call MatShellSetOperation(J,MATOP_MULT,mymult,ierr)
179:       else
180: !        If not doing matrix free then matrix operator, J,  and matrix used
181: !        to construct preconditioner, B, are the same
182:         J = B
183:       endif

185: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: !     Customize linear solver set runtime options
187: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
188: !
189: !     Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
190: !
191:        call KSPSetFromOptions(ksp,ierr)

193: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
194: !     Evaluate initial guess
195: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

197:        call FormInitialGuess(X,ierr)
198:        call ComputeFunction(X,F,ierr)
199:        call VecNorm(F,NORM_2,fnorm,ierr)
200:        ttol = fnorm*rtol
201:        if (.not. nooutput) then
202:          print*, 'Initial function norm ',fnorm
203:        endif

205: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
206: !     Solve nonlinear system with a user-defined method
207: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

209: !  This solver is a very simplistic inexact Newton method, with no
210: !  no damping strategies or bells and whistles. The intent of this code
211: !  is merely to demonstrate the repeated solution with KSP of linear
212: !  sytems with the same nonzero structure.
213: !
214: !  This is NOT the recommended approach for solving nonlinear problems
215: !  with PETSc!  We urge users to employ the SNES component for solving
216: !  nonlinear problems whenever possible with application codes, as it
217: !  offers many advantages over coding nonlinear solvers independently.

219:        do 10 i=0,max_nonlin_its

221: !  Compute the Jacobian matrix.  See the comments in this routine for
222: !  important information about setting the flag mat_flag.

224:          call ComputeJacobian(X,B,ierr)

226: !  Solve J Y = F, where J is the Jacobian matrix.
227: !    - First, set the KSP linear operators.  Here the matrix that
228: !      defines the linear system also serves as the preconditioning
229: !      matrix.
230: !    - Then solve the Newton system.

232:          call KSPSetOperators(ksp,J,B,ierr)
233:          call KSPSolve(ksp,F,Y,ierr)

235: !  Compute updated iterate

237:          call VecNorm(Y,NORM_2,ynorm,ierr)
238:          call VecAYPX(Y,mone,X,ierr)
239:          call VecCopy(Y,X,ierr)
240:          call VecNorm(X,NORM_2,xnorm,ierr)
241:          call KSPGetIterationNumber(ksp,lin_its,ierr)
242:          if (.not. nooutput) then
243:            print*,'linear solve iterations = ',lin_its,' xnorm = ',     &
244:      &         xnorm,' ynorm = ',ynorm
245:          endif

247: !  Evaluate nonlinear function at new location

249:          call ComputeFunction(X,F,ierr)
250:          call VecNorm(F,NORM_2,fnorm,ierr)
251:          if (.not. nooutput) then
252:            print*, 'Iteration ',i+1,' function norm',fnorm
253:          endif

255: !  Test for convergence

257:        if (fnorm .le. ttol) then
258:          if (.not. nooutput) then
259:            print*,'Converged: function norm ',fnorm,' tolerance ',ttol
260:          endif
261:          goto 20
262:        endif
263:  10   continue
264:  20   continue

266:       write(6,100) i+1
267:  100  format('Number of SNES iterations =',I2)

269: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
270: !     Free work space.  All PETSc objects should be destroyed when they
271: !     are no longer needed.
272: !  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

274:        call MatDestroy(B,ierr)
275:        if (usemf) then
276:          call MatDestroy(J,ierr)
277:        endif
278:        call VecDestroy(localX,ierr)
279:        call VecDestroy(X,ierr)
280:        call VecDestroy(Y,ierr)
281:        call VecDestroy(localF,ierr)
282:        call VecDestroy(F,ierr)
283:        call KSPDestroy(ksp,ierr)
284:        call DMDestroy(da,ierr)
285:        call PetscFinalize(ierr)
286:        end

288: ! -------------------------------------------------------------------
289: !
290: !   FormInitialGuess - Forms initial approximation.
291: !
292: !   Input Parameters:
293: !   X - vector
294: !
295: !   Output Parameter:
296: !   X - vector
297: !
298:       subroutine FormInitialGuess(X,ierr)
299:       implicit none

301: !     petscsys.h       - base PETSc routines   petscvec.h - vectors
302: !     petscmat.h - matrices
303: !     petscis.h     - index sets            petscksp.h - Krylov subspace methods
304: !     petscviewer.h - viewers               petscpc.h  - preconditioners

306: #include <finclude/petscsys.h>
307: #include <finclude/petscis.h>
308: #include <finclude/petscvec.h>
309: #include <finclude/petscmat.h>
310: #include <finclude/petscpc.h>
311: #include <finclude/petscksp.h>
312: #include <finclude/petscdm.h>
313: #include <finclude/petscdmda.h>
314:       PetscErrorCode    ierr
315:       PetscOffset      idx
316:       Vec       X,localX,localF
317:       PetscInt  i,j,row,mx
318:       PetscInt  my, xs,ys,xm
319:       PetscInt  ym,gxm,gym,gxs,gys
320:       PetscReal one,lambda,temp1,temp,hx,hy
321:       PetscScalar      xx(1)
322:       DM               da
323:       Mat              B
324:       common   /mycommon/ mx,my,B,localX,localF,da

326:       one    = 1.d0
327:       lambda = 6.d0
328:       hx     = one/(mx-1)
329:       hy     = one/(my-1)
330:       temp1  = lambda/(lambda + one)

332: !  Get a pointer to vector data.
333: !    - VecGetArray() returns a pointer to the data array.
334: !    - You MUST call VecRestoreArray() when you no longer need access to
335: !      the array.
336:        call VecGetArray(localX,xx,idx,ierr)

338: !  Get local grid boundaries (for 2-dimensional DMDA):
339: !    xs, ys   - starting grid indices (no ghost points)
340: !    xm, ym   - widths of local grid (no ghost points)
341: !    gxs, gys - starting grid indices (including ghost points)
342: !    gxm, gym - widths of local grid (including ghost points)

344:        call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,             &
345:      &      PETSC_NULL_INTEGER,ierr)
346:        call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,    &
347:      &      PETSC_NULL_INTEGER,ierr)

349: !  Compute initial guess over the locally owned part of the grid

351:       do 30 j=ys,ys+ym-1
352:         temp = (min(j,my-j-1))*hy
353:         do 40 i=xs,xs+xm-1
354:           row = i - gxs + (j - gys)*gxm + 1
355:           if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or.              &
356:      &        j .eq. my-1) then
357:             xx(idx+row) = 0.d0
358:             continue
359:           endif
360:           xx(idx+row) = temp1*sqrt(min((min(i,mx-i-1))*hx,temp))
361:  40     continue
362:  30   continue

364: !     Restore vector

366:        call VecRestoreArray(localX,xx,idx,ierr)

368: !     Insert values into global vector

370:        call DMLocalToGlobalBegin(da,localX,INSERT_VALUES,X,ierr)
371:        call DMLocalToGlobalEnd(da,localX,INSERT_VALUES,X,ierr)
372:        return
373:        end

375: ! -------------------------------------------------------------------
376: !
377: !   ComputeFunction - Evaluates nonlinear function, F(x).
378: !
379: !   Input Parameters:
380: !.  X - input vector
381: !
382: !   Output Parameter:
383: !.  F - function vector
384: !
385:       subroutine  ComputeFunction(X,F,ierr)
386:       implicit none

388: !     petscsys.h       - base PETSc routines   petscvec.h - vectors
389: !     petscmat.h - matrices
390: !     petscis.h     - index sets            petscksp.h - Krylov subspace methods
391: !     petscviewer.h - viewers               petscpc.h  - preconditioners

393: #include <finclude/petscsys.h>
394: #include <finclude/petscis.h>
395: #include <finclude/petscvec.h>
396: #include <finclude/petscmat.h>
397: #include <finclude/petscpc.h>
398: #include <finclude/petscksp.h>
399: #include <finclude/petscdm.h>
400: #include <finclude/petscdmda.h>

402:       Vec              X,F,localX,localF
403:       PetscInt         gys,gxm,gym
404:       PetscOffset      idx,idf
405:       PetscErrorCode ierr
406:       PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxs
407:       PetscReal two,one,lambda,hx
408:       PetscReal hy,hxdhy,hydhx,sc
409:       PetscScalar      u,uxx,uyy,xx(1),ff(1)
410:       DM               da
411:       Mat              B
412:       common   /mycommon/ mx,my,B,localX,localF,da

414:       two    = 2.d0
415:       one    = 1.d0
416:       lambda = 6.d0

418:       hx     = one/(mx-1)
419:       hy     = one/(my-1)
420:       sc     = hx*hy*lambda
421:       hxdhy  = hx/hy
422:       hydhx  = hy/hx

424: !  Scatter ghost points to local vector, using the 2-step process
425: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
426: !  By placing code between these two statements, computations can be
427: !  done while messages are in transition.
428: !
429:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
430:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)

432: !  Get pointers to vector data

434:       call VecGetArray(localX,xx,idx,ierr)
435:       call VecGetArray(localF,ff,idf,ierr)

437: !  Get local grid boundaries

439:       call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
440:      &     PETSC_NULL_INTEGER,ierr)
441:       call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
442:      &     PETSC_NULL_INTEGER,ierr)

444: !  Compute function over the locally owned part of the grid

446:       do 50 j=ys,ys+ym-1

448:         row = (j - gys)*gxm + xs - gxs
449:         do 60 i=xs,xs+xm-1
450:           row = row + 1

452:           if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or.              &
453:      &        j .eq. my-1) then
454:             ff(idf+row) = xx(idx+row)
455:             goto 60
456:           endif
457:           u   = xx(idx+row)
458:           uxx = (two*u - xx(idx+row-1) - xx(idx+row+1))*hydhx
459:           uyy = (two*u - xx(idx+row-gxm) - xx(idx+row+gxm))*hxdhy
460:           ff(idf+row) = uxx + uyy - sc*exp(u)
461:  60     continue
462:  50   continue

464: !  Restore vectors

466:        call VecRestoreArray(localX,xx,idx,ierr)
467:        call VecRestoreArray(localF,ff,idf,ierr)

469: !  Insert values into global vector

471:        call DMLocalToGlobalBegin(da,localF,INSERT_VALUES,F,ierr)
472:        call DMLocalToGlobalEnd(da,localF,INSERT_VALUES,F,ierr)
473:        return
474:        end

476: ! -------------------------------------------------------------------
477: !
478: !   ComputeJacobian - Evaluates Jacobian matrix.
479: !
480: !   Input Parameters:
481: !   x - input vector
482: !
483: !   Output Parameters:
484: !   jac - Jacobian matrix
485: !   flag - flag indicating matrix structure
486: !
487: !   Notes:
488: !   Due to grid point reordering with DMDAs, we must always work
489: !   with the local grid points, and then transform them to the new
490: !   global numbering with the 'ltog' mapping
491: !   We cannot work directly with the global numbers for the original
492: !   uniprocessor grid!
493: !
494:       subroutine ComputeJacobian(X,jac,ierr)
495:       implicit none

497: !     petscsys.h  - base PETSc routines   petscvec.h - vectors
498: !     petscmat.h - matrices
499: !     petscis.h     - index sets            petscksp.h - Krylov subspace methods
500: !     petscviewer.h - viewers               petscpc.h  - preconditioners

502: #include <finclude/petscsys.h>
503: #include <finclude/petscis.h>
504: #include <finclude/petscvec.h>
505: #include <finclude/petscmat.h>
506: #include <finclude/petscpc.h>
507: #include <finclude/petscksp.h>
508: #include <finclude/petscdm.h>
509: #include <finclude/petscdmda.h>

511:       Vec         X
512:       Mat         jac
513:       Vec         localX,localF
514:       DM          da
515:       PetscInt     ltog(1)
516:       PetscOffset idltog,idx
517:       PetscErrorCode ierr
518:       PetscInt xs,ys,xm,ym
519:       PetscInt gxs,gys,gxm,gym
520:       PetscInt grow(1),i,j
521:       PetscInt row,mx,my,ione
522:       PetscInt col(5),ifive
523:       PetscScalar two,one,lambda
524:       PetscScalar v(5),hx,hy,hxdhy
525:       PetscScalar hydhx,sc,xx(1)
526:       Mat         B
527:       ISLocalToGlobalMapping ltogm
528:       common   /mycommon/ mx,my,B,localX,localF,da

530:       ione   = 1
531:       ifive  = 5
532:       one    = 1.d0
533:       two    = 2.d0
534:       hx     = one/(mx-1)
535:       hy     = one/(my-1)
536:       sc     = hx*hy
537:       hxdhy  = hx/hy
538:       hydhx  = hy/hx
539:       lambda = 6.d0

541: !  Scatter ghost points to local vector, using the 2-step process
542: !     DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
543: !  By placing code between these two statements, computations can be
544: !  done while messages are in transition.

546:       call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
547:       call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)

549: !  Get pointer to vector data

551:       call VecGetArray(localX,xx,idx,ierr)

553: !  Get local grid boundaries

555:       call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym,              &
556:      &     PETSC_NULL_INTEGER,ierr)
557:       call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym,     &
558:      &                        PETSC_NULL_INTEGER,ierr)

560: !  Get the global node numbers for all local nodes, including ghost points

562:       call DMGetLocalToGlobalMapping(da,ltogm,ierr)
563:       call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)

565: !  Compute entries for the locally owned part of the Jacobian.
566: !   - Currently, all PETSc parallel matrix formats are partitioned by
567: !     contiguous chunks of rows across the processors. The 'grow'
568: !     parameter computed below specifies the global row number
569: !     corresponding to each local grid point.
570: !   - Each processor needs to insert only elements that it owns
571: !     locally (but any non-local elements will be sent to the
572: !     appropriate processor during matrix assembly).
573: !   - Always specify global row and columns of matrix entries.
574: !   - Here, we set all entries for a particular row at once.

576:       do 10 j=ys,ys+ym-1
577:         row = (j - gys)*gxm + xs - gxs
578:         do 20 i=xs,xs+xm-1
579:           row = row + 1
580:           grow(1) = ltog(idltog+row)
581:           if (i .eq. 0 .or. j .eq. 0 .or. i .eq. (mx-1) .or.            &
582:      &        j .eq. (my-1)) then
583:              call MatSetValues(jac,ione,grow,ione,grow,one,             &
584:      &                         INSERT_VALUES,ierr)
585:              go to 20
586:           endif
587:           v(1)   = -hxdhy
588:           col(1) = ltog(idltog+row - gxm)
589:           v(2)   = -hydhx
590:           col(2) = ltog(idltog+row - 1)
591:           v(3)   = two*(hydhx + hxdhy) - sc*lambda*exp(xx(idx+row))
592:           col(3) = grow(1)
593:           v(4)   = -hydhx
594:           col(4) = ltog(idltog+row + 1)
595:           v(5)   = -hxdhy
596:           col(5) = ltog(idltog+row + gxm)
597:           call MatSetValues(jac,ione,grow,ifive,col,v,INSERT_VALUES,       &
598:      &                      ierr)
599:  20     continue
600:  10   continue

602:       call ISLocalToGlobalMappingRestoreIndices(ltogm,ltog,idltog,ierr)

604: !  Assemble matrix, using the 2-step process:
605: !    MatAssemblyBegin(), MatAssemblyEnd().
606: !  By placing code between these two statements, computations can be
607: !  done while messages are in transition.

609:       call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
610:       call VecRestoreArray(localX,xx,idx,ierr)
611:       call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
612:       return
613:       end

616: ! -------------------------------------------------------------------
617: !
618: !   MyMult - user provided matrix multiply
619: !
620: !   Input Parameters:
621: !.  X - input vector
622: !
623: !   Output Parameter:
624: !.  F - function vector
625: !
626:       subroutine  MyMult(J,X,F,ierr)
627:       implicit none
628:       Mat     J,B
629:       Vec     X,F
630:       PetscErrorCode ierr
631:       PetscInt mx,my
632:       DM      da
633:       Vec     localX,localF

635:       common   /mycommon/ mx,my,B,localX,localF,da
636: !
637: !       Here we use the actual formed matrix B; users would
638: !     instead write their own matrix vector product routine
639: !
640:       call MatMult(B,X,F,ierr)
641:       return
642:       end

```