Actual source code: shell.c
petsc-3.14.6 2021-03-30
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: */
8: #include <petsc/private/matimpl.h>
9: #include <petsc/private/vecimpl.h>
11: struct _MatShellOps {
12: /* 3 */ PetscErrorCode (*mult)(Mat,Vec,Vec);
13: /* 5 */ PetscErrorCode (*multtranspose)(Mat,Vec,Vec);
14: /* 17 */ PetscErrorCode (*getdiagonal)(Mat,Vec);
15: /* 43 */ PetscErrorCode (*copy)(Mat,Mat,MatStructure);
16: /* 60 */ PetscErrorCode (*destroy)(Mat);
17: };
19: struct _n_MatShellMatFunctionList {
20: PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**);
21: PetscErrorCode (*numeric)(Mat,Mat,Mat,void*);
22: PetscErrorCode (*destroy)(void*);
23: MatProductType ptype;
24: char *composedname; /* string to identify routine with double dispatch */
25: char *resultname; /* result matrix type */
27: struct _n_MatShellMatFunctionList *next;
28: };
29: typedef struct _n_MatShellMatFunctionList *MatShellMatFunctionList;
31: typedef struct {
32: struct _MatShellOps ops[1];
34: /* The user will manage the scaling and shifts for the MATSHELL, not the default */
35: PetscBool managescalingshifts;
37: /* support for MatScale, MatShift and MatMultAdd */
38: PetscScalar vscale,vshift;
39: Vec dshift;
40: Vec left,right;
41: Vec left_work,right_work;
42: Vec left_add_work,right_add_work;
44: /* support for MatAXPY */
45: Mat axpy;
46: PetscScalar axpy_vscale;
47: Vec axpy_left,axpy_right;
48: PetscObjectState axpy_state;
50: /* support for ZeroRows/Columns operations */
51: IS zrows;
52: IS zcols;
53: Vec zvals;
54: Vec zvals_w;
55: VecScatter zvals_sct_r;
56: VecScatter zvals_sct_c;
58: /* MatMat operations */
59: MatShellMatFunctionList matmat;
61: /* user defined context */
62: void *ctx;
63: } Mat_Shell;
66: /*
67: Store and scale values on zeroed rows
68: xx = [x_1, 0], 0 on zeroed columns
69: */
70: static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx)
71: {
72: Mat_Shell *shell = (Mat_Shell*)A->data;
76: *xx = x;
77: if (shell->zrows) {
78: VecSet(shell->zvals_w,0.0);
79: VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
80: VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
81: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
82: }
83: if (shell->zcols) {
84: if (!shell->right_work) {
85: MatCreateVecs(A,&shell->right_work,NULL);
86: }
87: VecCopy(x,shell->right_work);
88: VecISSet(shell->right_work,shell->zcols,0.0);
89: *xx = shell->right_work;
90: }
91: return(0);
92: }
94: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
95: static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x)
96: {
97: Mat_Shell *shell = (Mat_Shell*)A->data;
101: if (shell->zrows) {
102: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
103: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
104: }
105: return(0);
106: }
108: /*
109: Store and scale values on zeroed rows
110: xx = [x_1, 0], 0 on zeroed rows
111: */
112: static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx)
113: {
114: Mat_Shell *shell = (Mat_Shell*)A->data;
118: *xx = NULL;
119: if (!shell->zrows) {
120: *xx = x;
121: } else {
122: if (!shell->left_work) {
123: MatCreateVecs(A,NULL,&shell->left_work);
124: }
125: VecCopy(x,shell->left_work);
126: VecSet(shell->zvals_w,0.0);
127: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
128: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
129: VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
130: VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
131: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
132: *xx = shell->left_work;
133: }
134: return(0);
135: }
137: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
138: static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x)
139: {
140: Mat_Shell *shell = (Mat_Shell*)A->data;
144: if (shell->zcols) {
145: VecISSet(x,shell->zcols,0.0);
146: }
147: if (shell->zrows) {
148: VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
149: VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
150: }
151: return(0);
152: }
154: /*
155: xx = diag(left)*x
156: */
157: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
158: {
159: Mat_Shell *shell = (Mat_Shell*)A->data;
163: *xx = NULL;
164: if (!shell->left) {
165: *xx = x;
166: } else {
167: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
168: VecPointwiseMult(shell->left_work,x,shell->left);
169: *xx = shell->left_work;
170: }
171: return(0);
172: }
174: /*
175: xx = diag(right)*x
176: */
177: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
178: {
179: Mat_Shell *shell = (Mat_Shell*)A->data;
183: *xx = NULL;
184: if (!shell->right) {
185: *xx = x;
186: } else {
187: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
188: VecPointwiseMult(shell->right_work,x,shell->right);
189: *xx = shell->right_work;
190: }
191: return(0);
192: }
194: /*
195: x = diag(left)*x
196: */
197: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
198: {
199: Mat_Shell *shell = (Mat_Shell*)A->data;
203: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
204: return(0);
205: }
207: /*
208: x = diag(right)*x
209: */
210: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
211: {
212: Mat_Shell *shell = (Mat_Shell*)A->data;
216: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
217: return(0);
218: }
220: /*
221: Y = vscale*Y + diag(dshift)*X + vshift*X
223: On input Y already contains A*x
224: */
225: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
226: {
227: Mat_Shell *shell = (Mat_Shell*)A->data;
231: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
232: PetscInt i,m;
233: const PetscScalar *x,*d;
234: PetscScalar *y;
235: VecGetLocalSize(X,&m);
236: VecGetArrayRead(shell->dshift,&d);
237: VecGetArrayRead(X,&x);
238: VecGetArray(Y,&y);
239: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
240: VecRestoreArrayRead(shell->dshift,&d);
241: VecRestoreArrayRead(X,&x);
242: VecRestoreArray(Y,&y);
243: } else {
244: VecScale(Y,shell->vscale);
245: }
246: if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
247: return(0);
248: }
250: PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx)
251: {
253: *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
254: return(0);
255: }
257: /*@
258: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
260: Not Collective
262: Input Parameter:
263: . mat - the matrix, should have been created with MatCreateShell()
265: Output Parameter:
266: . ctx - the user provided context
268: Level: advanced
270: Fortran Notes:
271: To use this from Fortran you must write a Fortran interface definition for this
272: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
274: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
275: @*/
276: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
277: {
283: PetscUseMethod(mat,"MatShellGetContext_C",(Mat,void*),(mat,ctx));
284: return(0);
285: }
287: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc)
288: {
290: Mat_Shell *shell = (Mat_Shell*)mat->data;
291: Vec x = NULL,b = NULL;
292: IS is1, is2;
293: const PetscInt *ridxs;
294: PetscInt *idxs,*gidxs;
295: PetscInt cum,rst,cst,i;
298: if (!shell->zvals) {
299: MatCreateVecs(mat,NULL,&shell->zvals);
300: }
301: if (!shell->zvals_w) {
302: VecDuplicate(shell->zvals,&shell->zvals_w);
303: }
304: MatGetOwnershipRange(mat,&rst,NULL);
305: MatGetOwnershipRangeColumn(mat,&cst,NULL);
307: /* Expand/create index set of zeroed rows */
308: PetscMalloc1(nr,&idxs);
309: for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
310: ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);
311: ISSort(is1);
312: VecISSet(shell->zvals,is1,diag);
313: if (shell->zrows) {
314: ISSum(shell->zrows,is1,&is2);
315: ISDestroy(&shell->zrows);
316: ISDestroy(&is1);
317: shell->zrows = is2;
318: } else shell->zrows = is1;
320: /* Create scatters for diagonal values communications */
321: VecScatterDestroy(&shell->zvals_sct_c);
322: VecScatterDestroy(&shell->zvals_sct_r);
324: /* row scatter: from/to left vector */
325: MatCreateVecs(mat,&x,&b);
326: VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);
328: /* col scatter: from right vector to left vector */
329: ISGetIndices(shell->zrows,&ridxs);
330: ISGetLocalSize(shell->zrows,&nr);
331: PetscMalloc1(nr,&gidxs);
332: for (i = 0, cum = 0; i < nr; i++) {
333: if (ridxs[i] >= mat->cmap->N) continue;
334: gidxs[cum] = ridxs[i];
335: cum++;
336: }
337: ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);
338: VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);
339: ISDestroy(&is1);
340: VecDestroy(&x);
341: VecDestroy(&b);
343: /* Expand/create index set of zeroed columns */
344: if (rc) {
345: PetscMalloc1(nc,&idxs);
346: for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
347: ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);
348: ISSort(is1);
349: if (shell->zcols) {
350: ISSum(shell->zcols,is1,&is2);
351: ISDestroy(&shell->zcols);
352: ISDestroy(&is1);
353: shell->zcols = is2;
354: } else shell->zcols = is1;
355: }
356: return(0);
357: }
359: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
360: {
361: Mat_Shell *shell = (Mat_Shell*)mat->data;
362: PetscInt nr, *lrows;
366: if (x && b) {
367: Vec xt;
368: PetscScalar *vals;
369: PetscInt *gcols,i,st,nl,nc;
371: PetscMalloc1(n,&gcols);
372: for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
374: MatCreateVecs(mat,&xt,NULL);
375: VecCopy(x,xt);
376: PetscCalloc1(nc,&vals);
377: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
378: PetscFree(vals);
379: VecAssemblyBegin(xt);
380: VecAssemblyEnd(xt);
381: VecAYPX(xt,-1.0,x); /* xt = [0, x2] */
383: VecGetOwnershipRange(xt,&st,NULL);
384: VecGetLocalSize(xt,&nl);
385: VecGetArray(xt,&vals);
386: for (i = 0; i < nl; i++) {
387: PetscInt g = i + st;
388: if (g > mat->rmap->N) continue;
389: if (PetscAbsScalar(vals[i]) == 0.0) continue;
390: VecSetValue(b,g,diag*vals[i],INSERT_VALUES);
391: }
392: VecRestoreArray(xt,&vals);
393: VecAssemblyBegin(b);
394: VecAssemblyEnd(b); /* b = [b1, x2 * diag] */
395: VecDestroy(&xt);
396: PetscFree(gcols);
397: }
398: PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);
399: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);
400: if (shell->axpy) {
401: MatZeroRows(shell->axpy,n,rows,0.0,NULL,NULL);
402: }
403: PetscFree(lrows);
404: return(0);
405: }
407: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
408: {
409: Mat_Shell *shell = (Mat_Shell*)mat->data;
410: PetscInt *lrows, *lcols;
411: PetscInt nr, nc;
412: PetscBool congruent;
416: if (x && b) {
417: Vec xt, bt;
418: PetscScalar *vals;
419: PetscInt *grows,*gcols,i,st,nl;
421: PetscMalloc2(n,&grows,n,&gcols);
422: for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
423: for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
424: PetscCalloc1(n,&vals);
426: MatCreateVecs(mat,&xt,&bt);
427: VecCopy(x,xt);
428: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
429: VecAssemblyBegin(xt);
430: VecAssemblyEnd(xt);
431: VecAXPY(xt,-1.0,x); /* xt = [0, -x2] */
432: MatMult(mat,xt,bt); /* bt = [-A12*x2,-A22*x2] */
433: VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* bt = [-A12*x2,0] */
434: VecAssemblyBegin(bt);
435: VecAssemblyEnd(bt);
436: VecAXPY(b,1.0,bt); /* b = [b1 - A12*x2, b2] */
437: VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* b = [b1 - A12*x2, 0] */
438: VecAssemblyBegin(bt);
439: VecAssemblyEnd(bt);
440: PetscFree(vals);
442: VecGetOwnershipRange(xt,&st,NULL);
443: VecGetLocalSize(xt,&nl);
444: VecGetArray(xt,&vals);
445: for (i = 0; i < nl; i++) {
446: PetscInt g = i + st;
447: if (g > mat->rmap->N) continue;
448: if (PetscAbsScalar(vals[i]) == 0.0) continue;
449: VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);
450: }
451: VecRestoreArray(xt,&vals);
452: VecAssemblyBegin(b);
453: VecAssemblyEnd(b); /* b = [b1 - A12*x2, x2 * diag] */
454: VecDestroy(&xt);
455: VecDestroy(&bt);
456: PetscFree2(grows,gcols);
457: }
458: PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);
459: MatHasCongruentLayouts(mat,&congruent);
460: if (congruent) {
461: nc = nr;
462: lcols = lrows;
463: } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
464: PetscInt i,nt,*t;
466: PetscMalloc1(n,&t);
467: for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
468: PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);
469: PetscFree(t);
470: }
471: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);
472: if (!congruent) {
473: PetscFree(lcols);
474: }
475: PetscFree(lrows);
476: if (shell->axpy) {
477: MatZeroRowsColumns(shell->axpy,n,rowscols,0.0,NULL,NULL);
478: }
479: return(0);
480: }
482: PetscErrorCode MatDestroy_Shell(Mat mat)
483: {
484: PetscErrorCode ierr;
485: Mat_Shell *shell = (Mat_Shell*)mat->data;
486: MatShellMatFunctionList matmat;
489: if (shell->ops->destroy) {
490: (*shell->ops->destroy)(mat);
491: }
492: PetscMemzero(shell->ops,sizeof(struct _MatShellOps));
493: VecDestroy(&shell->left);
494: VecDestroy(&shell->right);
495: VecDestroy(&shell->dshift);
496: VecDestroy(&shell->left_work);
497: VecDestroy(&shell->right_work);
498: VecDestroy(&shell->left_add_work);
499: VecDestroy(&shell->right_add_work);
500: VecDestroy(&shell->axpy_left);
501: VecDestroy(&shell->axpy_right);
502: MatDestroy(&shell->axpy);
503: VecDestroy(&shell->zvals_w);
504: VecDestroy(&shell->zvals);
505: VecScatterDestroy(&shell->zvals_sct_c);
506: VecScatterDestroy(&shell->zvals_sct_r);
507: ISDestroy(&shell->zrows);
508: ISDestroy(&shell->zcols);
510: matmat = shell->matmat;
511: while (matmat) {
512: MatShellMatFunctionList next = matmat->next;
514: PetscObjectComposeFunction((PetscObject)mat,matmat->composedname,NULL);
515: PetscFree(matmat->composedname);
516: PetscFree(matmat->resultname);
517: PetscFree(matmat);
518: matmat = next;
519: }
520: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
521: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
522: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);
523: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
524: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
525: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
526: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetMatProductOperation_C",NULL);
527: PetscFree(mat->data);
528: return(0);
529: }
531: typedef struct {
532: PetscErrorCode (*numeric)(Mat,Mat,Mat,void*);
533: PetscErrorCode (*destroy)(void*);
534: void *userdata;
535: Mat B;
536: Mat Bt;
537: Mat axpy;
538: } MatMatDataShell;
540: static PetscErrorCode DestroyMatMatDataShell(void *data)
541: {
542: MatMatDataShell *mmdata = (MatMatDataShell *)data;
546: if (mmdata->destroy) {
547: (*mmdata->destroy)(mmdata->userdata);
548: }
549: MatDestroy(&mmdata->B);
550: MatDestroy(&mmdata->Bt);
551: MatDestroy(&mmdata->axpy);
552: PetscFree(mmdata);
553: return(0);
554: }
556: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
557: {
558: PetscErrorCode ierr;
559: Mat_Product *product;
560: Mat A, B;
561: MatMatDataShell *mdata;
562: PetscScalar zero = 0.0;
565: MatCheckProduct(D,1);
566: product = D->product;
567: if (!product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty");
568: A = product->A;
569: B = product->B;
570: mdata = (MatMatDataShell*)product->data;
571: if (mdata->numeric) {
572: Mat_Shell *shell = (Mat_Shell*)A->data;
573: PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
574: PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
575: PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;
577: if (shell->managescalingshifts) {
578: if (shell->zcols || shell->zrows) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProduct not supported with zeroed rows/columns");
579: if (shell->right || shell->left) {
580: useBmdata = PETSC_TRUE;
581: if (!mdata->B) {
582: MatDuplicate(B,MAT_SHARE_NONZERO_PATTERN,&mdata->B);
583: } else {
584: newB = PETSC_FALSE;
585: }
586: MatCopy(B,mdata->B,SAME_NONZERO_PATTERN);
587: }
588: switch (product->type) {
589: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
590: if (shell->right) {
591: MatDiagonalScale(mdata->B,shell->right,NULL);
592: }
593: break;
594: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
595: if (shell->left) {
596: MatDiagonalScale(mdata->B,shell->left,NULL);
597: }
598: break;
599: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
600: if (shell->right) {
601: MatDiagonalScale(mdata->B,NULL,shell->right);
602: }
603: break;
604: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
605: if (shell->right && shell->left) {
606: PetscBool flg;
608: VecEqual(shell->right,shell->left,&flg);
609: if (!flg) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
610: }
611: if (shell->right) {
612: MatDiagonalScale(mdata->B,NULL,shell->right);
613: }
614: break;
615: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
616: if (shell->right && shell->left) {
617: PetscBool flg;
619: VecEqual(shell->right,shell->left,&flg);
620: if (!flg) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices because left scaling != from right scaling",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
621: }
622: if (shell->right) {
623: MatDiagonalScale(mdata->B,shell->right,NULL);
624: }
625: break;
626: default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
627: }
628: }
629: /* allow the user to call MatMat operations on D */
630: D->product = NULL;
631: D->ops->productsymbolic = NULL;
632: D->ops->productnumeric = NULL;
634: (*mdata->numeric)(A,useBmdata ? mdata->B : B,D,mdata->userdata);
636: /* clear any leftover user data and restore D pointers */
637: MatProductClear(D);
638: D->ops->productsymbolic = stashsym;
639: D->ops->productnumeric = stashnum;
640: D->product = product;
642: if (shell->managescalingshifts) {
643: MatScale(D,shell->vscale);
644: switch (product->type) {
645: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
646: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
647: if (shell->left) {
648: MatDiagonalScale(D,shell->left,NULL);
649: if (shell->dshift || shell->vshift != zero) {
650: if (!shell->left_work) {MatCreateVecs(A,NULL,&shell->left_work);}
651: if (shell->dshift) {
652: VecCopy(shell->dshift,shell->left_work);
653: VecShift(shell->left_work,shell->vshift);
654: VecPointwiseMult(shell->left_work,shell->left_work,shell->left);
655: } else {
656: VecSet(shell->left_work,shell->vshift);
657: }
658: if (product->type == MATPRODUCT_ABt) {
659: MatReuse reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
660: MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;
662: MatTranspose(mdata->B,reuse,&mdata->Bt);
663: MatDiagonalScale(mdata->Bt,shell->left_work,NULL);
664: MatAXPY(D,1.0,mdata->Bt,str);
665: } else {
666: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
668: MatDiagonalScale(mdata->B,shell->left_work,NULL);
669: MatAXPY(D,1.0,mdata->B,str);
670: }
671: }
672: }
673: break;
674: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
675: if (shell->right) {
676: MatDiagonalScale(D,shell->right,NULL);
677: if (shell->dshift || shell->vshift != zero) {
678: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
680: if (!shell->right_work) {MatCreateVecs(A,&shell->right_work,NULL);}
681: if (shell->dshift) {
682: VecCopy(shell->dshift,shell->right_work);
683: VecShift(shell->right_work,shell->vshift);
684: VecPointwiseMult(shell->right_work,shell->right_work,shell->right);
685: } else {
686: VecSet(shell->right_work,shell->vshift);
687: }
688: MatDiagonalScale(mdata->B,shell->right_work,NULL);
689: MatAXPY(D,1.0,mdata->B,str);
690: }
691: }
692: break;
693: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
694: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
695: if (shell->dshift || shell->vshift != zero) SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices with diagonal shift",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
696: break;
697: default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
698: }
699: if (shell->axpy && shell->axpy_vscale != zero) {
700: Mat X;
701: PetscObjectState axpy_state;
702: MatStructure str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */
704: MatShellGetContext(shell->axpy,(void *)&X);
705: PetscObjectStateGet((PetscObject)X,&axpy_state);
706: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
707: if (!mdata->axpy) {
708: str = DIFFERENT_NONZERO_PATTERN;
709: MatProductCreate(shell->axpy,B,NULL,&mdata->axpy);
710: MatProductSetType(mdata->axpy,product->type);
711: MatProductSetFromOptions(mdata->axpy);
712: MatProductSymbolic(mdata->axpy);
713: } else { /* May be that shell->axpy has changed */
714: PetscBool flg;
716: MatProductReplaceMats(shell->axpy,B,NULL,mdata->axpy);
717: MatHasOperation(mdata->axpy,MATOP_PRODUCTSYMBOLIC,&flg);
718: if (!flg) {
719: str = DIFFERENT_NONZERO_PATTERN;
720: MatProductSetFromOptions(mdata->axpy);
721: MatProductSymbolic(mdata->axpy);
722: }
723: }
724: MatProductNumeric(mdata->axpy);
725: MatAXPY(D,shell->axpy_vscale,mdata->axpy,str);
726: }
727: }
728: } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation");
729: return(0);
730: }
732: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
733: {
734: PetscErrorCode ierr;
735: Mat_Product *product;
736: Mat A,B;
737: MatShellMatFunctionList matmat;
738: Mat_Shell *shell;
739: PetscBool flg;
740: char composedname[256];
741: MatMatDataShell *mdata;
744: MatCheckProduct(D,1);
745: product = D->product;
746: if (product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty");
747: A = product->A;
748: B = product->B;
749: shell = (Mat_Shell*)A->data;
750: matmat = shell->matmat;
751: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
752: while (matmat) {
753: PetscStrcmp(composedname,matmat->composedname,&flg);
754: flg = (PetscBool)(flg && (matmat->ptype == product->type));
755: if (flg) break;
756: matmat = matmat->next;
757: }
758: if (!flg) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Composedname \"%s\" for product type %s not found",composedname,MatProductTypes[product->type]);
759: switch (product->type) {
760: case MATPRODUCT_AB:
761: MatSetSizes(D,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
762: break;
763: case MATPRODUCT_AtB:
764: MatSetSizes(D,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
765: break;
766: case MATPRODUCT_ABt:
767: MatSetSizes(D,A->rmap->n,B->rmap->n,A->rmap->N,B->rmap->N);
768: break;
769: case MATPRODUCT_RARt:
770: MatSetSizes(D,B->rmap->n,B->rmap->n,B->rmap->N,B->rmap->N);
771: break;
772: case MATPRODUCT_PtAP:
773: MatSetSizes(D,B->cmap->n,B->cmap->n,B->cmap->N,B->cmap->N);
774: break;
775: default: SETERRQ3(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProductSymbolic type %s not supported for %s and %s matrices",MatProductTypes[product->type],((PetscObject)A)->type_name,((PetscObject)B)->type_name);
776: }
777: /* respect users who passed in a matrix for which resultname is the base type */
778: if (matmat->resultname) {
779: PetscObjectBaseTypeCompare((PetscObject)D,matmat->resultname,&flg);
780: if (!flg) {
781: MatSetType(D,matmat->resultname);
782: }
783: }
784: /* If matrix type was not set or different, we need to reset this pointers */
785: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
786: D->ops->productnumeric = MatProductNumeric_Shell_X;
787: /* attach product data */
788: PetscNew(&mdata);
789: mdata->numeric = matmat->numeric;
790: mdata->destroy = matmat->destroy;
791: if (matmat->symbolic) {
792: (*matmat->symbolic)(A,B,D,&mdata->userdata);
793: } else { /* call general setup if symbolic operation not provided */
794: MatSetUp(D);
795: }
796: if (!D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product disappeared after user symbolic phase");
797: if (D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product data not empty after user symbolic phase");
798: D->product->data = mdata;
799: D->product->destroy = DestroyMatMatDataShell;
800: /* Be sure to reset these pointers if the user did something unexpected */
801: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
802: D->ops->productnumeric = MatProductNumeric_Shell_X;
803: return(0);
804: }
806: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
807: {
808: PetscErrorCode ierr;
809: Mat_Product *product;
810: Mat A,B;
811: MatShellMatFunctionList matmat;
812: Mat_Shell *shell;
813: PetscBool flg;
814: char composedname[256];
817: MatCheckProduct(D,1);
818: product = D->product;
819: A = product->A;
820: B = product->B;
821: MatIsShell(A,&flg);
822: if (!flg) return(0);
823: shell = (Mat_Shell*)A->data;
824: matmat = shell->matmat;
825: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
826: while (matmat) {
827: PetscStrcmp(composedname,matmat->composedname,&flg);
828: flg = (PetscBool)(flg && (matmat->ptype == product->type));
829: if (flg) break;
830: matmat = matmat->next;
831: }
832: if (flg) { D->ops->productsymbolic = MatProductSymbolic_Shell_X; }
833: else { PetscInfo2(D," symbolic product %s not registered for product type %s\n",composedname,MatProductTypes[product->type]); }
834: return(0);
835: }
837: static PetscErrorCode MatShellSetMatProductOperation_Private(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void*),char *composedname,const char *resultname)
838: {
839: PetscBool flg;
840: PetscErrorCode ierr;
841: Mat_Shell *shell;
842: MatShellMatFunctionList matmat;
845: if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing numeric routine");
846: if (!composedname) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing composed name");
848: /* add product callback */
849: shell = (Mat_Shell*)A->data;
850: matmat = shell->matmat;
851: if (!matmat) {
852: PetscNew(&shell->matmat);
853: matmat = shell->matmat;
854: } else {
855: MatShellMatFunctionList entry = matmat;
856: while (entry) {
857: PetscStrcmp(composedname,entry->composedname,&flg);
858: flg = (PetscBool)(flg && (entry->ptype == ptype));
859: if (flg) break;
860: matmat = entry;
861: entry = entry->next;
862: }
863: if (!flg) {
864: PetscNew(&matmat->next);
865: matmat = matmat->next;
866: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"This should not happen");
867: }
869: matmat->symbolic = symbolic;
870: matmat->numeric = numeric;
871: matmat->destroy = destroy;
872: matmat->ptype = ptype;
873: PetscFree(matmat->composedname);
874: PetscFree(matmat->resultname);
875: PetscStrallocpy(composedname,&matmat->composedname);
876: PetscStrallocpy(resultname,&matmat->resultname);
877: PetscInfo3(A,"Composing %s for product type %s with result %s\n",matmat->composedname,MatProductTypes[matmat->ptype],matmat->resultname ? matmat->resultname : "not specified");
878: PetscObjectComposeFunction((PetscObject)A,matmat->composedname,MatProductSetFromOptions_Shell_X);
879: return(0);
880: }
882: /*@C
883: MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a shell matrix.
885: Logically Collective on Mat
887: Input Parameters:
888: + A - the shell matrix
889: . ptype - the product type
890: . symbolic - the function for the symbolic phase (can be NULL)
891: . numeric - the function for the numerical phase
892: . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL)
893: . Btype - the matrix type for the matrix to be multiplied against
894: - Ctype - the matrix type for the result (can be NULL)
896: Level: advanced
898: Usage:
899: $ extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**);
900: $ extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*);
901: $ extern PetscErrorCode userdestroy(void*);
902: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
903: $ MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE);
904: $ [ create B of type SEQAIJ etc..]
905: $ MatProductCreate(A,B,NULL,&C);
906: $ MatProductSetType(C,MATPRODUCT_AB);
907: $ MatProductSetFromOptions(C);
908: $ MatProductSymbolic(C); -> actually runs the user defined symbolic operation
909: $ MatProductNumeric(C); -> actually runs the user defined numeric operation
910: $ [ use C = A*B ]
912: Notes:
913: MATPRODUCT_ABC is not supported yet. Not supported in Fortran.
914: If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL.
915: Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
916: PETSc will take care of calling the user-defined callbacks.
917: It is allowed to specify the same callbacks for different Btype matrix types.
918: The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence.
920: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp()
921: @*/
922: PetscErrorCode MatShellSetMatProductOperation(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void *),MatType Btype,MatType Ctype)
923: {
929: if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]);
930: if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4");
933: PetscTryMethod(A,"MatShellSetMatProductOperation_C",(Mat,MatProductType,PetscErrorCode(*)(Mat,Mat,Mat,void**),PetscErrorCode(*)(Mat,Mat,Mat,void*),PetscErrorCode(*)(void*),MatType,MatType),(A,ptype,symbolic,numeric,destroy,Btype,Ctype));
934: return(0);
935: }
937: PetscErrorCode MatShellSetMatProductOperation_Shell(Mat A,MatProductType ptype,PetscErrorCode (*symbolic)(Mat,Mat,Mat,void**),PetscErrorCode (*numeric)(Mat,Mat,Mat,void*),PetscErrorCode (*destroy)(void *),MatType Btype,MatType Ctype)
938: {
939: PetscBool flg;
941: char composedname[256];
942: MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
943: PetscMPIInt size;
947: while (Bnames) { /* user passed in the root name */
948: PetscStrcmp(Btype,Bnames->rname,&flg);
949: if (flg) break;
950: Bnames = Bnames->next;
951: }
952: while (Cnames) { /* user passed in the root name */
953: PetscStrcmp(Ctype,Cnames->rname,&flg);
954: if (flg) break;
955: Cnames = Cnames->next;
956: }
957: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
958: Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
959: Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
960: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,Btype);
961: MatShellSetMatProductOperation_Private(A,ptype,symbolic,numeric,destroy,composedname,Ctype);
962: return(0);
963: }
965: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
966: {
967: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
968: PetscErrorCode ierr;
969: PetscBool matflg;
970: MatShellMatFunctionList matmatA;
973: MatIsShell(B,&matflg);
974: if (!matflg) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix %s not derived from MATSHELL",((PetscObject)B)->type_name);
976: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
977: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
979: if (shellA->ops->copy) {
980: (*shellA->ops->copy)(A,B,str);
981: }
982: shellB->vscale = shellA->vscale;
983: shellB->vshift = shellA->vshift;
984: if (shellA->dshift) {
985: if (!shellB->dshift) {
986: VecDuplicate(shellA->dshift,&shellB->dshift);
987: }
988: VecCopy(shellA->dshift,shellB->dshift);
989: } else {
990: VecDestroy(&shellB->dshift);
991: }
992: if (shellA->left) {
993: if (!shellB->left) {
994: VecDuplicate(shellA->left,&shellB->left);
995: }
996: VecCopy(shellA->left,shellB->left);
997: } else {
998: VecDestroy(&shellB->left);
999: }
1000: if (shellA->right) {
1001: if (!shellB->right) {
1002: VecDuplicate(shellA->right,&shellB->right);
1003: }
1004: VecCopy(shellA->right,shellB->right);
1005: } else {
1006: VecDestroy(&shellB->right);
1007: }
1008: MatDestroy(&shellB->axpy);
1009: shellB->axpy_vscale = 0.0;
1010: shellB->axpy_state = 0;
1011: if (shellA->axpy) {
1012: PetscObjectReference((PetscObject)shellA->axpy);
1013: shellB->axpy = shellA->axpy;
1014: shellB->axpy_vscale = shellA->axpy_vscale;
1015: shellB->axpy_state = shellA->axpy_state;
1016: }
1017: if (shellA->zrows) {
1018: ISDuplicate(shellA->zrows,&shellB->zrows);
1019: if (shellA->zcols) {
1020: ISDuplicate(shellA->zcols,&shellB->zcols);
1021: }
1022: VecDuplicate(shellA->zvals,&shellB->zvals);
1023: VecCopy(shellA->zvals,shellB->zvals);
1024: VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
1025: PetscObjectReference((PetscObject)shellA->zvals_sct_r);
1026: PetscObjectReference((PetscObject)shellA->zvals_sct_c);
1027: shellB->zvals_sct_r = shellA->zvals_sct_r;
1028: shellB->zvals_sct_c = shellA->zvals_sct_c;
1029: }
1031: matmatA = shellA->matmat;
1032: if (matmatA) {
1033: while (matmatA->next) {
1034: MatShellSetMatProductOperation_Private(B,matmatA->ptype,matmatA->symbolic,matmatA->numeric,matmatA->destroy,matmatA->composedname,matmatA->resultname);
1035: matmatA = matmatA->next;
1036: }
1037: }
1038: return(0);
1039: }
1041: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
1042: {
1044: void *ctx;
1047: MatShellGetContext(mat,&ctx);
1048: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
1049: PetscObjectChangeTypeName((PetscObject)(*M),((PetscObject)mat)->type_name);
1050: if (op != MAT_DO_NOT_COPY_VALUES) {
1051: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
1052: }
1053: return(0);
1054: }
1056: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
1057: {
1058: Mat_Shell *shell = (Mat_Shell*)A->data;
1059: PetscErrorCode ierr;
1060: Vec xx;
1061: PetscObjectState instate,outstate;
1064: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
1065: MatShellPreZeroRight(A,x,&xx);
1066: MatShellPreScaleRight(A,xx,&xx);
1067: PetscObjectStateGet((PetscObject)y, &instate);
1068: (*shell->ops->mult)(A,xx,y);
1069: PetscObjectStateGet((PetscObject)y, &outstate);
1070: if (instate == outstate) {
1071: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1072: PetscObjectStateIncrease((PetscObject)y);
1073: }
1074: MatShellShiftAndScale(A,xx,y);
1075: MatShellPostScaleLeft(A,y);
1076: MatShellPostZeroLeft(A,y);
1078: if (shell->axpy) {
1079: Mat X;
1080: PetscObjectState axpy_state;
1082: MatShellGetContext(shell->axpy,(void *)&X);
1083: PetscObjectStateGet((PetscObject)X,&axpy_state);
1084: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1086: MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1087: VecCopy(x,shell->axpy_right);
1088: MatMult(shell->axpy,shell->axpy_right,shell->axpy_left);
1089: VecAXPY(y,shell->axpy_vscale,shell->axpy_left);
1090: }
1091: return(0);
1092: }
1094: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
1095: {
1096: Mat_Shell *shell = (Mat_Shell*)A->data;
1100: if (y == z) {
1101: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
1102: MatMult(A,x,shell->right_add_work);
1103: VecAXPY(z,1.0,shell->right_add_work);
1104: } else {
1105: MatMult(A,x,z);
1106: VecAXPY(z,1.0,y);
1107: }
1108: return(0);
1109: }
1111: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
1112: {
1113: Mat_Shell *shell = (Mat_Shell*)A->data;
1114: PetscErrorCode ierr;
1115: Vec xx;
1116: PetscObjectState instate,outstate;
1119: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
1120: MatShellPreZeroLeft(A,x,&xx);
1121: MatShellPreScaleLeft(A,xx,&xx);
1122: PetscObjectStateGet((PetscObject)y, &instate);
1123: (*shell->ops->multtranspose)(A,xx,y);
1124: PetscObjectStateGet((PetscObject)y, &outstate);
1125: if (instate == outstate) {
1126: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1127: PetscObjectStateIncrease((PetscObject)y);
1128: }
1129: MatShellShiftAndScale(A,xx,y);
1130: MatShellPostScaleRight(A,y);
1131: MatShellPostZeroRight(A,y);
1133: if (shell->axpy) {
1134: Mat X;
1135: PetscObjectState axpy_state;
1137: MatShellGetContext(shell->axpy,(void *)&X);
1138: PetscObjectStateGet((PetscObject)X,&axpy_state);
1139: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1140: MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1141: VecCopy(x,shell->axpy_left);
1142: MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);
1143: VecAXPY(y,shell->axpy_vscale,shell->axpy_right);
1144: }
1145: return(0);
1146: }
1148: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
1149: {
1150: Mat_Shell *shell = (Mat_Shell*)A->data;
1154: if (y == z) {
1155: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
1156: MatMultTranspose(A,x,shell->left_add_work);
1157: VecAXPY(z,1.0,shell->left_add_work);
1158: } else {
1159: MatMultTranspose(A,x,z);
1160: VecAXPY(z,1.0,y);
1161: }
1162: return(0);
1163: }
1165: /*
1166: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1167: */
1168: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
1169: {
1170: Mat_Shell *shell = (Mat_Shell*)A->data;
1174: if (shell->ops->getdiagonal) {
1175: (*shell->ops->getdiagonal)(A,v);
1176: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
1177: VecScale(v,shell->vscale);
1178: if (shell->dshift) {
1179: VecAXPY(v,1.0,shell->dshift);
1180: }
1181: VecShift(v,shell->vshift);
1182: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
1183: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
1184: if (shell->zrows) {
1185: VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
1186: VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
1187: }
1188: if (shell->axpy) {
1189: Mat X;
1190: PetscObjectState axpy_state;
1192: MatShellGetContext(shell->axpy,(void *)&X);
1193: PetscObjectStateGet((PetscObject)X,&axpy_state);
1194: if (shell->axpy_state != axpy_state) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ORDER,"Invalid AXPY state: cannot modify the X matrix passed to MatAXPY(Y,a,X,...)");
1195: MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);
1196: MatGetDiagonal(shell->axpy,shell->axpy_left);
1197: VecAXPY(v,shell->axpy_vscale,shell->axpy_left);
1198: }
1199: return(0);
1200: }
1202: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
1203: {
1204: Mat_Shell *shell = (Mat_Shell*)Y->data;
1206: PetscBool flg;
1209: MatHasCongruentLayouts(Y,&flg);
1210: if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
1211: if (shell->left || shell->right) {
1212: if (!shell->dshift) {
1213: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
1214: VecSet(shell->dshift,a);
1215: } else {
1216: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
1217: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
1218: VecShift(shell->dshift,a);
1219: }
1220: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
1221: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
1222: } else shell->vshift += a;
1223: if (shell->zrows) {
1224: VecShift(shell->zvals,a);
1225: }
1226: return(0);
1227: }
1229: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
1230: {
1231: Mat_Shell *shell = (Mat_Shell*)A->data;
1235: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
1236: if (shell->left || shell->right) {
1237: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
1238: if (shell->left && shell->right) {
1239: VecPointwiseDivide(shell->right_work,D,shell->left);
1240: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
1241: } else if (shell->left) {
1242: VecPointwiseDivide(shell->right_work,D,shell->left);
1243: } else {
1244: VecPointwiseDivide(shell->right_work,D,shell->right);
1245: }
1246: VecAXPY(shell->dshift,s,shell->right_work);
1247: } else {
1248: VecAXPY(shell->dshift,s,D);
1249: }
1250: return(0);
1251: }
1253: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
1254: {
1255: Mat_Shell *shell = (Mat_Shell*)A->data;
1256: Vec d;
1258: PetscBool flg;
1261: MatHasCongruentLayouts(A,&flg);
1262: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
1263: if (ins == INSERT_VALUES) {
1264: if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
1265: VecDuplicate(D,&d);
1266: MatGetDiagonal(A,d);
1267: MatDiagonalSet_Shell_Private(A,d,-1.);
1268: MatDiagonalSet_Shell_Private(A,D,1.);
1269: VecDestroy(&d);
1270: if (shell->zrows) {
1271: VecCopy(D,shell->zvals);
1272: }
1273: } else {
1274: MatDiagonalSet_Shell_Private(A,D,1.);
1275: if (shell->zrows) {
1276: VecAXPY(shell->zvals,1.0,D);
1277: }
1278: }
1279: return(0);
1280: }
1282: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
1283: {
1284: Mat_Shell *shell = (Mat_Shell*)Y->data;
1288: shell->vscale *= a;
1289: shell->vshift *= a;
1290: if (shell->dshift) {
1291: VecScale(shell->dshift,a);
1292: }
1293: shell->axpy_vscale *= a;
1294: if (shell->zrows) {
1295: VecScale(shell->zvals,a);
1296: }
1297: return(0);
1298: }
1300: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
1301: {
1302: Mat_Shell *shell = (Mat_Shell*)Y->data;
1306: if (left) {
1307: if (!shell->left) {
1308: VecDuplicate(left,&shell->left);
1309: VecCopy(left,shell->left);
1310: } else {
1311: VecPointwiseMult(shell->left,shell->left,left);
1312: }
1313: if (shell->zrows) {
1314: VecPointwiseMult(shell->zvals,shell->zvals,left);
1315: }
1316: }
1317: if (right) {
1318: if (!shell->right) {
1319: VecDuplicate(right,&shell->right);
1320: VecCopy(right,shell->right);
1321: } else {
1322: VecPointwiseMult(shell->right,shell->right,right);
1323: }
1324: if (shell->zrows) {
1325: if (!shell->left_work) {
1326: MatCreateVecs(Y,NULL,&shell->left_work);
1327: }
1328: VecSet(shell->zvals_w,1.0);
1329: VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
1330: VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
1331: VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
1332: }
1333: }
1334: if (shell->axpy) {
1335: MatDiagonalScale(shell->axpy,left,right);
1336: }
1337: return(0);
1338: }
1340: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
1341: {
1342: Mat_Shell *shell = (Mat_Shell*)Y->data;
1346: if (t == MAT_FINAL_ASSEMBLY) {
1347: shell->vshift = 0.0;
1348: shell->vscale = 1.0;
1349: shell->axpy_vscale = 0.0;
1350: shell->axpy_state = 0;
1351: VecDestroy(&shell->dshift);
1352: VecDestroy(&shell->left);
1353: VecDestroy(&shell->right);
1354: MatDestroy(&shell->axpy);
1355: VecDestroy(&shell->axpy_left);
1356: VecDestroy(&shell->axpy_right);
1357: VecScatterDestroy(&shell->zvals_sct_c);
1358: VecScatterDestroy(&shell->zvals_sct_r);
1359: ISDestroy(&shell->zrows);
1360: ISDestroy(&shell->zcols);
1361: }
1362: return(0);
1363: }
1365: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
1366: {
1368: *missing = PETSC_FALSE;
1369: return(0);
1370: }
1372: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
1373: {
1374: Mat_Shell *shell = (Mat_Shell*)Y->data;
1378: if (X == Y) {
1379: MatScale(Y,1.0 + a);
1380: return(0);
1381: }
1382: if (!shell->axpy) {
1383: MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);
1384: shell->axpy_vscale = a;
1385: PetscObjectStateGet((PetscObject)X,&shell->axpy_state);
1386: } else {
1387: MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);
1388: }
1389: return(0);
1390: }
1392: static struct _MatOps MatOps_Values = {NULL,
1393: NULL,
1394: NULL,
1395: NULL,
1396: /* 4*/ MatMultAdd_Shell,
1397: NULL,
1398: MatMultTransposeAdd_Shell,
1399: NULL,
1400: NULL,
1401: NULL,
1402: /*10*/ NULL,
1403: NULL,
1404: NULL,
1405: NULL,
1406: NULL,
1407: /*15*/ NULL,
1408: NULL,
1409: NULL,
1410: MatDiagonalScale_Shell,
1411: NULL,
1412: /*20*/ NULL,
1413: MatAssemblyEnd_Shell,
1414: NULL,
1415: NULL,
1416: /*24*/ MatZeroRows_Shell,
1417: NULL,
1418: NULL,
1419: NULL,
1420: NULL,
1421: /*29*/ NULL,
1422: NULL,
1423: NULL,
1424: NULL,
1425: NULL,
1426: /*34*/ MatDuplicate_Shell,
1427: NULL,
1428: NULL,
1429: NULL,
1430: NULL,
1431: /*39*/ MatAXPY_Shell,
1432: NULL,
1433: NULL,
1434: NULL,
1435: MatCopy_Shell,
1436: /*44*/ NULL,
1437: MatScale_Shell,
1438: MatShift_Shell,
1439: MatDiagonalSet_Shell,
1440: MatZeroRowsColumns_Shell,
1441: /*49*/ NULL,
1442: NULL,
1443: NULL,
1444: NULL,
1445: NULL,
1446: /*54*/ NULL,
1447: NULL,
1448: NULL,
1449: NULL,
1450: NULL,
1451: /*59*/ NULL,
1452: MatDestroy_Shell,
1453: NULL,
1454: MatConvertFrom_Shell,
1455: NULL,
1456: /*64*/ NULL,
1457: NULL,
1458: NULL,
1459: NULL,
1460: NULL,
1461: /*69*/ NULL,
1462: NULL,
1463: MatConvert_Shell,
1464: NULL,
1465: NULL,
1466: /*74*/ NULL,
1467: NULL,
1468: NULL,
1469: NULL,
1470: NULL,
1471: /*79*/ NULL,
1472: NULL,
1473: NULL,
1474: NULL,
1475: NULL,
1476: /*84*/ NULL,
1477: NULL,
1478: NULL,
1479: NULL,
1480: NULL,
1481: /*89*/ NULL,
1482: NULL,
1483: NULL,
1484: NULL,
1485: NULL,
1486: /*94*/ NULL,
1487: NULL,
1488: NULL,
1489: NULL,
1490: NULL,
1491: /*99*/ NULL,
1492: NULL,
1493: NULL,
1494: NULL,
1495: NULL,
1496: /*104*/ NULL,
1497: NULL,
1498: NULL,
1499: NULL,
1500: NULL,
1501: /*109*/ NULL,
1502: NULL,
1503: NULL,
1504: NULL,
1505: MatMissingDiagonal_Shell,
1506: /*114*/ NULL,
1507: NULL,
1508: NULL,
1509: NULL,
1510: NULL,
1511: /*119*/ NULL,
1512: NULL,
1513: NULL,
1514: NULL,
1515: NULL,
1516: /*124*/ NULL,
1517: NULL,
1518: NULL,
1519: NULL,
1520: NULL,
1521: /*129*/ NULL,
1522: NULL,
1523: NULL,
1524: NULL,
1525: NULL,
1526: /*134*/ NULL,
1527: NULL,
1528: NULL,
1529: NULL,
1530: NULL,
1531: /*139*/ NULL,
1532: NULL,
1533: NULL
1534: };
1536: PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx)
1537: {
1538: Mat_Shell *shell = (Mat_Shell*)mat->data;
1541: shell->ctx = ctx;
1542: return(0);
1543: }
1545: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1546: {
1550: PetscFree(mat->defaultvectype);
1551: PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
1552: return(0);
1553: }
1555: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1556: {
1557: Mat_Shell *shell = (Mat_Shell*)A->data;
1560: shell->managescalingshifts = PETSC_FALSE;
1561: A->ops->diagonalset = NULL;
1562: A->ops->diagonalscale = NULL;
1563: A->ops->scale = NULL;
1564: A->ops->shift = NULL;
1565: A->ops->axpy = NULL;
1566: return(0);
1567: }
1569: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1570: {
1571: Mat_Shell *shell = (Mat_Shell*)mat->data;
1574: switch (op) {
1575: case MATOP_DESTROY:
1576: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1577: break;
1578: case MATOP_VIEW:
1579: if (!mat->ops->viewnative) {
1580: mat->ops->viewnative = mat->ops->view;
1581: }
1582: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1583: break;
1584: case MATOP_COPY:
1585: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1586: break;
1587: case MATOP_DIAGONAL_SET:
1588: case MATOP_DIAGONAL_SCALE:
1589: case MATOP_SHIFT:
1590: case MATOP_SCALE:
1591: case MATOP_AXPY:
1592: case MATOP_ZERO_ROWS:
1593: case MATOP_ZERO_ROWS_COLUMNS:
1594: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1595: (((void(**)(void))mat->ops)[op]) = f;
1596: break;
1597: case MATOP_GET_DIAGONAL:
1598: if (shell->managescalingshifts) {
1599: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1600: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1601: } else {
1602: shell->ops->getdiagonal = NULL;
1603: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1604: }
1605: break;
1606: case MATOP_MULT:
1607: if (shell->managescalingshifts) {
1608: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1609: mat->ops->mult = MatMult_Shell;
1610: } else {
1611: shell->ops->mult = NULL;
1612: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1613: }
1614: break;
1615: case MATOP_MULT_TRANSPOSE:
1616: if (shell->managescalingshifts) {
1617: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1618: mat->ops->multtranspose = MatMultTranspose_Shell;
1619: } else {
1620: shell->ops->multtranspose = NULL;
1621: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1622: }
1623: break;
1624: default:
1625: (((void(**)(void))mat->ops)[op]) = f;
1626: break;
1627: }
1628: return(0);
1629: }
1631: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1632: {
1633: Mat_Shell *shell = (Mat_Shell*)mat->data;
1636: switch (op) {
1637: case MATOP_DESTROY:
1638: *f = (void (*)(void))shell->ops->destroy;
1639: break;
1640: case MATOP_VIEW:
1641: *f = (void (*)(void))mat->ops->view;
1642: break;
1643: case MATOP_COPY:
1644: *f = (void (*)(void))shell->ops->copy;
1645: break;
1646: case MATOP_DIAGONAL_SET:
1647: case MATOP_DIAGONAL_SCALE:
1648: case MATOP_SHIFT:
1649: case MATOP_SCALE:
1650: case MATOP_AXPY:
1651: case MATOP_ZERO_ROWS:
1652: case MATOP_ZERO_ROWS_COLUMNS:
1653: *f = (((void (**)(void))mat->ops)[op]);
1654: break;
1655: case MATOP_GET_DIAGONAL:
1656: if (shell->ops->getdiagonal)
1657: *f = (void (*)(void))shell->ops->getdiagonal;
1658: else
1659: *f = (((void (**)(void))mat->ops)[op]);
1660: break;
1661: case MATOP_MULT:
1662: if (shell->ops->mult)
1663: *f = (void (*)(void))shell->ops->mult;
1664: else
1665: *f = (((void (**)(void))mat->ops)[op]);
1666: break;
1667: case MATOP_MULT_TRANSPOSE:
1668: if (shell->ops->multtranspose)
1669: *f = (void (*)(void))shell->ops->multtranspose;
1670: else
1671: *f = (((void (**)(void))mat->ops)[op]);
1672: break;
1673: default:
1674: *f = (((void (**)(void))mat->ops)[op]);
1675: }
1676: return(0);
1677: }
1679: /*MC
1680: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
1682: Level: advanced
1684: .seealso: MatCreateShell()
1685: M*/
1687: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1688: {
1689: Mat_Shell *b;
1693: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1695: PetscNewLog(A,&b);
1696: A->data = (void*)b;
1698: b->ctx = NULL;
1699: b->vshift = 0.0;
1700: b->vscale = 1.0;
1701: b->managescalingshifts = PETSC_TRUE;
1702: A->assembled = PETSC_TRUE;
1703: A->preallocated = PETSC_FALSE;
1705: PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1706: PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1707: PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
1708: PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1709: PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1710: PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1711: PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell);
1712: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1713: return(0);
1714: }
1716: /*@C
1717: MatCreateShell - Creates a new matrix class for use with a user-defined
1718: private data storage format.
1720: Collective
1722: Input Parameters:
1723: + comm - MPI communicator
1724: . m - number of local rows (must be given)
1725: . n - number of local columns (must be given)
1726: . M - number of global rows (may be PETSC_DETERMINE)
1727: . N - number of global columns (may be PETSC_DETERMINE)
1728: - ctx - pointer to data needed by the shell matrix routines
1730: Output Parameter:
1731: . A - the matrix
1733: Level: advanced
1735: Usage:
1736: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1737: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
1738: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1739: $ [ Use matrix for operations that have been set ]
1740: $ MatDestroy(mat);
1742: Notes:
1743: The shell matrix type is intended to provide a simple class to use
1744: with KSP (such as, for use with matrix-free methods). You should not
1745: use the shell type if you plan to define a complete matrix class.
1747: Fortran Notes:
1748: To use this from Fortran with a ctx you must write an interface definition for this
1749: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
1750: in as the ctx argument.
1752: PETSc requires that matrices and vectors being used for certain
1753: operations are partitioned accordingly. For example, when
1754: creating a shell matrix, A, that supports parallel matrix-vector
1755: products using MatMult(A,x,y) the user should set the number
1756: of local matrix rows to be the number of local elements of the
1757: corresponding result vector, y. Note that this is information is
1758: required for use of the matrix interface routines, even though
1759: the shell matrix may not actually be physically partitioned.
1760: For example,
1762: $
1763: $ Vec x, y
1764: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1765: $ Mat A
1766: $
1767: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1768: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1769: $ VecGetLocalSize(y,&m);
1770: $ VecGetLocalSize(x,&n);
1771: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1772: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1773: $ MatMult(A,x,y);
1774: $ MatDestroy(&A);
1775: $ VecDestroy(&y);
1776: $ VecDestroy(&x);
1777: $
1780: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these
1781: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
1784: For rectangular matrices do all the scalings and shifts make sense?
1786: Developers Notes:
1787: Regarding shifting and scaling. The general form is
1789: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1791: The order you apply the operations is important. For example if you have a dshift then
1792: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1793: you get s*vscale*A + diag(shift)
1795: A is the user provided function.
1797: KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
1798: for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
1799: an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
1800: each time the MATSHELL matrix has changed.
1802: Matrix product operations (i.e. MatMat, MatTranposeMat etc) can be specified using MatShellSetMatProductOperation()
1804: Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
1805: with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().
1807: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
1808: @*/
1809: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1810: {
1814: MatCreate(comm,A);
1815: MatSetSizes(*A,m,n,M,N);
1816: MatSetType(*A,MATSHELL);
1817: MatShellSetContext(*A,ctx);
1818: MatSetUp(*A);
1819: return(0);
1820: }
1822: /*@
1823: MatShellSetContext - sets the context for a shell matrix
1825: Logically Collective on Mat
1827: Input Parameters:
1828: + mat - the shell matrix
1829: - ctx - the context
1831: Level: advanced
1833: Fortran Notes:
1834: To use this from Fortran you must write a Fortran interface definition for this
1835: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1837: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1838: @*/
1839: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
1840: {
1845: PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1846: return(0);
1847: }
1849: /*@C
1850: MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs()
1852: Logically collective
1854: Input Parameters:
1855: + mat - the shell matrix
1856: - vtype - type to use for creating vectors
1858: Notes:
1860: Level: advanced
1862: .seealso: MatCreateVecs()
1863: @*/
1864: PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype)
1865: {
1869: PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1870: return(0);
1871: }
1873: /*@
1874: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
1875: after MatCreateShell()
1877: Logically Collective on Mat
1879: Input Parameter:
1880: . mat - the shell matrix
1882: Level: advanced
1884: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1885: @*/
1886: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1887: {
1892: PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1893: return(0);
1894: }
1896: /*@C
1897: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
1899: Logically Collective on Mat
1901: Input Parameters:
1902: + mat - the shell matrix
1903: . f - the function
1904: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1905: - ctx - an optional context for the function
1907: Output Parameter:
1908: . flg - PETSC_TRUE if the multiply is likely correct
1910: Options Database:
1911: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1913: Level: advanced
1915: Fortran Notes:
1916: Not supported from Fortran
1918: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1919: @*/
1920: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1921: {
1923: PetscInt m,n;
1924: Mat mf,Dmf,Dmat,Ddiff;
1925: PetscReal Diffnorm,Dmfnorm;
1926: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1930: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1931: MatGetLocalSize(mat,&m,&n);
1932: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1933: MatMFFDSetFunction(mf,f,ctx);
1934: MatMFFDSetBase(mf,base,NULL);
1936: MatComputeOperator(mf,MATAIJ,&Dmf);
1937: MatComputeOperator(mat,MATAIJ,&Dmat);
1939: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1940: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1941: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1942: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1943: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1944: flag = PETSC_FALSE;
1945: if (v) {
1946: 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));
1947: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1948: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1949: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1950: }
1951: } else if (v) {
1952: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1953: }
1954: if (flg) *flg = flag;
1955: MatDestroy(&Ddiff);
1956: MatDestroy(&mf);
1957: MatDestroy(&Dmf);
1958: MatDestroy(&Dmat);
1959: return(0);
1960: }
1962: /*@C
1963: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
1965: Logically Collective on Mat
1967: Input Parameters:
1968: + mat - the shell matrix
1969: . f - the function
1970: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1971: - ctx - an optional context for the function
1973: Output Parameter:
1974: . flg - PETSC_TRUE if the multiply is likely correct
1976: Options Database:
1977: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1979: Level: advanced
1981: Fortran Notes:
1982: Not supported from Fortran
1984: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1985: @*/
1986: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1987: {
1989: Vec x,y,z;
1990: PetscInt m,n,M,N;
1991: Mat mf,Dmf,Dmat,Ddiff;
1992: PetscReal Diffnorm,Dmfnorm;
1993: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1997: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1998: MatCreateVecs(mat,&x,&y);
1999: VecDuplicate(y,&z);
2000: MatGetLocalSize(mat,&m,&n);
2001: MatGetSize(mat,&M,&N);
2002: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
2003: MatMFFDSetFunction(mf,f,ctx);
2004: MatMFFDSetBase(mf,base,NULL);
2005: MatComputeOperator(mf,MATAIJ,&Dmf);
2006: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
2007: MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);
2009: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
2010: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
2011: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
2012: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
2013: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
2014: flag = PETSC_FALSE;
2015: if (v) {
2016: 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));
2017: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2018: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2019: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2020: }
2021: } else if (v) {
2022: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
2023: }
2024: if (flg) *flg = flag;
2025: MatDestroy(&mf);
2026: MatDestroy(&Dmat);
2027: MatDestroy(&Ddiff);
2028: MatDestroy(&Dmf);
2029: VecDestroy(&x);
2030: VecDestroy(&y);
2031: VecDestroy(&z);
2032: return(0);
2033: }
2035: /*@C
2036: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
2038: Logically Collective on Mat
2040: Input Parameters:
2041: + mat - the shell matrix
2042: . op - the name of the operation
2043: - g - the function that provides the operation.
2045: Level: advanced
2047: Usage:
2048: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
2049: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
2050: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
2052: Notes:
2053: See the file include/petscmat.h for a complete list of matrix
2054: operations, which all have the form MATOP_<OPERATION>, where
2055: <OPERATION> is the name (in all capital letters) of the
2056: user interface routine (e.g., MatMult() -> MATOP_MULT).
2058: All user-provided functions (except for MATOP_DESTROY) should have the same calling
2059: sequence as the usual matrix interface routines, since they
2060: are intended to be accessed via the usual matrix interface
2061: routines, e.g.,
2062: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
2064: In particular each function MUST return an error code of 0 on success and
2065: nonzero on failure.
2067: Within each user-defined routine, the user should call
2068: MatShellGetContext() to obtain the user-defined context that was
2069: set by MatCreateShell().
2071: Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation()
2073: Fortran Notes:
2074: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
2075: generate a matrix. See src/mat/tests/ex120f.F
2077: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
2078: @*/
2079: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
2080: {
2085: PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
2086: return(0);
2087: }
2089: /*@C
2090: MatShellGetOperation - Gets a matrix function for a shell matrix.
2092: Not Collective
2094: Input Parameters:
2095: + mat - the shell matrix
2096: - op - the name of the operation
2098: Output Parameter:
2099: . g - the function that provides the operation.
2101: Level: advanced
2103: Notes:
2104: See the file include/petscmat.h for a complete list of matrix
2105: operations, which all have the form MATOP_<OPERATION>, where
2106: <OPERATION> is the name (in all capital letters) of the
2107: user interface routine (e.g., MatMult() -> MATOP_MULT).
2109: All user-provided functions have the same calling
2110: sequence as the usual matrix interface routines, since they
2111: are intended to be accessed via the usual matrix interface
2112: routines, e.g.,
2113: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
2115: Within each user-defined routine, the user should call
2116: MatShellGetContext() to obtain the user-defined context that was
2117: set by MatCreateShell().
2119: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
2120: @*/
2121: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
2122: {
2127: PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
2128: return(0);
2129: }
2131: /*@
2132: MatIsShell - Inquires if a matrix is derived from MATSHELL
2134: Input Parameter:
2135: . mat - the matrix
2137: Output Parameter:
2138: . flg - the boolean value
2140: Level: developer
2142: Notes: in the future, we should allow the object type name to be changed still using the MatShell data structure for other matrices (i.e. MATTRANSPOSEMAT, MATSCHURCOMPLEMENT etc)
2144: .seealso: MatCreateShell()
2145: @*/
2146: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2147: {
2151: *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2152: return(0);
2153: }