Actual source code: ex14f.F
petsc-3.7.5 2017-01-01
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 <petsc/finclude/petscsys.h>
61: #include <petsc/finclude/petscis.h>
62: #include <petsc/finclude/petscvec.h>
63: #include <petsc/finclude/petscmat.h>
64: #include <petsc/finclude/petscpc.h>
65: #include <petsc/finclude/petscksp.h>
66: #include <petsc/finclude/petscdm.h>
67: #include <petsc/finclude/petscdmda.h>
69: MPI_Comm comm
70: Vec X,Y,F,localX
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,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.0
94: rtol = 1.e-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_OBJECT,PETSC_NULL_CHARACTER, &
109: & '-mx',mx,flg,ierr)
110: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
111: & '-my',my,flg,ierr)
112: N = mx*my
114: nooutput = .false.
115: call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
116: & '-no_output',nooutput,ierr)
118: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119: ! Create linear solver context
120: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: call KSPCreate(comm,ksp,ierr)
124: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125: ! Create vector data structures
126: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: !
129: ! Create distributed array (DMDA) to manage parallel grid and vectors
130: !
131: Nx = PETSC_DECIDE
132: Ny = PETSC_DECIDE
133: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
134: & '-Nx',Nx,flg,ierr)
135: call PetscOptionsGetInt(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
136: & '-Ny',Ny,flg,ierr)
137: call DMDACreate2d(comm,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE, &
138: & DMDA_STENCIL_STAR,mx,my,Nx,Ny,one,one, &
139: & PETSC_NULL_INTEGER,PETSC_NULL_INTEGER,da,ierr)
141: !
142: ! Extract global and local vectors from DMDA then duplicate for remaining
143: ! vectors that are the same types
144: !
145: call DMCreateGlobalVector(da,X,ierr)
146: call DMCreateLocalVector(da,localX,ierr)
147: call VecDuplicate(X,F,ierr)
148: call VecDuplicate(X,Y,ierr)
151: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
152: ! Create matrix data structure for Jacobian
153: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
154: !
155: ! Note: For the parallel case, vectors and matrices MUST be partitioned
156: ! accordingly. When using distributed arrays (DMDAs) to create vectors,
157: ! the DMDAs determine the problem partitioning. We must explicitly
158: ! specify the local matrix dimensions upon its creation for compatibility
159: ! with the vector distribution.
160: !
161: ! Note: Here we only approximately preallocate storage space for the
162: ! Jacobian. See the users manual for a discussion of better techniques
163: ! for preallocating matrix memory.
164: !
165: call VecGetLocalSize(X,m,ierr)
166: call MatCreateAIJ(comm,m,m,N,N,ifive,PETSC_NULL_INTEGER,ithree, &
167: & PETSC_NULL_INTEGER,B,ierr)
169: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170: ! if usemf is on then matrix vector product is done via matrix free
171: ! approach. Note this is just an example, and not realistic because
172: ! we still use the actual formed matrix, but in reality one would
173: ! provide their own subroutine that would directly do the matrix
174: ! vector product and not call MatMult()
175: ! Note: we put B into a common block so it will be visible to the
176: ! mymult() routine
177: usemf = .false.
178: call PetscOptionsHasName(PETSC_NULL_OBJECT,PETSC_NULL_CHARACTER, &
179: & '-mf',usemf,ierr)
180: if (usemf) then
181: call MatCreateShell(comm,m,m,N,N,PETSC_NULL_INTEGER,J,ierr)
182: call MatShellSetOperation(J,MATOP_MULT,mymult,ierr)
183: else
184: ! If not doing matrix free then matrix operator, J, and matrix used
185: ! to construct preconditioner, B, are the same
186: J = B
187: endif
189: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
190: ! Customize linear solver set runtime options
191: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: !
193: ! Set runtime options (e.g., -ksp_monitor -ksp_rtol <rtol> -ksp_type <type>)
194: !
195: call KSPSetFromOptions(ksp,ierr)
197: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: ! Evaluate initial guess
199: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201: call FormInitialGuess(X,ierr)
202: call ComputeFunction(X,F,ierr)
203: call VecNorm(F,NORM_2,fnorm,ierr)
204: ttol = fnorm*rtol
205: if (.not. nooutput) then
206: print*, 'Initial function norm ',fnorm
207: endif
209: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210: ! Solve nonlinear system with a user-defined method
211: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213: ! This solver is a very simplistic inexact Newton method, with no
214: ! no damping strategies or bells and whistles. The intent of this code
215: ! is merely to demonstrate the repeated solution with KSP of linear
216: ! sytems with the same nonzero structure.
217: !
218: ! This is NOT the recommended approach for solving nonlinear problems
219: ! with PETSc! We urge users to employ the SNES component for solving
220: ! nonlinear problems whenever possible with application codes, as it
221: ! offers many advantages over coding nonlinear solvers independently.
223: do 10 i=0,max_nonlin_its
225: ! Compute the Jacobian matrix. See the comments in this routine for
226: ! important information about setting the flag mat_flag.
228: call ComputeJacobian(X,B,ierr)
230: ! Solve J Y = F, where J is the Jacobian matrix.
231: ! - First, set the KSP linear operators. Here the matrix that
232: ! defines the linear system also serves as the preconditioning
233: ! matrix.
234: ! - Then solve the Newton system.
236: call KSPSetOperators(ksp,J,B,ierr)
237: call KSPSolve(ksp,F,Y,ierr)
239: ! Compute updated iterate
241: call VecNorm(Y,NORM_2,ynorm,ierr)
242: call VecAYPX(Y,mone,X,ierr)
243: call VecCopy(Y,X,ierr)
244: call VecNorm(X,NORM_2,xnorm,ierr)
245: call KSPGetIterationNumber(ksp,lin_its,ierr)
246: if (.not. nooutput) then
247: print*,'linear solve iterations = ',lin_its,' xnorm = ', &
248: & xnorm,' ynorm = ',ynorm
249: endif
251: ! Evaluate nonlinear function at new location
253: call ComputeFunction(X,F,ierr)
254: call VecNorm(F,NORM_2,fnorm,ierr)
255: if (.not. nooutput) then
256: print*, 'Iteration ',i+1,' function norm',fnorm
257: endif
259: ! Test for convergence
261: if (fnorm .le. ttol) then
262: if (.not. nooutput) then
263: print*,'Converged: function norm ',fnorm,' tolerance ',ttol
264: endif
265: goto 20
266: endif
267: 10 continue
268: 20 continue
270: write(6,100) i+1
271: 100 format('Number of SNES iterations =',I2)
273: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: ! Free work space. All PETSc objects should be destroyed when they
275: ! are no longer needed.
276: ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
278: call MatDestroy(B,ierr)
279: if (usemf) then
280: call MatDestroy(J,ierr)
281: endif
282: call VecDestroy(localX,ierr)
283: call VecDestroy(X,ierr)
284: call VecDestroy(Y,ierr)
285: call VecDestroy(F,ierr)
286: call KSPDestroy(ksp,ierr)
287: call DMDestroy(da,ierr)
288: call PetscFinalize(ierr)
289: end
291: ! -------------------------------------------------------------------
292: !
293: ! FormInitialGuess - Forms initial approximation.
294: !
295: ! Input Parameters:
296: ! X - vector
297: !
298: ! Output Parameter:
299: ! X - vector
300: !
301: subroutine FormInitialGuess(X,ierr)
302: implicit none
304: ! petscsys.h - base PETSc routines petscvec.h - vectors
305: ! petscmat.h - matrices
306: ! petscis.h - index sets petscksp.h - Krylov subspace methods
307: ! petscviewer.h - viewers petscpc.h - preconditioners
309: #include <petsc/finclude/petscsys.h>
310: #include <petsc/finclude/petscis.h>
311: #include <petsc/finclude/petscvec.h>
312: #include <petsc/finclude/petscmat.h>
313: #include <petsc/finclude/petscpc.h>
314: #include <petsc/finclude/petscksp.h>
315: #include <petsc/finclude/petscdm.h>
316: #include <petsc/finclude/petscdmda.h>
317: PetscErrorCode ierr
318: PetscOffset idx
319: Vec X,localX
320: PetscInt i,j,row,mx
321: PetscInt my, xs,ys,xm
322: PetscInt ym
323: PetscReal one,lambda,temp1,temp,hx,hy
324: PetscScalar xx(2)
325: DM da
326: Mat B
327: common /mycommon/ mx,my,B,localX,da
329: one = 1.0
330: lambda = 6.0
331: hx = one/(mx-1)
332: hy = one/(my-1)
333: temp1 = lambda/(lambda + one)
335: ! Get a pointer to vector data.
336: ! - VecGetArray() returns a pointer to the data array.
337: ! - You MUST call VecRestoreArray() when you no longer need access to
338: ! the array.
339: call VecGetArray(X,xx,idx,ierr)
341: ! Get local grid boundaries (for 2-dimensional DMDA):
342: ! xs, ys - starting grid indices (no ghost points)
343: ! xm, ym - widths of local grid (no ghost points)
345: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
346: & PETSC_NULL_INTEGER,ierr)
348: ! Compute initial guess over the locally owned part of the grid
350: do 30 j=ys,ys+ym-1
351: temp = (min(j,my-j-1))*hy
352: do 40 i=xs,xs+xm-1
353: row = i - xs + (j - ys)*xm + 1
354: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. &
355: & j .eq. my-1) then
356: xx(idx+row) = 0.0
357: continue
358: endif
359: xx(idx+row) = temp1*sqrt(min((min(i,mx-i-1))*hx,temp))
360: 40 continue
361: 30 continue
363: ! Restore vector
365: call VecRestoreArray(X,xx,idx,ierr)
366: return
367: end
369: ! -------------------------------------------------------------------
370: !
371: ! ComputeFunction - Evaluates nonlinear function, F(x).
372: !
373: ! Input Parameters:
374: !. X - input vector
375: !
376: ! Output Parameter:
377: !. F - function vector
378: !
379: subroutine ComputeFunction(X,F,ierr)
380: implicit none
382: ! petscsys.h - base PETSc routines petscvec.h - vectors
383: ! petscmat.h - matrices
384: ! petscis.h - index sets petscksp.h - Krylov subspace methods
385: ! petscviewer.h - viewers petscpc.h - preconditioners
387: #include <petsc/finclude/petscsys.h>
388: #include <petsc/finclude/petscis.h>
389: #include <petsc/finclude/petscvec.h>
390: #include <petsc/finclude/petscmat.h>
391: #include <petsc/finclude/petscpc.h>
392: #include <petsc/finclude/petscksp.h>
393: #include <petsc/finclude/petscdm.h>
394: #include <petsc/finclude/petscdmda.h>
396: Vec X,F,localX
397: PetscInt gys,gxm,gym
398: PetscOffset idx,idf
399: PetscErrorCode ierr
400: PetscInt i,j,row,mx,my,xs,ys,xm,ym,gxs
401: PetscInt rowf
402: PetscReal two,one,lambda,hx
403: PetscReal hy,hxdhy,hydhx,sc
404: PetscScalar u,uxx,uyy,xx(2),ff(2)
405: DM da
406: Mat B
407: common /mycommon/ mx,my,B,localX,da
409: two = 2.0
410: one = 1.0
411: lambda = 6.0
413: hx = one/(mx-1)
414: hy = one/(my-1)
415: sc = hx*hy*lambda
416: hxdhy = hx/hy
417: hydhx = hy/hx
419: ! Scatter ghost points to local vector, using the 2-step process
420: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
421: ! By placing code between these two statements, computations can be
422: ! done while messages are in transition.
423: !
424: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
425: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
427: ! Get pointers to vector data
429: call VecGetArray(localX,xx,idx,ierr)
430: call VecGetArray(F,ff,idf,ierr)
432: ! Get local grid boundaries
434: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
435: & PETSC_NULL_INTEGER,ierr)
436: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
437: & PETSC_NULL_INTEGER,ierr)
439: ! Compute function over the locally owned part of the grid
440: rowf = 0
441: do 50 j=ys,ys+ym-1
443: row = (j - gys)*gxm + xs - gxs
444: do 60 i=xs,xs+xm-1
445: row = row + 1
446: rowf = rowf + 1
448: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. mx-1 .or. &
449: & j .eq. my-1) then
450: ff(idf+rowf) = xx(idx+row)
451: goto 60
452: endif
453: u = xx(idx+row)
454: uxx = (two*u - xx(idx+row-1) - xx(idx+row+1))*hydhx
455: uyy = (two*u - xx(idx+row-gxm) - xx(idx+row+gxm))*hxdhy
456: ff(idf+rowf) = uxx + uyy - sc*exp(u)
457: 60 continue
458: 50 continue
460: ! Restore vectors
462: call VecRestoreArray(localX,xx,idx,ierr)
463: call VecRestoreArray(F,ff,idf,ierr)
464: return
465: end
467: ! -------------------------------------------------------------------
468: !
469: ! ComputeJacobian - Evaluates Jacobian matrix.
470: !
471: ! Input Parameters:
472: ! x - input vector
473: !
474: ! Output Parameters:
475: ! jac - Jacobian matrix
476: ! flag - flag indicating matrix structure
477: !
478: ! Notes:
479: ! Due to grid point reordering with DMDAs, we must always work
480: ! with the local grid points, and then transform them to the new
481: ! global numbering with the 'ltog' mapping
482: ! We cannot work directly with the global numbers for the original
483: ! uniprocessor grid!
484: !
485: subroutine ComputeJacobian(X,jac,ierr)
486: implicit none
488: ! petscsys.h - base PETSc routines petscvec.h - vectors
489: ! petscmat.h - matrices
490: ! petscis.h - index sets petscksp.h - Krylov subspace methods
491: ! petscviewer.h - viewers petscpc.h - preconditioners
493: #include <petsc/finclude/petscsys.h>
494: #include <petsc/finclude/petscis.h>
495: #include <petsc/finclude/petscvec.h>
496: #include <petsc/finclude/petscmat.h>
497: #include <petsc/finclude/petscpc.h>
498: #include <petsc/finclude/petscksp.h>
499: #include <petsc/finclude/petscdm.h>
500: #include <petsc/finclude/petscdmda.h>
502: Vec X
503: Mat jac
504: Vec localX
505: DM da
506: PetscInt ltog(2)
507: PetscOffset idltog,idx
508: PetscErrorCode ierr
509: PetscInt xs,ys,xm,ym
510: PetscInt gxs,gys,gxm,gym
511: PetscInt grow(1),i,j
512: PetscInt row,mx,my,ione
513: PetscInt col(5),ifive
514: PetscScalar two,one,lambda
515: PetscScalar v(5),hx,hy,hxdhy
516: PetscScalar hydhx,sc,xx(2)
517: Mat B
518: ISLocalToGlobalMapping ltogm
519: common /mycommon/ mx,my,B,localX,da
521: ione = 1
522: ifive = 5
523: one = 1.0
524: two = 2.0
525: hx = one/(mx-1)
526: hy = one/(my-1)
527: sc = hx*hy
528: hxdhy = hx/hy
529: hydhx = hy/hx
530: lambda = 6.0
532: ! Scatter ghost points to local vector, using the 2-step process
533: ! DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
534: ! By placing code between these two statements, computations can be
535: ! done while messages are in transition.
537: call DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX,ierr)
538: call DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX,ierr)
540: ! Get pointer to vector data
542: call VecGetArray(localX,xx,idx,ierr)
544: ! Get local grid boundaries
546: call DMDAGetCorners(da,xs,ys,PETSC_NULL_INTEGER,xm,ym, &
547: & PETSC_NULL_INTEGER,ierr)
548: call DMDAGetGhostCorners(da,gxs,gys,PETSC_NULL_INTEGER,gxm,gym, &
549: & PETSC_NULL_INTEGER,ierr)
551: ! Get the global node numbers for all local nodes, including ghost points
553: call DMGetLocalToGlobalMapping(da,ltogm,ierr)
554: call ISLocalToGlobalMappingGetIndices(ltogm,ltog,idltog,ierr)
556: ! Compute entries for the locally owned part of the Jacobian.
557: ! - Currently, all PETSc parallel matrix formats are partitioned by
558: ! contiguous chunks of rows across the processors. The 'grow'
559: ! parameter computed below specifies the global row number
560: ! corresponding to each local grid point.
561: ! - Each processor needs to insert only elements that it owns
562: ! locally (but any non-local elements will be sent to the
563: ! appropriate processor during matrix assembly).
564: ! - Always specify global row and columns of matrix entries.
565: ! - Here, we set all entries for a particular row at once.
567: do 10 j=ys,ys+ym-1
568: row = (j - gys)*gxm + xs - gxs
569: do 20 i=xs,xs+xm-1
570: row = row + 1
571: grow(1) = ltog(idltog+row)
572: if (i .eq. 0 .or. j .eq. 0 .or. i .eq. (mx-1) .or. &
573: & j .eq. (my-1)) then
574: call MatSetValues(jac,ione,grow,ione,grow,one, &
575: & INSERT_VALUES,ierr)
576: go to 20
577: endif
578: v(1) = -hxdhy
579: col(1) = ltog(idltog+row - gxm)
580: v(2) = -hydhx
581: col(2) = ltog(idltog+row - 1)
582: v(3) = two*(hydhx + hxdhy) - sc*lambda*exp(xx(idx+row))
583: col(3) = grow(1)
584: v(4) = -hydhx
585: col(4) = ltog(idltog+row + 1)
586: v(5) = -hxdhy
587: col(5) = ltog(idltog+row + gxm)
588: call MatSetValues(jac,ione,grow,ifive,col,v,INSERT_VALUES, &
589: & ierr)
590: 20 continue
591: 10 continue
593: call ISLocalToGlobalMappingRestoreIndices(ltogm,ltog,idltog,ierr)
595: ! Assemble matrix, using the 2-step process:
596: ! MatAssemblyBegin(), MatAssemblyEnd().
597: ! By placing code between these two statements, computations can be
598: ! done while messages are in transition.
600: call MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY,ierr)
601: call VecRestoreArray(localX,xx,idx,ierr)
602: call MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY,ierr)
603: return
604: end
607: ! -------------------------------------------------------------------
608: !
609: ! MyMult - user provided matrix multiply
610: !
611: ! Input Parameters:
612: !. X - input vector
613: !
614: ! Output Parameter:
615: !. F - function vector
616: !
617: subroutine MyMult(J,X,F,ierr)
618: implicit none
619: Mat J,B
620: Vec X,F
621: PetscErrorCode ierr
622: PetscInt mx,my
623: DM da
624: Vec localX
626: common /mycommon/ mx,my,B,localX,da
627: !
628: ! Here we use the actual formed matrix B; users would
629: ! instead write their own matrix vector product routine
630: !
631: call MatMult(B,X,F,ierr)
632: return
633: end