Actual source code: ex14f.F
petsc-3.5.4 2015-05-23
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