xref: /petsc/src/mat/impls/shell/shell.c (revision 92fabd1f57b177038eb4af3534af2066536c5d0d) 
1 
2 /*
3    This provides a simple shell for Fortran (and C programmers) to
4   create a very simple matrix class for use with KSP without coding
5   much of anything.
6 */
7 
8 #include <petsc/private/matimpl.h>        /*I "petscmat.h" I*/
9 
10 struct _MatShellOps {
11   /*   0 */
12   PetscErrorCode (*mult)(Mat,Vec,Vec);
13   /*   5 */
14   PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
15   /*  10 */
16   /*  15 */
17   PetscErrorCode (*getdiagonal)(Mat,Vec);
18   /*  20 */
19   /*  24 */
20   /*  29 */
21   /*  34 */
22   /*  39 */
23   PetscErrorCode (*axpy)(Mat,PetscScalar,Mat,MatStructure);
24   PetscErrorCode (*copy)(Mat,Mat,MatStructure);
25   /*  44 */
26   PetscErrorCode (*diagonalset)(Mat,Vec,InsertMode);
27   /*  49 */
28   /*  54 */
29   /*  59 */
30   PetscErrorCode (*destroy)(Mat);
31   /*  64 */
32   /*  69 */
33   /*  74 */
34   /*  79 */
35   /*  84 */
36   /*  89 */
37   /*  94 */
38   /*  99 */
39   /* 104 */
40   /* 109 */
41   /* 114 */
42   /* 119 */
43   /* 124 */
44   /* 129 */
45   /* 134 */
46   /* 139 */
47   /* 144 */
48 };
49 
50 typedef struct {
51   struct _MatShellOps ops[1];
52 
53   PetscScalar vscale,vshift;
54   Vec         dshift;
55   Vec         left,right;
56   Vec         left_work,right_work;
57   Vec         left_add_work,right_add_work;
58   Mat         axpy;
59   PetscScalar axpy_vscale;
60   PetscBool   managescalingshifts;                   /* The user will manage the scaling and shifts for the MATSHELL, not the default */
61   void        *ctx;
62 } Mat_Shell;
63 
64 /*
65       xx = diag(left)*x
66 */
67 static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
68 {
69   Mat_Shell      *shell = (Mat_Shell*)A->data;
70   PetscErrorCode ierr;
71 
72   PetscFunctionBegin;
73   *xx = NULL;
74   if (!shell->left) {
75     *xx = x;
76   } else {
77     if (!shell->left_work) {ierr = VecDuplicate(shell->left,&shell->left_work);CHKERRQ(ierr);}
78     ierr = VecPointwiseMult(shell->left_work,x,shell->left);CHKERRQ(ierr);
79     *xx  = shell->left_work;
80   }
81   PetscFunctionReturn(0);
82 }
83 
84 /*
85      xx = diag(right)*x
86 */
87 static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
88 {
89   Mat_Shell      *shell = (Mat_Shell*)A->data;
90   PetscErrorCode ierr;
91 
92   PetscFunctionBegin;
93   *xx = NULL;
94   if (!shell->right) {
95     *xx = x;
96   } else {
97     if (!shell->right_work) {ierr = VecDuplicate(shell->right,&shell->right_work);CHKERRQ(ierr);}
98     ierr = VecPointwiseMult(shell->right_work,x,shell->right);CHKERRQ(ierr);
99     *xx  = shell->right_work;
100   }
101   PetscFunctionReturn(0);
102 }
103 
104 /*
105     x = diag(left)*x
106 */
107 static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
108 {
109   Mat_Shell      *shell = (Mat_Shell*)A->data;
110   PetscErrorCode ierr;
111 
112   PetscFunctionBegin;
113   if (shell->left) {ierr = VecPointwiseMult(x,x,shell->left);CHKERRQ(ierr);}
114   PetscFunctionReturn(0);
115 }
116 
117 /*
118     x = diag(right)*x
119 */
120 static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
121 {
122   Mat_Shell      *shell = (Mat_Shell*)A->data;
123   PetscErrorCode ierr;
124 
125   PetscFunctionBegin;
126   if (shell->right) {ierr = VecPointwiseMult(x,x,shell->right);CHKERRQ(ierr);}
127   PetscFunctionReturn(0);
128 }
129 
130 /*
131          Y = vscale*Y + diag(dshift)*X + vshift*X
132 
133          On input Y already contains A*x
134 */
135 static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
136 {
137   Mat_Shell      *shell = (Mat_Shell*)A->data;
138   PetscErrorCode ierr;
139 
140   PetscFunctionBegin;
141   if (shell->dshift) {          /* get arrays because there is no VecPointwiseMultAdd() */
142     PetscInt          i,m;
143     const PetscScalar *x,*d;
144     PetscScalar       *y;
145     ierr = VecGetLocalSize(X,&m);CHKERRQ(ierr);
146     ierr = VecGetArrayRead(shell->dshift,&d);CHKERRQ(ierr);
147     ierr = VecGetArrayRead(X,&x);CHKERRQ(ierr);
148     ierr = VecGetArray(Y,&y);CHKERRQ(ierr);
149     for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
150     ierr = VecRestoreArrayRead(shell->dshift,&d);CHKERRQ(ierr);
151     ierr = VecRestoreArrayRead(X,&x);CHKERRQ(ierr);
152     ierr = VecRestoreArray(Y,&y);CHKERRQ(ierr);
153   } else {
154     ierr = VecScale(Y,shell->vscale);CHKERRQ(ierr);
155   }
156   if (shell->vshift != 0.0) {ierr = VecAXPY(Y,shell->vshift,X);CHKERRQ(ierr);} /* if test is for non-square matrices */
157   PetscFunctionReturn(0);
158 }
159 
160 /*@
161     MatShellGetContext - Returns the user-provided context associated with a shell matrix.
162 
163     Not Collective
164 
165     Input Parameter:
166 .   mat - the matrix, should have been created with MatCreateShell()
167 
168     Output Parameter:
169 .   ctx - the user provided context
170 
171     Level: advanced
172 
173    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
174     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
175 
176 .keywords: matrix, shell, get, context
177 
178 .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
179 @*/
180 PetscErrorCode  MatShellGetContext(Mat mat,void *ctx)
181 {
182   PetscErrorCode ierr;
183   PetscBool      flg;
184 
185   PetscFunctionBegin;
186   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
187   PetscValidPointer(ctx,2);
188   ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr);
189   if (flg) *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
190   else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot get context from non-shell matrix");
191   PetscFunctionReturn(0);
192 }
193 
194 PetscErrorCode MatDestroy_Shell(Mat mat)
195 {
196   PetscErrorCode ierr;
197   Mat_Shell      *shell = (Mat_Shell*)mat->data;
198 
199   PetscFunctionBegin;
200   if (shell->ops->destroy) {
201     ierr = (*shell->ops->destroy)(mat);CHKERRQ(ierr);
202   }
203   ierr = VecDestroy(&shell->left);CHKERRQ(ierr);
204   ierr = VecDestroy(&shell->right);CHKERRQ(ierr);
205   ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr);
206   ierr = VecDestroy(&shell->left_work);CHKERRQ(ierr);
207   ierr = VecDestroy(&shell->right_work);CHKERRQ(ierr);
208   ierr = VecDestroy(&shell->left_add_work);CHKERRQ(ierr);
209   ierr = VecDestroy(&shell->right_add_work);CHKERRQ(ierr);
210   ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr);
211   ierr = PetscFree(mat->data);CHKERRQ(ierr);
212   PetscFunctionReturn(0);
213 }
214 
215 PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
216 {
217   Mat_Shell       *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
218   PetscErrorCode  ierr;
219   PetscBool       matflg;
220 
221   PetscFunctionBegin;
222   ierr = PetscObjectTypeCompare((PetscObject)B,MATSHELL,&matflg);CHKERRQ(ierr);
223   if (!matflg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_NOTSAMETYPE,"Output matrix must be a MATSHELL");
224 
225   ierr = PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));CHKERRQ(ierr);
226   ierr = PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));CHKERRQ(ierr);
227 
228   if (shellA->ops->copy) {
229     ierr = (*shellA->ops->copy)(A,B,str);CHKERRQ(ierr);
230   }
231   shellB->vscale = shellA->vscale;
232   shellB->vshift = shellA->vshift;
233   if (shellA->dshift) {
234     if (!shellB->dshift) {
235       ierr = VecDuplicate(shellA->dshift,&shellB->dshift);CHKERRQ(ierr);
236     }
237     ierr = VecCopy(shellA->dshift,shellB->dshift);CHKERRQ(ierr);
238   } else {
239     ierr = VecDestroy(&shellB->dshift);CHKERRQ(ierr);
240   }
241   if (shellA->left) {
242     if (!shellB->left) {
243       ierr = VecDuplicate(shellA->left,&shellB->left);CHKERRQ(ierr);
244     }
245     ierr = VecCopy(shellA->left,shellB->left);CHKERRQ(ierr);
246   } else {
247     ierr = VecDestroy(&shellB->left);CHKERRQ(ierr);
248   }
249   if (shellA->right) {
250     if (!shellB->right) {
251       ierr = VecDuplicate(shellA->right,&shellB->right);CHKERRQ(ierr);
252     }
253     ierr = VecCopy(shellA->right,shellB->right);CHKERRQ(ierr);
254   } else {
255     ierr = VecDestroy(&shellB->right);CHKERRQ(ierr);
256   }
257   PetscFunctionReturn(0);
258 }
259 
260 PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
261 {
262   PetscErrorCode ierr;
263   void           *ctx;
264 
265   PetscFunctionBegin;
266   ierr = MatShellGetContext(mat,&ctx);CHKERRQ(ierr);
267   ierr = MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);CHKERRQ(ierr);
268   ierr = MatCopy(mat,*M,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
269   PetscFunctionReturn(0);
270 }
271 
272 PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
273 {
274   Mat_Shell        *shell = (Mat_Shell*)A->data;
275   PetscErrorCode   ierr;
276   Vec              xx;
277   PetscObjectState instate,outstate;
278 
279   PetscFunctionBegin;
280   ierr = MatShellPreScaleRight(A,x,&xx);CHKERRQ(ierr);
281   ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr);
282   ierr = (*shell->ops->mult)(A,xx,y);CHKERRQ(ierr);
283   ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr);
284   if (instate == outstate) {
285     /* increase the state of the output vector since the user did not update its state themself as should have been done */
286     ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr);
287   }
288   ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr);
289   ierr = MatShellPostScaleLeft(A,y);CHKERRQ(ierr);
290 
291   if (shell->axpy) {
292     if (!shell->left_work) {ierr = MatCreateVecs(A,&shell->left_work,NULL);CHKERRQ(ierr);}
293     ierr = MatMult(shell->axpy,x,shell->left_work);CHKERRQ(ierr);
294     ierr = VecAXPY(y,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr);
295   }
296   PetscFunctionReturn(0);
297 }
298 
299 PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
300 {
301   Mat_Shell      *shell = (Mat_Shell*)A->data;
302   PetscErrorCode ierr;
303 
304   PetscFunctionBegin;
305   if (y == z) {
306     if (!shell->right_add_work) {ierr = VecDuplicate(z,&shell->right_add_work);CHKERRQ(ierr);}
307     ierr = MatMult(A,x,shell->right_add_work);CHKERRQ(ierr);
308     ierr = VecAXPY(z,1.0,shell->right_add_work);CHKERRQ(ierr);
309   } else {
310     ierr = MatMult(A,x,z);CHKERRQ(ierr);
311     ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr);
312   }
313   PetscFunctionReturn(0);
314 }
315 
316 PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
317 {
318   Mat_Shell        *shell = (Mat_Shell*)A->data;
319   PetscErrorCode   ierr;
320   Vec              xx;
321   PetscObjectState instate,outstate;
322 
323   PetscFunctionBegin;
324   ierr = MatShellPreScaleLeft(A,x,&xx);CHKERRQ(ierr);
325   ierr = PetscObjectStateGet((PetscObject)y, &instate);CHKERRQ(ierr);
326   if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
327   ierr = (*shell->ops->multtranspose)(A,xx,y);CHKERRQ(ierr);
328   ierr = PetscObjectStateGet((PetscObject)y, &outstate);CHKERRQ(ierr);
329   if (instate == outstate) {
330     /* increase the state of the output vector since the user did not update its state themself as should have been done */
331     ierr = PetscObjectStateIncrease((PetscObject)y);CHKERRQ(ierr);
332   }
333   ierr = MatShellShiftAndScale(A,xx,y);CHKERRQ(ierr);
334   ierr = MatShellPostScaleRight(A,y);CHKERRQ(ierr);
335   PetscFunctionReturn(0);
336 }
337 
338 PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
339 {
340   Mat_Shell      *shell = (Mat_Shell*)A->data;
341   PetscErrorCode ierr;
342 
343   PetscFunctionBegin;
344   if (y == z) {
345     if (!shell->left_add_work) {ierr = VecDuplicate(z,&shell->left_add_work);CHKERRQ(ierr);}
346     ierr = MatMultTranspose(A,x,shell->left_add_work);CHKERRQ(ierr);
347     ierr = VecWAXPY(z,1.0,shell->left_add_work,y);CHKERRQ(ierr);
348   } else {
349     ierr = MatMultTranspose(A,x,z);CHKERRQ(ierr);
350     ierr = VecAXPY(z,1.0,y);CHKERRQ(ierr);
351   }
352   PetscFunctionReturn(0);
353 }
354 
355 /*
356           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
357 */
358 PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
359 {
360   Mat_Shell      *shell = (Mat_Shell*)A->data;
361   PetscErrorCode ierr;
362 
363   PetscFunctionBegin;
364   if (shell->ops->getdiagonal) {
365     ierr = (*shell->ops->getdiagonal)(A,v);CHKERRQ(ierr);
366   } else {
367     ierr = VecSet(v,0.0);CHKERRQ(ierr);
368   }
369   ierr = VecScale(v,shell->vscale);CHKERRQ(ierr);
370   if (shell->dshift) {
371     ierr = VecAXPY(v,1.0,shell->dshift);CHKERRQ(ierr);
372   }
373   ierr = VecShift(v,shell->vshift);CHKERRQ(ierr);
374   if (shell->left)  {ierr = VecPointwiseMult(v,v,shell->left);CHKERRQ(ierr);}
375   if (shell->right) {ierr = VecPointwiseMult(v,v,shell->right);CHKERRQ(ierr);}
376   if (shell->axpy) {
377     if (!shell->left_work) {ierr = VecDuplicate(v,&shell->left_work);CHKERRQ(ierr);}
378     ierr = MatGetDiagonal(shell->axpy,shell->left_work);CHKERRQ(ierr);
379     ierr = VecAXPY(v,shell->axpy_vscale,shell->left_work);CHKERRQ(ierr);
380   }
381   PetscFunctionReturn(0);
382 }
383 
384 PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
385 {
386   Mat_Shell      *shell = (Mat_Shell*)Y->data;
387   PetscErrorCode ierr;
388 
389   PetscFunctionBegin;
390   if (shell->left || shell->right) {
391     if (!shell->dshift) {
392       ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr);
393       ierr = VecSet(shell->dshift,a);CHKERRQ(ierr);
394     } else {
395       if (shell->left)  {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);}
396       if (shell->right) {ierr = VecPointwiseMult(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);}
397       ierr = VecShift(shell->dshift,a);CHKERRQ(ierr);
398     }
399     if (shell->left)  {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);CHKERRQ(ierr);}
400     if (shell->right) {ierr = VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);CHKERRQ(ierr);}
401   } else shell->vshift += a;
402   PetscFunctionReturn(0);
403 }
404 
405 PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
406 {
407   Mat_Shell      *shell = (Mat_Shell*)A->data;
408   PetscErrorCode ierr;
409 
410   PetscFunctionBegin;
411   if (ins == INSERT_VALUES) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE, "Operation not supported with INSERT_VALUES");
412   if (!shell->dshift) {ierr = VecDuplicate(D,&shell->dshift);CHKERRQ(ierr);}
413   if (shell->left || shell->right) {
414     if (!shell->right_work) {ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);CHKERRQ(ierr);}
415     if (shell->left && shell->right)  {
416       ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr);
417       ierr = VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);CHKERRQ(ierr);
418     } else if (shell->left) {
419       ierr = VecPointwiseDivide(shell->right_work,D,shell->left);CHKERRQ(ierr);
420     } else {
421       ierr = VecPointwiseDivide(shell->right_work,D,shell->right);CHKERRQ(ierr);
422     }
423     if (!shell->dshift) {
424       ierr = VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);CHKERRQ(ierr);
425       ierr = VecCopy(shell->dshift,shell->right_work);CHKERRQ(ierr);
426     } else {
427       ierr = VecAXPY(shell->dshift,1.0,shell->right_work);CHKERRQ(ierr);
428     }
429   } else {
430     ierr = VecAXPY(shell->dshift,1.0,D);CHKERRQ(ierr);
431   }
432   PetscFunctionReturn(0);
433 }
434 
435 PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
436 {
437   Mat_Shell      *shell = (Mat_Shell*)Y->data;
438   PetscErrorCode ierr;
439 
440   PetscFunctionBegin;
441   shell->vscale *= a;
442   shell->vshift *= a;
443   if (shell->dshift) {
444     ierr = VecScale(shell->dshift,a);CHKERRQ(ierr);
445   }
446   shell->axpy_vscale *= a;
447   PetscFunctionReturn(0);
448 }
449 
450 static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
451 {
452   Mat_Shell      *shell = (Mat_Shell*)Y->data;
453   PetscErrorCode ierr;
454 
455   PetscFunctionBegin;
456   if (shell->axpy) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot diagonal scale MATSHELL after MatAXPY operation");
457   if (left) {
458     if (!shell->left) {
459       ierr = VecDuplicate(left,&shell->left);CHKERRQ(ierr);
460       ierr = VecCopy(left,shell->left);CHKERRQ(ierr);
461     } else {
462       ierr = VecPointwiseMult(shell->left,shell->left,left);CHKERRQ(ierr);
463     }
464   }
465   if (right) {
466     if (!shell->right) {
467       ierr = VecDuplicate(right,&shell->right);CHKERRQ(ierr);
468       ierr = VecCopy(right,shell->right);CHKERRQ(ierr);
469     } else {
470       ierr = VecPointwiseMult(shell->right,shell->right,right);CHKERRQ(ierr);
471     }
472   }
473   PetscFunctionReturn(0);
474 }
475 
476 PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
477 {
478   Mat_Shell      *shell = (Mat_Shell*)Y->data;
479   PetscErrorCode ierr;
480 
481   PetscFunctionBegin;
482   if (t == MAT_FINAL_ASSEMBLY) {
483     shell->vshift = 0.0;
484     shell->vscale = 1.0;
485     ierr = VecDestroy(&shell->dshift);CHKERRQ(ierr);
486     ierr = VecDestroy(&shell->left);CHKERRQ(ierr);
487     ierr = VecDestroy(&shell->right);CHKERRQ(ierr);
488     ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr);
489   }
490   PetscFunctionReturn(0);
491 }
492 
493 PETSC_INTERN PetscErrorCode MatConvert_Shell(Mat, MatType,MatReuse,Mat*);
494 
495 static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool  *missing,PetscInt *d)
496 {
497   PetscFunctionBegin;
498   *missing = PETSC_FALSE;
499   PetscFunctionReturn(0);
500 }
501 
502 PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
503 {
504   Mat_Shell      *shell = (Mat_Shell*)Y->data;
505   PetscErrorCode ierr;
506 
507   PetscFunctionBegin;
508   ierr = PetscObjectReference((PetscObject)X);CHKERRQ(ierr);
509   ierr = MatDestroy(&shell->axpy);CHKERRQ(ierr);
510   shell->axpy        = X;
511   shell->axpy_vscale = a;
512   PetscFunctionReturn(0);
513 }
514 
515 static struct _MatOps MatOps_Values = {0,
516                                        0,
517                                        0,
518                                        0,
519                                 /* 4*/ MatMultAdd_Shell,
520                                        0,
521                                        MatMultTransposeAdd_Shell,
522                                        0,
523                                        0,
524                                        0,
525                                 /*10*/ 0,
526                                        0,
527                                        0,
528                                        0,
529                                        0,
530                                 /*15*/ 0,
531                                        0,
532                                        MatGetDiagonal_Shell,
533                                        MatDiagonalScale_Shell,
534                                        0,
535                                 /*20*/ 0,
536                                        MatAssemblyEnd_Shell,
537                                        0,
538                                        0,
539                                 /*24*/ 0,
540                                        0,
541                                        0,
542                                        0,
543                                        0,
544                                 /*29*/ 0,
545                                        0,
546                                        0,
547                                        0,
548                                        0,
549                                 /*34*/ MatDuplicate_Shell,
550                                        0,
551                                        0,
552                                        0,
553                                        0,
554                                 /*39*/ MatAXPY_Shell,
555                                        0,
556                                        0,
557                                        0,
558                                        MatCopy_Shell,
559                                 /*44*/ 0,
560                                        MatScale_Shell,
561                                        MatShift_Shell,
562                                        MatDiagonalSet_Shell,
563                                        0,
564                                 /*49*/ 0,
565                                        0,
566                                        0,
567                                        0,
568                                        0,
569                                 /*54*/ 0,
570                                        0,
571                                        0,
572                                        0,
573                                        0,
574                                 /*59*/ 0,
575                                        MatDestroy_Shell,
576                                        0,
577                                        0,
578                                        0,
579                                 /*64*/ 0,
580                                        0,
581                                        0,
582                                        0,
583                                        0,
584                                 /*69*/ 0,
585                                        0,
586                                        MatConvert_Shell,
587                                        0,
588                                        0,
589                                 /*74*/ 0,
590                                        0,
591                                        0,
592                                        0,
593                                        0,
594                                 /*79*/ 0,
595                                        0,
596                                        0,
597                                        0,
598                                        0,
599                                 /*84*/ 0,
600                                        0,
601                                        0,
602                                        0,
603                                        0,
604                                 /*89*/ 0,
605                                        0,
606                                        0,
607                                        0,
608                                        0,
609                                 /*94*/ 0,
610                                        0,
611                                        0,
612                                        0,
613                                        0,
614                                 /*99*/ 0,
615                                        0,
616                                        0,
617                                        0,
618                                        0,
619                                /*104*/ 0,
620                                        0,
621                                        0,
622                                        0,
623                                        0,
624                                /*109*/ 0,
625                                        0,
626                                        0,
627                                        0,
628                                        MatMissingDiagonal_Shell,
629                                /*114*/ 0,
630                                        0,
631                                        0,
632                                        0,
633                                        0,
634                                /*119*/ 0,
635                                        0,
636                                        0,
637                                        0,
638                                        0,
639                                /*124*/ 0,
640                                        0,
641                                        0,
642                                        0,
643                                        0,
644                                /*129*/ 0,
645                                        0,
646                                        0,
647                                        0,
648                                        0,
649                                /*134*/ 0,
650                                        0,
651                                        0,
652                                        0,
653                                        0,
654                                /*139*/ 0,
655                                        0,
656                                        0
657 };
658 
659 /*MC
660    MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
661 
662   Level: advanced
663 
664 .seealso: MatCreateShell()
665 M*/
666 
667 PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
668 {
669   Mat_Shell      *b;
670   PetscErrorCode ierr;
671 
672   PetscFunctionBegin;
673   ierr = PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));CHKERRQ(ierr);
674 
675   ierr    = PetscNewLog(A,&b);CHKERRQ(ierr);
676   A->data = (void*)b;
677 
678   ierr = PetscLayoutSetUp(A->rmap);CHKERRQ(ierr);
679   ierr = PetscLayoutSetUp(A->cmap);CHKERRQ(ierr);
680 
681   b->ctx                 = 0;
682   b->vshift              = 0.0;
683   b->vscale              = 1.0;
684   b->managescalingshifts = PETSC_TRUE;
685   A->assembled           = PETSC_TRUE;
686   A->preallocated        = PETSC_FALSE;
687 
688   ierr = PetscObjectChangeTypeName((PetscObject)A,MATSHELL);CHKERRQ(ierr);
689   PetscFunctionReturn(0);
690 }
691 
692 /*@C
693    MatCreateShell - Creates a new matrix class for use with a user-defined
694    private data storage format.
695 
696   Collective on MPI_Comm
697 
698    Input Parameters:
699 +  comm - MPI communicator
700 .  m - number of local rows (must be given)
701 .  n - number of local columns (must be given)
702 .  M - number of global rows (may be PETSC_DETERMINE)
703 .  N - number of global columns (may be PETSC_DETERMINE)
704 -  ctx - pointer to data needed by the shell matrix routines
705 
706    Output Parameter:
707 .  A - the matrix
708 
709    Level: advanced
710 
711   Usage:
712 $    extern int mult(Mat,Vec,Vec);
713 $    MatCreateShell(comm,m,n,M,N,ctx,&mat);
714 $    MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
715 $    [ Use matrix for operations that have been set ]
716 $    MatDestroy(mat);
717 
718    Notes:
719    The shell matrix type is intended to provide a simple class to use
720    with KSP (such as, for use with matrix-free methods). You should not
721    use the shell type if you plan to define a complete matrix class.
722 
723    Fortran Notes: To use this from Fortran with a ctx you must write an interface definition for this
724     function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
725     in as the ctx argument.
726 
727    PETSc requires that matrices and vectors being used for certain
728    operations are partitioned accordingly.  For example, when
729    creating a shell matrix, A, that supports parallel matrix-vector
730    products using MatMult(A,x,y) the user should set the number
731    of local matrix rows to be the number of local elements of the
732    corresponding result vector, y. Note that this is information is
733    required for use of the matrix interface routines, even though
734    the shell matrix may not actually be physically partitioned.
735    For example,
736 
737    MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), and MatScale() internally so these
738    operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
739 
740 $
741 $     Vec x, y
742 $     extern int mult(Mat,Vec,Vec);
743 $     Mat A
744 $
745 $     VecCreateMPI(comm,PETSC_DECIDE,M,&y);
746 $     VecCreateMPI(comm,PETSC_DECIDE,N,&x);
747 $     VecGetLocalSize(y,&m);
748 $     VecGetLocalSize(x,&n);
749 $     MatCreateShell(comm,m,n,M,N,ctx,&A);
750 $     MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
751 $     MatMult(A,x,y);
752 $     MatDestroy(A);
753 $     VecDestroy(y); VecDestroy(x);
754 $
755 
756     For rectangular matrices do all the scalings and shifts make sense?
757 
758     Developers Notes: Regarding shifting and scaling. The general form is
759 
760           diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
761 
762       The order you apply the operations is important. For example if you have a dshift then
763       apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
764       you get s*vscale*A + diag(shift)
765 
766           A is the user provided function.
767 
768 .keywords: matrix, shell, create
769 
770 .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts()
771 @*/
772 PetscErrorCode  MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
773 {
774   PetscErrorCode ierr;
775 
776   PetscFunctionBegin;
777   ierr = MatCreate(comm,A);CHKERRQ(ierr);
778   ierr = MatSetSizes(*A,m,n,M,N);CHKERRQ(ierr);
779   ierr = MatSetType(*A,MATSHELL);CHKERRQ(ierr);
780   ierr = MatShellSetContext(*A,ctx);CHKERRQ(ierr);
781   ierr = MatSetUp(*A);CHKERRQ(ierr);
782   PetscFunctionReturn(0);
783 }
784 
785 /*@
786     MatShellSetContext - sets the context for a shell matrix
787 
788    Logically Collective on Mat
789 
790     Input Parameters:
791 +   mat - the shell matrix
792 -   ctx - the context
793 
794    Level: advanced
795 
796    Fortran Notes: To use this from Fortran you must write a Fortran interface definition for this
797     function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
798 
799 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
800 @*/
801 PetscErrorCode  MatShellSetContext(Mat mat,void *ctx)
802 {
803   Mat_Shell      *shell = (Mat_Shell*)mat->data;
804   PetscErrorCode ierr;
805   PetscBool      flg;
806 
807   PetscFunctionBegin;
808   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
809   ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr);
810   if (flg) {
811     shell->ctx = ctx;
812   } else SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Cannot attach context to non-shell matrix");
813   PetscFunctionReturn(0);
814 }
815 
816 /*@
817     MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
818           after MatCreateShell()
819 
820    Logically Collective on Mat
821 
822     Input Parameter:
823 .   mat - the shell matrix
824 
825   Level: advanced
826 
827 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
828 @*/
829 PetscErrorCode MatShellSetManageScalingShifts(Mat A)
830 {
831   PetscErrorCode ierr;
832   Mat_Shell      *shell = (Mat_Shell*)A->data;
833   PetscBool      flg;
834 
835   PetscFunctionBegin;
836   PetscValidHeaderSpecific(A,MAT_CLASSID,1);
837   ierr = PetscObjectTypeCompare((PetscObject)A,MATSHELL,&flg);CHKERRQ(ierr);
838   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
839   shell->managescalingshifts = PETSC_FALSE;
840   A->ops->diagonalscale = 0;
841   A->ops->scale         = 0;
842   A->ops->shift         = 0;
843   A->ops->diagonalset   = 0;
844   PetscFunctionReturn(0);
845 }
846 
847 /*@C
848     MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
849 
850    Logically Collective on Mat
851 
852     Input Parameters:
853 +   mat - the shell matrix
854 .   f - the function
855 .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
856 -   ctx - an optional context for the function
857 
858    Output Parameter:
859 .   flg - PETSC_TRUE if the multiply is likely correct
860 
861    Options Database:
862 .   -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
863 
864    Level: advanced
865 
866    Fortran Notes: Not supported from Fortran
867 
868 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
869 @*/
870 PetscErrorCode  MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
871 {
872   PetscErrorCode ierr;
873   PetscInt       m,n;
874   Mat            mf,Dmf,Dmat,Ddiff;
875   PetscReal      Diffnorm,Dmfnorm;
876   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;
877 
878   PetscFunctionBegin;
879   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
880   ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);CHKERRQ(ierr);
881   ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr);
882   ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);CHKERRQ(ierr);
883   ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr);
884   ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr);
885 
886   ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr);
887   ierr = MatComputeExplicitOperator(mat,&Dmat);CHKERRQ(ierr);
888 
889   ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr);
890   ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
891   ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr);
892   ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr);
893   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
894     flag = PETSC_FALSE;
895     if (v) {
896       ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));
897       ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr);
898       ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr);
899       ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");CHKERRQ(ierr);
900     }
901   } else if (v) {
902     ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
903   }
904   if (flg) *flg = flag;
905   ierr = MatDestroy(&Ddiff);CHKERRQ(ierr);
906   ierr = MatDestroy(&mf);CHKERRQ(ierr);
907   ierr = MatDestroy(&Dmf);CHKERRQ(ierr);
908   ierr = MatDestroy(&Dmat);CHKERRQ(ierr);
909   PetscFunctionReturn(0);
910 }
911 
912 /*@C
913     MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
914 
915    Logically Collective on Mat
916 
917     Input Parameters:
918 +   mat - the shell matrix
919 .   f - the function
920 .   base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
921 -   ctx - an optional context for the function
922 
923    Output Parameter:
924 .   flg - PETSC_TRUE if the multiply is likely correct
925 
926    Options Database:
927 .   -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
928 
929    Level: advanced
930 
931    Fortran Notes: Not supported from Fortran
932 
933 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
934 @*/
935 PetscErrorCode  MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
936 {
937   PetscErrorCode ierr;
938   Vec            x,y,z;
939   PetscInt       m,n,M,N;
940   Mat            mf,Dmf,Dmat,Ddiff;
941   PetscReal      Diffnorm,Dmfnorm;
942   PetscBool      v = PETSC_FALSE, flag = PETSC_TRUE;
943 
944   PetscFunctionBegin;
945   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
946   ierr = PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);CHKERRQ(ierr);
947   ierr = MatCreateVecs(mat,&x,&y);CHKERRQ(ierr);
948   ierr = VecDuplicate(y,&z);CHKERRQ(ierr);
949   ierr = MatGetLocalSize(mat,&m,&n);CHKERRQ(ierr);
950   ierr = MatGetSize(mat,&M,&N);CHKERRQ(ierr);
951   ierr = MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);CHKERRQ(ierr);
952   ierr = MatMFFDSetFunction(mf,f,ctx);CHKERRQ(ierr);
953   ierr = MatMFFDSetBase(mf,base,NULL);CHKERRQ(ierr);
954   ierr = MatComputeExplicitOperator(mf,&Dmf);CHKERRQ(ierr);
955   ierr = MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);CHKERRQ(ierr);
956   ierr = MatComputeExplicitOperatorTranspose(mat,&Dmat);CHKERRQ(ierr);
957 
958   ierr = MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);CHKERRQ(ierr);
959   ierr = MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);CHKERRQ(ierr);
960   ierr = MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);CHKERRQ(ierr);
961   ierr = MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);CHKERRQ(ierr);
962   if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
963     flag = PETSC_FALSE;
964     if (v) {
965       ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce different results.\n Norm Ratio %g Difference results followed by finite difference one\n",(double)(Diffnorm/Dmfnorm));
966       ierr = MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr);
967       ierr = MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr);
968       ierr = MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");CHKERRQ(ierr);
969     }
970   } else if (v) {
971     ierr = PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
972   }
973   if (flg) *flg = flag;
974   ierr = MatDestroy(&mf);CHKERRQ(ierr);
975   ierr = MatDestroy(&Dmat);CHKERRQ(ierr);
976   ierr = MatDestroy(&Ddiff);CHKERRQ(ierr);
977   ierr = MatDestroy(&Dmf);CHKERRQ(ierr);
978   ierr = VecDestroy(&x);CHKERRQ(ierr);
979   ierr = VecDestroy(&y);CHKERRQ(ierr);
980   ierr = VecDestroy(&z);CHKERRQ(ierr);
981   PetscFunctionReturn(0);
982 }
983 
984 /*@C
985     MatShellSetOperation - Allows user to set a matrix operation for
986                            a shell matrix.
987     MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
988 
989    Logically Collective on Mat
990 
991     Input Parameters:
992 +   mat - the shell matrix
993 .   op - the name of the operation
994 -   f - the function that provides the operation.
995 
996    Level: advanced
997 
998     Usage:
999 $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1000 $      ierr = MatCreateShell(comm,m,n,M,N,ctx,&A);
1001 $      ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
1002 
1003     Notes:
1004     See the file include/petscmat.h for a complete list of matrix
1005     operations, which all have the form MATOP_<OPERATION>, where
1006     <OPERATION> is the name (in all capital letters) of the
1007     user interface routine (e.g., MatMult() -> MATOP_MULT).
1008 
1009     All user-provided functions (except for MATOP_DESTROY) should have the same calling
1010     sequence as the usual matrix interface routines, since they
1011     are intended to be accessed via the usual matrix interface
1012     routines, e.g.,
1013 $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1014 
1015     In particular each function MUST return an error code of 0 on success and
1016     nonzero on failure.
1017 
1018     Within each user-defined routine, the user should call
1019     MatShellGetContext() to obtain the user-defined context that was
1020     set by MatCreateShell().
1021 
1022     Fortran Notes: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
1023        generate a matrix. See src/mat/examples/tests/ex120f.F
1024 
1025     Use MatSetOperation() to set an operation for any matrix type
1026 
1027 .keywords: matrix, shell, set, operation
1028 
1029 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts()
1030 @*/
1031 PetscErrorCode  MatShellSetOperation(Mat mat,MatOperation op,void (*f)(void))
1032 {
1033   PetscErrorCode ierr;
1034   PetscBool      flg;
1035   Mat_Shell      *shell = (Mat_Shell*)mat->data;
1036 
1037   PetscFunctionBegin;
1038   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
1039   ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr);
1040   if (!flg) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_SUP,"Can only use with MATSHELL matrices");
1041   switch (op) {
1042   case MATOP_DESTROY:
1043     shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1044     break;
1045   case MATOP_DIAGONAL_SET:
1046   case MATOP_DIAGONAL_SCALE:
1047   case MATOP_SHIFT:
1048   case MATOP_SCALE:
1049   case MATOP_AXPY:
1050     if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1051     (((void(**)(void))mat->ops)[op]) = f;
1052     break;
1053   case MATOP_GET_DIAGONAL:
1054     if (shell->managescalingshifts) shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1055     else mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1056   case MATOP_VIEW:
1057     if (!mat->ops->viewnative) {
1058       mat->ops->viewnative = mat->ops->view;
1059     }
1060     mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1061     break;
1062   case MATOP_MULT:
1063     if (shell->managescalingshifts){
1064       shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1065       mat->ops->mult   = MatMult_Shell;
1066     } else mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1067     break;
1068   case MATOP_MULT_TRANSPOSE:
1069     if (shell->managescalingshifts) {
1070       shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1071       mat->ops->multtranspose   = MatMultTranspose_Shell;
1072     } else mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1073     break;
1074   default:
1075     (((void(**)(void))mat->ops)[op]) = f;
1076   }
1077   PetscFunctionReturn(0);
1078 }
1079 
1080 /*@C
1081     MatShellGetOperation - Gets a matrix function for a shell matrix.
1082 
1083     Not Collective
1084 
1085     Input Parameters:
1086 +   mat - the shell matrix
1087 -   op - the name of the operation
1088 
1089     Output Parameter:
1090 .   f - the function that provides the operation.
1091 
1092     Level: advanced
1093 
1094     Notes:
1095     See the file include/petscmat.h for a complete list of matrix
1096     operations, which all have the form MATOP_<OPERATION>, where
1097     <OPERATION> is the name (in all capital letters) of the
1098     user interface routine (e.g., MatMult() -> MATOP_MULT).
1099 
1100     All user-provided functions have the same calling
1101     sequence as the usual matrix interface routines, since they
1102     are intended to be accessed via the usual matrix interface
1103     routines, e.g.,
1104 $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1105 
1106     Within each user-defined routine, the user should call
1107     MatShellGetContext() to obtain the user-defined context that was
1108     set by MatCreateShell().
1109 
1110 .keywords: matrix, shell, set, operation
1111 
1112 .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
1113 @*/
1114 PetscErrorCode  MatShellGetOperation(Mat mat,MatOperation op,void(**f)(void))
1115 {
1116   PetscErrorCode ierr;
1117   PetscBool      flg;
1118 
1119   PetscFunctionBegin;
1120   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
1121   switch (op) {
1122   case MATOP_DESTROY:
1123     ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr);
1124     if (flg) {
1125       Mat_Shell *shell = (Mat_Shell*)mat->data;
1126       *f = (void (*)(void))shell->ops->destroy;
1127     } else {
1128       *f = (void (*)(void))mat->ops->destroy;
1129     }
1130     break;
1131   case MATOP_DIAGONAL_SET:
1132     ierr = PetscObjectTypeCompare((PetscObject)mat,MATSHELL,&flg);CHKERRQ(ierr);
1133     if (flg) {
1134       Mat_Shell *shell = (Mat_Shell*)mat->data;
1135       *f = (void (*)(void))shell->ops->diagonalset;
1136     } else {
1137       *f = (void (*)(void))mat->ops->diagonalset;
1138     }
1139     break;
1140   case MATOP_VIEW:
1141     *f = (void (*)(void))mat->ops->view;
1142     break;
1143   default:
1144     *f = (((void (**)(void))mat->ops)[op]);
1145   }
1146   PetscFunctionReturn(0);
1147 }
1148 
1149 /*@C
1150     MatSetOperation - Allows user to set a matrix operation for any matrix type
1151 
1152    Logically Collective on Mat
1153 
1154     Input Parameters:
1155 +   mat - the shell matrix
1156 .   op - the name of the operation
1157 -   f - the function that provides the operation.
1158 
1159    Level: developer
1160 
1161     Usage:
1162 $      extern PetscErrorCode usermult(Mat,Vec,Vec);
1163 $      ierr = MatCreateXXX(comm,...&A);
1164 $      ierr = MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
1165 
1166     Notes:
1167     See the file include/petscmat.h for a complete list of matrix
1168     operations, which all have the form MATOP_<OPERATION>, where
1169     <OPERATION> is the name (in all capital letters) of the
1170     user interface routine (e.g., MatMult() -> MATOP_MULT).
1171 
1172     All user-provided functions (except for MATOP_DESTROY) should have the same calling
1173     sequence as the usual matrix interface routines, since they
1174     are intended to be accessed via the usual matrix interface
1175     routines, e.g.,
1176 $       MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
1177 
1178     In particular each function MUST return an error code of 0 on success and
1179     nonzero on failure.
1180 
1181     This routine is distinct from MatShellSetOperation() in that it can be called on any matrix type.
1182 
1183 .keywords: matrix, shell, set, operation
1184 
1185 .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1186 @*/
1187 PetscErrorCode  MatSetOperation(Mat mat,MatOperation op,void (*f)(void))
1188 {
1189   PetscFunctionBegin;
1190   PetscValidHeaderSpecific(mat,MAT_CLASSID,1);
1191   (((void(**)(void))mat->ops)[op]) = f;
1192   PetscFunctionReturn(0);
1193 }
1194