Actual source code: shell.c
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;
65: /*
66: Store and scale values on zeroed rows
67: xx = [x_1, 0], 0 on zeroed columns
68: */
69: static PetscErrorCode MatShellPreZeroRight(Mat A,Vec x,Vec *xx)
70: {
71: Mat_Shell *shell = (Mat_Shell*)A->data;
75: *xx = x;
76: if (shell->zrows) {
77: VecSet(shell->zvals_w,0.0);
78: VecScatterBegin(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
79: VecScatterEnd(shell->zvals_sct_c,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
80: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
81: }
82: if (shell->zcols) {
83: if (!shell->right_work) {
84: MatCreateVecs(A,&shell->right_work,NULL);
85: }
86: VecCopy(x,shell->right_work);
87: VecISSet(shell->right_work,shell->zcols,0.0);
88: *xx = shell->right_work;
89: }
90: return(0);
91: }
93: /* Insert properly diagonally scaled values stored in MatShellPreZeroRight */
94: static PetscErrorCode MatShellPostZeroLeft(Mat A,Vec x)
95: {
96: Mat_Shell *shell = (Mat_Shell*)A->data;
100: if (shell->zrows) {
101: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
102: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,x,INSERT_VALUES,SCATTER_REVERSE);
103: }
104: return(0);
105: }
107: /*
108: Store and scale values on zeroed rows
109: xx = [x_1, 0], 0 on zeroed rows
110: */
111: static PetscErrorCode MatShellPreZeroLeft(Mat A,Vec x,Vec *xx)
112: {
113: Mat_Shell *shell = (Mat_Shell*)A->data;
117: *xx = NULL;
118: if (!shell->zrows) {
119: *xx = x;
120: } else {
121: if (!shell->left_work) {
122: MatCreateVecs(A,NULL,&shell->left_work);
123: }
124: VecCopy(x,shell->left_work);
125: VecSet(shell->zvals_w,0.0);
126: VecScatterBegin(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
127: VecScatterEnd(shell->zvals_sct_r,shell->zvals_w,shell->left_work,INSERT_VALUES,SCATTER_REVERSE);
128: VecScatterBegin(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
129: VecScatterEnd(shell->zvals_sct_r,x,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
130: VecPointwiseMult(shell->zvals_w,shell->zvals_w,shell->zvals);
131: *xx = shell->left_work;
132: }
133: return(0);
134: }
136: /* Zero zero-columns contributions, sum contributions from properly scaled values stored in MatShellPreZeroLeft */
137: static PetscErrorCode MatShellPostZeroRight(Mat A,Vec x)
138: {
139: Mat_Shell *shell = (Mat_Shell*)A->data;
143: if (shell->zcols) {
144: VecISSet(x,shell->zcols,0.0);
145: }
146: if (shell->zrows) {
147: VecScatterBegin(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
148: VecScatterEnd(shell->zvals_sct_c,shell->zvals_w,x,ADD_VALUES,SCATTER_REVERSE);
149: }
150: return(0);
151: }
153: /*
154: xx = diag(left)*x
155: */
156: static PetscErrorCode MatShellPreScaleLeft(Mat A,Vec x,Vec *xx)
157: {
158: Mat_Shell *shell = (Mat_Shell*)A->data;
162: *xx = NULL;
163: if (!shell->left) {
164: *xx = x;
165: } else {
166: if (!shell->left_work) {VecDuplicate(shell->left,&shell->left_work);}
167: VecPointwiseMult(shell->left_work,x,shell->left);
168: *xx = shell->left_work;
169: }
170: return(0);
171: }
173: /*
174: xx = diag(right)*x
175: */
176: static PetscErrorCode MatShellPreScaleRight(Mat A,Vec x,Vec *xx)
177: {
178: Mat_Shell *shell = (Mat_Shell*)A->data;
182: *xx = NULL;
183: if (!shell->right) {
184: *xx = x;
185: } else {
186: if (!shell->right_work) {VecDuplicate(shell->right,&shell->right_work);}
187: VecPointwiseMult(shell->right_work,x,shell->right);
188: *xx = shell->right_work;
189: }
190: return(0);
191: }
193: /*
194: x = diag(left)*x
195: */
196: static PetscErrorCode MatShellPostScaleLeft(Mat A,Vec x)
197: {
198: Mat_Shell *shell = (Mat_Shell*)A->data;
202: if (shell->left) {VecPointwiseMult(x,x,shell->left);}
203: return(0);
204: }
206: /*
207: x = diag(right)*x
208: */
209: static PetscErrorCode MatShellPostScaleRight(Mat A,Vec x)
210: {
211: Mat_Shell *shell = (Mat_Shell*)A->data;
215: if (shell->right) {VecPointwiseMult(x,x,shell->right);}
216: return(0);
217: }
219: /*
220: Y = vscale*Y + diag(dshift)*X + vshift*X
222: On input Y already contains A*x
223: */
224: static PetscErrorCode MatShellShiftAndScale(Mat A,Vec X,Vec Y)
225: {
226: Mat_Shell *shell = (Mat_Shell*)A->data;
230: if (shell->dshift) { /* get arrays because there is no VecPointwiseMultAdd() */
231: PetscInt i,m;
232: const PetscScalar *x,*d;
233: PetscScalar *y;
234: VecGetLocalSize(X,&m);
235: VecGetArrayRead(shell->dshift,&d);
236: VecGetArrayRead(X,&x);
237: VecGetArray(Y,&y);
238: for (i=0; i<m; i++) y[i] = shell->vscale*y[i] + d[i]*x[i];
239: VecRestoreArrayRead(shell->dshift,&d);
240: VecRestoreArrayRead(X,&x);
241: VecRestoreArray(Y,&y);
242: } else {
243: VecScale(Y,shell->vscale);
244: }
245: if (shell->vshift != 0.0) {VecAXPY(Y,shell->vshift,X);} /* if test is for non-square matrices */
246: return(0);
247: }
249: PetscErrorCode MatShellGetContext_Shell(Mat mat,void *ctx)
250: {
252: *(void**)ctx = ((Mat_Shell*)(mat->data))->ctx;
253: return(0);
254: }
256: /*@
257: MatShellGetContext - Returns the user-provided context associated with a shell matrix.
259: Not Collective
261: Input Parameter:
262: . mat - the matrix, should have been created with MatCreateShell()
264: Output Parameter:
265: . ctx - the user provided context
267: Level: advanced
269: Fortran Notes:
270: To use this from Fortran you must write a Fortran interface definition for this
271: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
273: .seealso: MatCreateShell(), MatShellSetOperation(), MatShellSetContext()
274: @*/
275: PetscErrorCode MatShellGetContext(Mat mat,void *ctx)
276: {
282: PetscUseMethod(mat,"MatShellGetContext_C",(Mat,void*),(mat,ctx));
283: return(0);
284: }
286: static PetscErrorCode MatZeroRowsColumns_Local_Shell(Mat mat,PetscInt nr,PetscInt rows[],PetscInt nc,PetscInt cols[],PetscScalar diag,PetscBool rc)
287: {
289: Mat_Shell *shell = (Mat_Shell*)mat->data;
290: Vec x = NULL,b = NULL;
291: IS is1, is2;
292: const PetscInt *ridxs;
293: PetscInt *idxs,*gidxs;
294: PetscInt cum,rst,cst,i;
297: if (!shell->zvals) {
298: MatCreateVecs(mat,NULL,&shell->zvals);
299: }
300: if (!shell->zvals_w) {
301: VecDuplicate(shell->zvals,&shell->zvals_w);
302: }
303: MatGetOwnershipRange(mat,&rst,NULL);
304: MatGetOwnershipRangeColumn(mat,&cst,NULL);
306: /* Expand/create index set of zeroed rows */
307: PetscMalloc1(nr,&idxs);
308: for (i = 0; i < nr; i++) idxs[i] = rows[i] + rst;
309: ISCreateGeneral(PETSC_COMM_SELF,nr,idxs,PETSC_OWN_POINTER,&is1);
310: ISSort(is1);
311: VecISSet(shell->zvals,is1,diag);
312: if (shell->zrows) {
313: ISSum(shell->zrows,is1,&is2);
314: ISDestroy(&shell->zrows);
315: ISDestroy(&is1);
316: shell->zrows = is2;
317: } else shell->zrows = is1;
319: /* Create scatters for diagonal values communications */
320: VecScatterDestroy(&shell->zvals_sct_c);
321: VecScatterDestroy(&shell->zvals_sct_r);
323: /* row scatter: from/to left vector */
324: MatCreateVecs(mat,&x,&b);
325: VecScatterCreate(b,shell->zrows,shell->zvals_w,shell->zrows,&shell->zvals_sct_r);
327: /* col scatter: from right vector to left vector */
328: ISGetIndices(shell->zrows,&ridxs);
329: ISGetLocalSize(shell->zrows,&nr);
330: PetscMalloc1(nr,&gidxs);
331: for (i = 0, cum = 0; i < nr; i++) {
332: if (ridxs[i] >= mat->cmap->N) continue;
333: gidxs[cum] = ridxs[i];
334: cum++;
335: }
336: ISCreateGeneral(PetscObjectComm((PetscObject)mat),cum,gidxs,PETSC_OWN_POINTER,&is1);
337: VecScatterCreate(x,is1,shell->zvals_w,is1,&shell->zvals_sct_c);
338: ISDestroy(&is1);
339: VecDestroy(&x);
340: VecDestroy(&b);
342: /* Expand/create index set of zeroed columns */
343: if (rc) {
344: PetscMalloc1(nc,&idxs);
345: for (i = 0; i < nc; i++) idxs[i] = cols[i] + cst;
346: ISCreateGeneral(PETSC_COMM_SELF,nc,idxs,PETSC_OWN_POINTER,&is1);
347: ISSort(is1);
348: if (shell->zcols) {
349: ISSum(shell->zcols,is1,&is2);
350: ISDestroy(&shell->zcols);
351: ISDestroy(&is1);
352: shell->zcols = is2;
353: } else shell->zcols = is1;
354: }
355: return(0);
356: }
358: static PetscErrorCode MatZeroRows_Shell(Mat mat,PetscInt n,const PetscInt rows[],PetscScalar diag,Vec x,Vec b)
359: {
360: Mat_Shell *shell = (Mat_Shell*)mat->data;
361: PetscInt nr, *lrows;
365: if (x && b) {
366: Vec xt;
367: PetscScalar *vals;
368: PetscInt *gcols,i,st,nl,nc;
370: PetscMalloc1(n,&gcols);
371: for (i = 0, nc = 0; i < n; i++) if (rows[i] < mat->cmap->N) gcols[nc++] = rows[i];
373: MatCreateVecs(mat,&xt,NULL);
374: VecCopy(x,xt);
375: PetscCalloc1(nc,&vals);
376: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
377: PetscFree(vals);
378: VecAssemblyBegin(xt);
379: VecAssemblyEnd(xt);
380: VecAYPX(xt,-1.0,x); /* xt = [0, x2] */
382: VecGetOwnershipRange(xt,&st,NULL);
383: VecGetLocalSize(xt,&nl);
384: VecGetArray(xt,&vals);
385: for (i = 0; i < nl; i++) {
386: PetscInt g = i + st;
387: if (g > mat->rmap->N) continue;
388: if (PetscAbsScalar(vals[i]) == 0.0) continue;
389: VecSetValue(b,g,diag*vals[i],INSERT_VALUES);
390: }
391: VecRestoreArray(xt,&vals);
392: VecAssemblyBegin(b);
393: VecAssemblyEnd(b); /* b = [b1, x2 * diag] */
394: VecDestroy(&xt);
395: PetscFree(gcols);
396: }
397: PetscLayoutMapLocal(mat->rmap,n,rows,&nr,&lrows,NULL);
398: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,0,NULL,diag,PETSC_FALSE);
399: if (shell->axpy) {
400: MatZeroRows(shell->axpy,n,rows,0.0,NULL,NULL);
401: }
402: PetscFree(lrows);
403: return(0);
404: }
406: static PetscErrorCode MatZeroRowsColumns_Shell(Mat mat,PetscInt n,const PetscInt rowscols[],PetscScalar diag,Vec x,Vec b)
407: {
408: Mat_Shell *shell = (Mat_Shell*)mat->data;
409: PetscInt *lrows, *lcols;
410: PetscInt nr, nc;
411: PetscBool congruent;
415: if (x && b) {
416: Vec xt, bt;
417: PetscScalar *vals;
418: PetscInt *grows,*gcols,i,st,nl;
420: PetscMalloc2(n,&grows,n,&gcols);
421: for (i = 0, nr = 0; i < n; i++) if (rowscols[i] < mat->rmap->N) grows[nr++] = rowscols[i];
422: for (i = 0, nc = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) gcols[nc++] = rowscols[i];
423: PetscCalloc1(n,&vals);
425: MatCreateVecs(mat,&xt,&bt);
426: VecCopy(x,xt);
427: VecSetValues(xt,nc,gcols,vals,INSERT_VALUES); /* xt = [x1, 0] */
428: VecAssemblyBegin(xt);
429: VecAssemblyEnd(xt);
430: VecAXPY(xt,-1.0,x); /* xt = [0, -x2] */
431: MatMult(mat,xt,bt); /* bt = [-A12*x2,-A22*x2] */
432: VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* bt = [-A12*x2,0] */
433: VecAssemblyBegin(bt);
434: VecAssemblyEnd(bt);
435: VecAXPY(b,1.0,bt); /* b = [b1 - A12*x2, b2] */
436: VecSetValues(bt,nr,grows,vals,INSERT_VALUES); /* b = [b1 - A12*x2, 0] */
437: VecAssemblyBegin(bt);
438: VecAssemblyEnd(bt);
439: PetscFree(vals);
441: VecGetOwnershipRange(xt,&st,NULL);
442: VecGetLocalSize(xt,&nl);
443: VecGetArray(xt,&vals);
444: for (i = 0; i < nl; i++) {
445: PetscInt g = i + st;
446: if (g > mat->rmap->N) continue;
447: if (PetscAbsScalar(vals[i]) == 0.0) continue;
448: VecSetValue(b,g,-diag*vals[i],INSERT_VALUES);
449: }
450: VecRestoreArray(xt,&vals);
451: VecAssemblyBegin(b);
452: VecAssemblyEnd(b); /* b = [b1 - A12*x2, x2 * diag] */
453: VecDestroy(&xt);
454: VecDestroy(&bt);
455: PetscFree2(grows,gcols);
456: }
457: PetscLayoutMapLocal(mat->rmap,n,rowscols,&nr,&lrows,NULL);
458: MatHasCongruentLayouts(mat,&congruent);
459: if (congruent) {
460: nc = nr;
461: lcols = lrows;
462: } else { /* MatZeroRowsColumns implicitly assumes the rowscols indices are for a square matrix, here we handle a more general case */
463: PetscInt i,nt,*t;
465: PetscMalloc1(n,&t);
466: for (i = 0, nt = 0; i < n; i++) if (rowscols[i] < mat->cmap->N) t[nt++] = rowscols[i];
467: PetscLayoutMapLocal(mat->cmap,nt,t,&nc,&lcols,NULL);
468: PetscFree(t);
469: }
470: MatZeroRowsColumns_Local_Shell(mat,nr,lrows,nc,lcols,diag,PETSC_TRUE);
471: if (!congruent) {
472: PetscFree(lcols);
473: }
474: PetscFree(lrows);
475: if (shell->axpy) {
476: MatZeroRowsColumns(shell->axpy,n,rowscols,0.0,NULL,NULL);
477: }
478: return(0);
479: }
481: PetscErrorCode MatDestroy_Shell(Mat mat)
482: {
483: PetscErrorCode ierr;
484: Mat_Shell *shell = (Mat_Shell*)mat->data;
485: MatShellMatFunctionList matmat;
488: if (shell->ops->destroy) {
489: (*shell->ops->destroy)(mat);
490: }
491: PetscMemzero(shell->ops,sizeof(struct _MatShellOps));
492: VecDestroy(&shell->left);
493: VecDestroy(&shell->right);
494: VecDestroy(&shell->dshift);
495: VecDestroy(&shell->left_work);
496: VecDestroy(&shell->right_work);
497: VecDestroy(&shell->left_add_work);
498: VecDestroy(&shell->right_add_work);
499: VecDestroy(&shell->axpy_left);
500: VecDestroy(&shell->axpy_right);
501: MatDestroy(&shell->axpy);
502: VecDestroy(&shell->zvals_w);
503: VecDestroy(&shell->zvals);
504: VecScatterDestroy(&shell->zvals_sct_c);
505: VecScatterDestroy(&shell->zvals_sct_r);
506: ISDestroy(&shell->zrows);
507: ISDestroy(&shell->zcols);
509: matmat = shell->matmat;
510: while (matmat) {
511: MatShellMatFunctionList next = matmat->next;
513: PetscObjectComposeFunction((PetscObject)mat,matmat->composedname,NULL);
514: PetscFree(matmat->composedname);
515: PetscFree(matmat->resultname);
516: PetscFree(matmat);
517: matmat = next;
518: }
519: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetContext_C",NULL);
520: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetContext_C",NULL);
521: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetVecType_C",NULL);
522: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetManageScalingShifts_C",NULL);
523: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetOperation_C",NULL);
524: PetscObjectComposeFunction((PetscObject)mat,"MatShellGetOperation_C",NULL);
525: PetscObjectComposeFunction((PetscObject)mat,"MatShellSetMatProductOperation_C",NULL);
526: PetscFree(mat->data);
527: return(0);
528: }
530: typedef struct {
531: PetscErrorCode (*numeric)(Mat,Mat,Mat,void*);
532: PetscErrorCode (*destroy)(void*);
533: void *userdata;
534: Mat B;
535: Mat Bt;
536: Mat axpy;
537: } MatMatDataShell;
539: static PetscErrorCode DestroyMatMatDataShell(void *data)
540: {
541: MatMatDataShell *mmdata = (MatMatDataShell *)data;
545: if (mmdata->destroy) {
546: (*mmdata->destroy)(mmdata->userdata);
547: }
548: MatDestroy(&mmdata->B);
549: MatDestroy(&mmdata->Bt);
550: MatDestroy(&mmdata->axpy);
551: PetscFree(mmdata);
552: return(0);
553: }
555: static PetscErrorCode MatProductNumeric_Shell_X(Mat D)
556: {
557: PetscErrorCode ierr;
558: Mat_Product *product;
559: Mat A, B;
560: MatMatDataShell *mdata;
561: PetscScalar zero = 0.0;
564: MatCheckProduct(D,1);
565: product = D->product;
566: if (!product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data empty");
567: A = product->A;
568: B = product->B;
569: mdata = (MatMatDataShell*)product->data;
570: if (mdata->numeric) {
571: Mat_Shell *shell = (Mat_Shell*)A->data;
572: PetscErrorCode (*stashsym)(Mat) = D->ops->productsymbolic;
573: PetscErrorCode (*stashnum)(Mat) = D->ops->productnumeric;
574: PetscBool useBmdata = PETSC_FALSE, newB = PETSC_TRUE;
576: if (shell->managescalingshifts) {
577: if (shell->zcols || shell->zrows) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_SUP,"MatProduct not supported with zeroed rows/columns");
578: if (shell->right || shell->left) {
579: useBmdata = PETSC_TRUE;
580: if (!mdata->B) {
581: MatDuplicate(B,MAT_SHARE_NONZERO_PATTERN,&mdata->B);
582: } else {
583: newB = PETSC_FALSE;
584: }
585: MatCopy(B,mdata->B,SAME_NONZERO_PATTERN);
586: }
587: switch (product->type) {
588: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
589: if (shell->right) {
590: MatDiagonalScale(mdata->B,shell->right,NULL);
591: }
592: break;
593: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
594: if (shell->left) {
595: MatDiagonalScale(mdata->B,shell->left,NULL);
596: }
597: break;
598: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
599: if (shell->right) {
600: MatDiagonalScale(mdata->B,NULL,shell->right);
601: }
602: break;
603: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
604: if (shell->right && shell->left) {
605: PetscBool flg;
607: VecEqual(shell->right,shell->left,&flg);
608: 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);
609: }
610: if (shell->right) {
611: MatDiagonalScale(mdata->B,NULL,shell->right);
612: }
613: break;
614: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
615: if (shell->right && shell->left) {
616: PetscBool flg;
618: VecEqual(shell->right,shell->left,&flg);
619: 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);
620: }
621: if (shell->right) {
622: MatDiagonalScale(mdata->B,shell->right,NULL);
623: }
624: break;
625: 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);
626: }
627: }
628: /* allow the user to call MatMat operations on D */
629: D->product = NULL;
630: D->ops->productsymbolic = NULL;
631: D->ops->productnumeric = NULL;
633: (*mdata->numeric)(A,useBmdata ? mdata->B : B,D,mdata->userdata);
635: /* clear any leftover user data and restore D pointers */
636: MatProductClear(D);
637: D->ops->productsymbolic = stashsym;
638: D->ops->productnumeric = stashnum;
639: D->product = product;
641: if (shell->managescalingshifts) {
642: MatScale(D,shell->vscale);
643: switch (product->type) {
644: case MATPRODUCT_AB: /* s L A R B + v L R B + L D R B */
645: case MATPRODUCT_ABt: /* s L A R B^t + v L R B^t + L D R B^t */
646: if (shell->left) {
647: MatDiagonalScale(D,shell->left,NULL);
648: if (shell->dshift || shell->vshift != zero) {
649: if (!shell->left_work) {MatCreateVecs(A,NULL,&shell->left_work);}
650: if (shell->dshift) {
651: VecCopy(shell->dshift,shell->left_work);
652: VecShift(shell->left_work,shell->vshift);
653: VecPointwiseMult(shell->left_work,shell->left_work,shell->left);
654: } else {
655: VecSet(shell->left_work,shell->vshift);
656: }
657: if (product->type == MATPRODUCT_ABt) {
658: MatReuse reuse = mdata->Bt ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX;
659: MatStructure str = mdata->Bt ? SUBSET_NONZERO_PATTERN : DIFFERENT_NONZERO_PATTERN;
661: MatTranspose(mdata->B,reuse,&mdata->Bt);
662: MatDiagonalScale(mdata->Bt,shell->left_work,NULL);
663: MatAXPY(D,1.0,mdata->Bt,str);
664: } else {
665: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
667: MatDiagonalScale(mdata->B,shell->left_work,NULL);
668: MatAXPY(D,1.0,mdata->B,str);
669: }
670: }
671: }
672: break;
673: case MATPRODUCT_AtB: /* s R A^t L B + v R L B + R D L B */
674: if (shell->right) {
675: MatDiagonalScale(D,shell->right,NULL);
676: if (shell->dshift || shell->vshift != zero) {
677: MatStructure str = newB ? DIFFERENT_NONZERO_PATTERN : SUBSET_NONZERO_PATTERN;
679: if (!shell->right_work) {MatCreateVecs(A,&shell->right_work,NULL);}
680: if (shell->dshift) {
681: VecCopy(shell->dshift,shell->right_work);
682: VecShift(shell->right_work,shell->vshift);
683: VecPointwiseMult(shell->right_work,shell->right_work,shell->right);
684: } else {
685: VecSet(shell->right_work,shell->vshift);
686: }
687: MatDiagonalScale(mdata->B,shell->right_work,NULL);
688: MatAXPY(D,1.0,mdata->B,str);
689: }
690: }
691: break;
692: case MATPRODUCT_PtAP: /* s B^t L A R B + v B^t L R B + B^t L D R B */
693: case MATPRODUCT_RARt: /* s B L A R B^t + v B L R B^t + B L D R B^t */
694: 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);
695: break;
696: 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);
697: }
698: if (shell->axpy && shell->axpy_vscale != zero) {
699: Mat X;
700: PetscObjectState axpy_state;
701: MatStructure str = DIFFERENT_NONZERO_PATTERN; /* not sure it is safe to ever use SUBSET_NONZERO_PATTERN */
703: MatShellGetContext(shell->axpy,&X);
704: PetscObjectStateGet((PetscObject)X,&axpy_state);
705: 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,...)");
706: if (!mdata->axpy) {
707: str = DIFFERENT_NONZERO_PATTERN;
708: MatProductCreate(shell->axpy,B,NULL,&mdata->axpy);
709: MatProductSetType(mdata->axpy,product->type);
710: MatProductSetFromOptions(mdata->axpy);
711: MatProductSymbolic(mdata->axpy);
712: } else { /* May be that shell->axpy has changed */
713: PetscBool flg;
715: MatProductReplaceMats(shell->axpy,B,NULL,mdata->axpy);
716: MatHasOperation(mdata->axpy,MATOP_PRODUCTSYMBOLIC,&flg);
717: if (!flg) {
718: str = DIFFERENT_NONZERO_PATTERN;
719: MatProductSetFromOptions(mdata->axpy);
720: MatProductSymbolic(mdata->axpy);
721: }
722: }
723: MatProductNumeric(mdata->axpy);
724: MatAXPY(D,shell->axpy_vscale,mdata->axpy,str);
725: }
726: }
727: } else SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Missing numeric operation");
728: return(0);
729: }
731: static PetscErrorCode MatProductSymbolic_Shell_X(Mat D)
732: {
733: PetscErrorCode ierr;
734: Mat_Product *product;
735: Mat A,B;
736: MatShellMatFunctionList matmat;
737: Mat_Shell *shell;
738: PetscBool flg;
739: char composedname[256];
740: MatMatDataShell *mdata;
743: MatCheckProduct(D,1);
744: product = D->product;
745: if (product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_PLIB,"Product data not empty");
746: A = product->A;
747: B = product->B;
748: shell = (Mat_Shell*)A->data;
749: matmat = shell->matmat;
750: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
751: while (matmat) {
752: PetscStrcmp(composedname,matmat->composedname,&flg);
753: flg = (PetscBool)(flg && (matmat->ptype == product->type));
754: if (flg) break;
755: matmat = matmat->next;
756: }
757: if (!flg) SETERRQ2(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Composedname \"%s\" for product type %s not found",composedname,MatProductTypes[product->type]);
758: switch (product->type) {
759: case MATPRODUCT_AB:
760: MatSetSizes(D,A->rmap->n,B->cmap->n,A->rmap->N,B->cmap->N);
761: break;
762: case MATPRODUCT_AtB:
763: MatSetSizes(D,A->cmap->n,B->cmap->n,A->cmap->N,B->cmap->N);
764: break;
765: case MATPRODUCT_ABt:
766: MatSetSizes(D,A->rmap->n,B->rmap->n,A->rmap->N,B->rmap->N);
767: break;
768: case MATPRODUCT_RARt:
769: MatSetSizes(D,B->rmap->n,B->rmap->n,B->rmap->N,B->rmap->N);
770: break;
771: case MATPRODUCT_PtAP:
772: MatSetSizes(D,B->cmap->n,B->cmap->n,B->cmap->N,B->cmap->N);
773: break;
774: 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);
775: }
776: /* respect users who passed in a matrix for which resultname is the base type */
777: if (matmat->resultname) {
778: PetscObjectBaseTypeCompare((PetscObject)D,matmat->resultname,&flg);
779: if (!flg) {
780: MatSetType(D,matmat->resultname);
781: }
782: }
783: /* If matrix type was not set or different, we need to reset this pointers */
784: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
785: D->ops->productnumeric = MatProductNumeric_Shell_X;
786: /* attach product data */
787: PetscNew(&mdata);
788: mdata->numeric = matmat->numeric;
789: mdata->destroy = matmat->destroy;
790: if (matmat->symbolic) {
791: (*matmat->symbolic)(A,B,D,&mdata->userdata);
792: } else { /* call general setup if symbolic operation not provided */
793: MatSetUp(D);
794: }
795: if (!D->product) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product disappeared after user symbolic phase");
796: if (D->product->data) SETERRQ(PetscObjectComm((PetscObject)D),PETSC_ERR_COR,"Product data not empty after user symbolic phase");
797: D->product->data = mdata;
798: D->product->destroy = DestroyMatMatDataShell;
799: /* Be sure to reset these pointers if the user did something unexpected */
800: D->ops->productsymbolic = MatProductSymbolic_Shell_X;
801: D->ops->productnumeric = MatProductNumeric_Shell_X;
802: return(0);
803: }
805: static PetscErrorCode MatProductSetFromOptions_Shell_X(Mat D)
806: {
807: PetscErrorCode ierr;
808: Mat_Product *product;
809: Mat A,B;
810: MatShellMatFunctionList matmat;
811: Mat_Shell *shell;
812: PetscBool flg;
813: char composedname[256];
816: MatCheckProduct(D,1);
817: product = D->product;
818: A = product->A;
819: B = product->B;
820: MatIsShell(A,&flg);
821: if (!flg) return(0);
822: shell = (Mat_Shell*)A->data;
823: matmat = shell->matmat;
824: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,((PetscObject)B)->type_name);
825: while (matmat) {
826: PetscStrcmp(composedname,matmat->composedname,&flg);
827: flg = (PetscBool)(flg && (matmat->ptype == product->type));
828: if (flg) break;
829: matmat = matmat->next;
830: }
831: if (flg) { D->ops->productsymbolic = MatProductSymbolic_Shell_X; }
832: else { PetscInfo2(D," symbolic product %s not registered for product type %s\n",composedname,MatProductTypes[product->type]); }
833: return(0);
834: }
836: 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)
837: {
838: PetscBool flg;
839: PetscErrorCode ierr;
840: Mat_Shell *shell;
841: MatShellMatFunctionList matmat;
844: if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing numeric routine");
845: if (!composedname) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_PLIB,"Missing composed name");
847: /* add product callback */
848: shell = (Mat_Shell*)A->data;
849: matmat = shell->matmat;
850: if (!matmat) {
851: PetscNew(&shell->matmat);
852: matmat = shell->matmat;
853: } else {
854: MatShellMatFunctionList entry = matmat;
855: while (entry) {
856: PetscStrcmp(composedname,entry->composedname,&flg);
857: flg = (PetscBool)(flg && (entry->ptype == ptype));
858: if (flg) goto set;
859: matmat = entry;
860: entry = entry->next;
861: }
862: PetscNew(&matmat->next);
863: matmat = matmat->next;
864: }
866: set:
867: matmat->symbolic = symbolic;
868: matmat->numeric = numeric;
869: matmat->destroy = destroy;
870: matmat->ptype = ptype;
871: PetscFree(matmat->composedname);
872: PetscFree(matmat->resultname);
873: PetscStrallocpy(composedname,&matmat->composedname);
874: PetscStrallocpy(resultname,&matmat->resultname);
875: PetscInfo3(A,"Composing %s for product type %s with result %s\n",matmat->composedname,MatProductTypes[matmat->ptype],matmat->resultname ? matmat->resultname : "not specified");
876: PetscObjectComposeFunction((PetscObject)A,matmat->composedname,MatProductSetFromOptions_Shell_X);
877: return(0);
878: }
880: /*@C
881: MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a shell matrix.
883: Logically Collective on Mat
885: Input Parameters:
886: + A - the shell matrix
887: . ptype - the product type
888: . symbolic - the function for the symbolic phase (can be NULL)
889: . numeric - the function for the numerical phase
890: . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL)
891: . Btype - the matrix type for the matrix to be multiplied against
892: - Ctype - the matrix type for the result (can be NULL)
894: Level: advanced
896: Usage:
897: $ extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**);
898: $ extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*);
899: $ extern PetscErrorCode userdestroy(void*);
900: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
901: $ MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE);
902: $ [ create B of type SEQAIJ etc..]
903: $ MatProductCreate(A,B,NULL,&C);
904: $ MatProductSetType(C,MATPRODUCT_AB);
905: $ MatProductSetFromOptions(C);
906: $ MatProductSymbolic(C); -> actually runs the user defined symbolic operation
907: $ MatProductNumeric(C); -> actually runs the user defined numeric operation
908: $ [ use C = A*B ]
910: Notes:
911: MATPRODUCT_ABC is not supported yet. Not supported in Fortran.
912: If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL.
913: Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
914: PETSc will take care of calling the user-defined callbacks.
915: It is allowed to specify the same callbacks for different Btype matrix types.
916: The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence.
918: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp()
919: @*/
920: 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)
921: {
927: if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]);
928: if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4");
931: 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));
932: return(0);
933: }
935: 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)
936: {
937: PetscBool flg;
939: char composedname[256];
940: MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
941: PetscMPIInt size;
945: while (Bnames) { /* user passed in the root name */
946: PetscStrcmp(Btype,Bnames->rname,&flg);
947: if (flg) break;
948: Bnames = Bnames->next;
949: }
950: while (Cnames) { /* user passed in the root name */
951: PetscStrcmp(Ctype,Cnames->rname,&flg);
952: if (flg) break;
953: Cnames = Cnames->next;
954: }
955: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
956: Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
957: Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
958: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,Btype);
959: MatShellSetMatProductOperation_Private(A,ptype,symbolic,numeric,destroy,composedname,Ctype);
960: return(0);
961: }
963: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
964: {
965: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
966: PetscErrorCode ierr;
967: PetscBool matflg;
968: MatShellMatFunctionList matmatA;
971: MatIsShell(B,&matflg);
972: if (!matflg) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix %s not derived from MATSHELL",((PetscObject)B)->type_name);
974: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
975: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
977: if (shellA->ops->copy) {
978: (*shellA->ops->copy)(A,B,str);
979: }
980: shellB->vscale = shellA->vscale;
981: shellB->vshift = shellA->vshift;
982: if (shellA->dshift) {
983: if (!shellB->dshift) {
984: VecDuplicate(shellA->dshift,&shellB->dshift);
985: }
986: VecCopy(shellA->dshift,shellB->dshift);
987: } else {
988: VecDestroy(&shellB->dshift);
989: }
990: if (shellA->left) {
991: if (!shellB->left) {
992: VecDuplicate(shellA->left,&shellB->left);
993: }
994: VecCopy(shellA->left,shellB->left);
995: } else {
996: VecDestroy(&shellB->left);
997: }
998: if (shellA->right) {
999: if (!shellB->right) {
1000: VecDuplicate(shellA->right,&shellB->right);
1001: }
1002: VecCopy(shellA->right,shellB->right);
1003: } else {
1004: VecDestroy(&shellB->right);
1005: }
1006: MatDestroy(&shellB->axpy);
1007: shellB->axpy_vscale = 0.0;
1008: shellB->axpy_state = 0;
1009: if (shellA->axpy) {
1010: PetscObjectReference((PetscObject)shellA->axpy);
1011: shellB->axpy = shellA->axpy;
1012: shellB->axpy_vscale = shellA->axpy_vscale;
1013: shellB->axpy_state = shellA->axpy_state;
1014: }
1015: if (shellA->zrows) {
1016: ISDuplicate(shellA->zrows,&shellB->zrows);
1017: if (shellA->zcols) {
1018: ISDuplicate(shellA->zcols,&shellB->zcols);
1019: }
1020: VecDuplicate(shellA->zvals,&shellB->zvals);
1021: VecCopy(shellA->zvals,shellB->zvals);
1022: VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
1023: PetscObjectReference((PetscObject)shellA->zvals_sct_r);
1024: PetscObjectReference((PetscObject)shellA->zvals_sct_c);
1025: shellB->zvals_sct_r = shellA->zvals_sct_r;
1026: shellB->zvals_sct_c = shellA->zvals_sct_c;
1027: }
1029: matmatA = shellA->matmat;
1030: if (matmatA) {
1031: while (matmatA->next) {
1032: MatShellSetMatProductOperation_Private(B,matmatA->ptype,matmatA->symbolic,matmatA->numeric,matmatA->destroy,matmatA->composedname,matmatA->resultname);
1033: matmatA = matmatA->next;
1034: }
1035: }
1036: return(0);
1037: }
1039: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
1040: {
1042: void *ctx;
1045: MatShellGetContext(mat,&ctx);
1046: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
1047: PetscObjectChangeTypeName((PetscObject)(*M),((PetscObject)mat)->type_name);
1048: if (op != MAT_DO_NOT_COPY_VALUES) {
1049: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
1050: }
1051: return(0);
1052: }
1054: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
1055: {
1056: Mat_Shell *shell = (Mat_Shell*)A->data;
1057: PetscErrorCode ierr;
1058: Vec xx;
1059: PetscObjectState instate,outstate;
1062: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
1063: MatShellPreZeroRight(A,x,&xx);
1064: MatShellPreScaleRight(A,xx,&xx);
1065: PetscObjectStateGet((PetscObject)y, &instate);
1066: (*shell->ops->mult)(A,xx,y);
1067: PetscObjectStateGet((PetscObject)y, &outstate);
1068: if (instate == outstate) {
1069: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1070: PetscObjectStateIncrease((PetscObject)y);
1071: }
1072: MatShellShiftAndScale(A,xx,y);
1073: MatShellPostScaleLeft(A,y);
1074: MatShellPostZeroLeft(A,y);
1076: if (shell->axpy) {
1077: Mat X;
1078: PetscObjectState axpy_state;
1080: MatShellGetContext(shell->axpy,&X);
1081: PetscObjectStateGet((PetscObject)X,&axpy_state);
1082: 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,...)");
1084: MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1085: VecCopy(x,shell->axpy_right);
1086: MatMult(shell->axpy,shell->axpy_right,shell->axpy_left);
1087: VecAXPY(y,shell->axpy_vscale,shell->axpy_left);
1088: }
1089: return(0);
1090: }
1092: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
1093: {
1094: Mat_Shell *shell = (Mat_Shell*)A->data;
1098: if (y == z) {
1099: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
1100: MatMult(A,x,shell->right_add_work);
1101: VecAXPY(z,1.0,shell->right_add_work);
1102: } else {
1103: MatMult(A,x,z);
1104: VecAXPY(z,1.0,y);
1105: }
1106: return(0);
1107: }
1109: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
1110: {
1111: Mat_Shell *shell = (Mat_Shell*)A->data;
1112: PetscErrorCode ierr;
1113: Vec xx;
1114: PetscObjectState instate,outstate;
1117: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
1118: MatShellPreZeroLeft(A,x,&xx);
1119: MatShellPreScaleLeft(A,xx,&xx);
1120: PetscObjectStateGet((PetscObject)y, &instate);
1121: (*shell->ops->multtranspose)(A,xx,y);
1122: PetscObjectStateGet((PetscObject)y, &outstate);
1123: if (instate == outstate) {
1124: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1125: PetscObjectStateIncrease((PetscObject)y);
1126: }
1127: MatShellShiftAndScale(A,xx,y);
1128: MatShellPostScaleRight(A,y);
1129: MatShellPostZeroRight(A,y);
1131: if (shell->axpy) {
1132: Mat X;
1133: PetscObjectState axpy_state;
1135: MatShellGetContext(shell->axpy,&X);
1136: PetscObjectStateGet((PetscObject)X,&axpy_state);
1137: 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,...)");
1138: MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1139: VecCopy(x,shell->axpy_left);
1140: MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);
1141: VecAXPY(y,shell->axpy_vscale,shell->axpy_right);
1142: }
1143: return(0);
1144: }
1146: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
1147: {
1148: Mat_Shell *shell = (Mat_Shell*)A->data;
1152: if (y == z) {
1153: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
1154: MatMultTranspose(A,x,shell->left_add_work);
1155: VecAXPY(z,1.0,shell->left_add_work);
1156: } else {
1157: MatMultTranspose(A,x,z);
1158: VecAXPY(z,1.0,y);
1159: }
1160: return(0);
1161: }
1163: /*
1164: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1165: */
1166: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
1167: {
1168: Mat_Shell *shell = (Mat_Shell*)A->data;
1172: if (shell->ops->getdiagonal) {
1173: (*shell->ops->getdiagonal)(A,v);
1174: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
1175: VecScale(v,shell->vscale);
1176: if (shell->dshift) {
1177: VecAXPY(v,1.0,shell->dshift);
1178: }
1179: VecShift(v,shell->vshift);
1180: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
1181: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
1182: if (shell->zrows) {
1183: VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
1184: VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
1185: }
1186: if (shell->axpy) {
1187: Mat X;
1188: PetscObjectState axpy_state;
1190: MatShellGetContext(shell->axpy,&X);
1191: PetscObjectStateGet((PetscObject)X,&axpy_state);
1192: 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,...)");
1193: MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);
1194: MatGetDiagonal(shell->axpy,shell->axpy_left);
1195: VecAXPY(v,shell->axpy_vscale,shell->axpy_left);
1196: }
1197: return(0);
1198: }
1200: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
1201: {
1202: Mat_Shell *shell = (Mat_Shell*)Y->data;
1204: PetscBool flg;
1207: MatHasCongruentLayouts(Y,&flg);
1208: if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
1209: if (shell->left || shell->right) {
1210: if (!shell->dshift) {
1211: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
1212: VecSet(shell->dshift,a);
1213: } else {
1214: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
1215: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
1216: VecShift(shell->dshift,a);
1217: }
1218: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
1219: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
1220: } else shell->vshift += a;
1221: if (shell->zrows) {
1222: VecShift(shell->zvals,a);
1223: }
1224: return(0);
1225: }
1227: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
1228: {
1229: Mat_Shell *shell = (Mat_Shell*)A->data;
1233: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
1234: if (shell->left || shell->right) {
1235: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
1236: if (shell->left && shell->right) {
1237: VecPointwiseDivide(shell->right_work,D,shell->left);
1238: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
1239: } else if (shell->left) {
1240: VecPointwiseDivide(shell->right_work,D,shell->left);
1241: } else {
1242: VecPointwiseDivide(shell->right_work,D,shell->right);
1243: }
1244: VecAXPY(shell->dshift,s,shell->right_work);
1245: } else {
1246: VecAXPY(shell->dshift,s,D);
1247: }
1248: return(0);
1249: }
1251: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
1252: {
1253: Mat_Shell *shell = (Mat_Shell*)A->data;
1254: Vec d;
1256: PetscBool flg;
1259: MatHasCongruentLayouts(A,&flg);
1260: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
1261: if (ins == INSERT_VALUES) {
1262: if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
1263: VecDuplicate(D,&d);
1264: MatGetDiagonal(A,d);
1265: MatDiagonalSet_Shell_Private(A,d,-1.);
1266: MatDiagonalSet_Shell_Private(A,D,1.);
1267: VecDestroy(&d);
1268: if (shell->zrows) {
1269: VecCopy(D,shell->zvals);
1270: }
1271: } else {
1272: MatDiagonalSet_Shell_Private(A,D,1.);
1273: if (shell->zrows) {
1274: VecAXPY(shell->zvals,1.0,D);
1275: }
1276: }
1277: return(0);
1278: }
1280: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
1281: {
1282: Mat_Shell *shell = (Mat_Shell*)Y->data;
1286: shell->vscale *= a;
1287: shell->vshift *= a;
1288: if (shell->dshift) {
1289: VecScale(shell->dshift,a);
1290: }
1291: shell->axpy_vscale *= a;
1292: if (shell->zrows) {
1293: VecScale(shell->zvals,a);
1294: }
1295: return(0);
1296: }
1298: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
1299: {
1300: Mat_Shell *shell = (Mat_Shell*)Y->data;
1304: if (left) {
1305: if (!shell->left) {
1306: VecDuplicate(left,&shell->left);
1307: VecCopy(left,shell->left);
1308: } else {
1309: VecPointwiseMult(shell->left,shell->left,left);
1310: }
1311: if (shell->zrows) {
1312: VecPointwiseMult(shell->zvals,shell->zvals,left);
1313: }
1314: }
1315: if (right) {
1316: if (!shell->right) {
1317: VecDuplicate(right,&shell->right);
1318: VecCopy(right,shell->right);
1319: } else {
1320: VecPointwiseMult(shell->right,shell->right,right);
1321: }
1322: if (shell->zrows) {
1323: if (!shell->left_work) {
1324: MatCreateVecs(Y,NULL,&shell->left_work);
1325: }
1326: VecSet(shell->zvals_w,1.0);
1327: VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
1328: VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
1329: VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
1330: }
1331: }
1332: if (shell->axpy) {
1333: MatDiagonalScale(shell->axpy,left,right);
1334: }
1335: return(0);
1336: }
1338: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
1339: {
1340: Mat_Shell *shell = (Mat_Shell*)Y->data;
1344: if (t == MAT_FINAL_ASSEMBLY) {
1345: shell->vshift = 0.0;
1346: shell->vscale = 1.0;
1347: shell->axpy_vscale = 0.0;
1348: shell->axpy_state = 0;
1349: VecDestroy(&shell->dshift);
1350: VecDestroy(&shell->left);
1351: VecDestroy(&shell->right);
1352: MatDestroy(&shell->axpy);
1353: VecDestroy(&shell->axpy_left);
1354: VecDestroy(&shell->axpy_right);
1355: VecScatterDestroy(&shell->zvals_sct_c);
1356: VecScatterDestroy(&shell->zvals_sct_r);
1357: ISDestroy(&shell->zrows);
1358: ISDestroy(&shell->zcols);
1359: }
1360: return(0);
1361: }
1363: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
1364: {
1366: *missing = PETSC_FALSE;
1367: return(0);
1368: }
1370: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
1371: {
1372: Mat_Shell *shell = (Mat_Shell*)Y->data;
1376: if (X == Y) {
1377: MatScale(Y,1.0 + a);
1378: return(0);
1379: }
1380: if (!shell->axpy) {
1381: MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);
1382: shell->axpy_vscale = a;
1383: PetscObjectStateGet((PetscObject)X,&shell->axpy_state);
1384: } else {
1385: MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);
1386: }
1387: return(0);
1388: }
1390: static struct _MatOps MatOps_Values = {NULL,
1391: NULL,
1392: NULL,
1393: NULL,
1394: /* 4*/ MatMultAdd_Shell,
1395: NULL,
1396: MatMultTransposeAdd_Shell,
1397: NULL,
1398: NULL,
1399: NULL,
1400: /*10*/ NULL,
1401: NULL,
1402: NULL,
1403: NULL,
1404: NULL,
1405: /*15*/ NULL,
1406: NULL,
1407: NULL,
1408: MatDiagonalScale_Shell,
1409: NULL,
1410: /*20*/ NULL,
1411: MatAssemblyEnd_Shell,
1412: NULL,
1413: NULL,
1414: /*24*/ MatZeroRows_Shell,
1415: NULL,
1416: NULL,
1417: NULL,
1418: NULL,
1419: /*29*/ NULL,
1420: NULL,
1421: NULL,
1422: NULL,
1423: NULL,
1424: /*34*/ MatDuplicate_Shell,
1425: NULL,
1426: NULL,
1427: NULL,
1428: NULL,
1429: /*39*/ MatAXPY_Shell,
1430: NULL,
1431: NULL,
1432: NULL,
1433: MatCopy_Shell,
1434: /*44*/ NULL,
1435: MatScale_Shell,
1436: MatShift_Shell,
1437: MatDiagonalSet_Shell,
1438: MatZeroRowsColumns_Shell,
1439: /*49*/ NULL,
1440: NULL,
1441: NULL,
1442: NULL,
1443: NULL,
1444: /*54*/ NULL,
1445: NULL,
1446: NULL,
1447: NULL,
1448: NULL,
1449: /*59*/ NULL,
1450: MatDestroy_Shell,
1451: NULL,
1452: MatConvertFrom_Shell,
1453: NULL,
1454: /*64*/ NULL,
1455: NULL,
1456: NULL,
1457: NULL,
1458: NULL,
1459: /*69*/ NULL,
1460: NULL,
1461: MatConvert_Shell,
1462: NULL,
1463: NULL,
1464: /*74*/ NULL,
1465: NULL,
1466: NULL,
1467: NULL,
1468: NULL,
1469: /*79*/ NULL,
1470: NULL,
1471: NULL,
1472: NULL,
1473: NULL,
1474: /*84*/ NULL,
1475: NULL,
1476: NULL,
1477: NULL,
1478: NULL,
1479: /*89*/ NULL,
1480: NULL,
1481: NULL,
1482: NULL,
1483: NULL,
1484: /*94*/ NULL,
1485: NULL,
1486: NULL,
1487: NULL,
1488: NULL,
1489: /*99*/ NULL,
1490: NULL,
1491: NULL,
1492: NULL,
1493: NULL,
1494: /*104*/ NULL,
1495: NULL,
1496: NULL,
1497: NULL,
1498: NULL,
1499: /*109*/ NULL,
1500: NULL,
1501: NULL,
1502: NULL,
1503: MatMissingDiagonal_Shell,
1504: /*114*/ NULL,
1505: NULL,
1506: NULL,
1507: NULL,
1508: NULL,
1509: /*119*/ NULL,
1510: NULL,
1511: NULL,
1512: NULL,
1513: NULL,
1514: /*124*/ NULL,
1515: NULL,
1516: NULL,
1517: NULL,
1518: NULL,
1519: /*129*/ NULL,
1520: NULL,
1521: NULL,
1522: NULL,
1523: NULL,
1524: /*134*/ NULL,
1525: NULL,
1526: NULL,
1527: NULL,
1528: NULL,
1529: /*139*/ NULL,
1530: NULL,
1531: NULL
1532: };
1534: PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx)
1535: {
1536: Mat_Shell *shell = (Mat_Shell*)mat->data;
1539: shell->ctx = ctx;
1540: return(0);
1541: }
1543: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1544: {
1548: PetscFree(mat->defaultvectype);
1549: PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
1550: return(0);
1551: }
1553: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1554: {
1555: Mat_Shell *shell = (Mat_Shell*)A->data;
1558: shell->managescalingshifts = PETSC_FALSE;
1559: A->ops->diagonalset = NULL;
1560: A->ops->diagonalscale = NULL;
1561: A->ops->scale = NULL;
1562: A->ops->shift = NULL;
1563: A->ops->axpy = NULL;
1564: return(0);
1565: }
1567: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1568: {
1569: Mat_Shell *shell = (Mat_Shell*)mat->data;
1572: switch (op) {
1573: case MATOP_DESTROY:
1574: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1575: break;
1576: case MATOP_VIEW:
1577: if (!mat->ops->viewnative) {
1578: mat->ops->viewnative = mat->ops->view;
1579: }
1580: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1581: break;
1582: case MATOP_COPY:
1583: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1584: break;
1585: case MATOP_DIAGONAL_SET:
1586: case MATOP_DIAGONAL_SCALE:
1587: case MATOP_SHIFT:
1588: case MATOP_SCALE:
1589: case MATOP_AXPY:
1590: case MATOP_ZERO_ROWS:
1591: case MATOP_ZERO_ROWS_COLUMNS:
1592: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1593: (((void(**)(void))mat->ops)[op]) = f;
1594: break;
1595: case MATOP_GET_DIAGONAL:
1596: if (shell->managescalingshifts) {
1597: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1598: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1599: } else {
1600: shell->ops->getdiagonal = NULL;
1601: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1602: }
1603: break;
1604: case MATOP_MULT:
1605: if (shell->managescalingshifts) {
1606: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1607: mat->ops->mult = MatMult_Shell;
1608: } else {
1609: shell->ops->mult = NULL;
1610: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1611: }
1612: break;
1613: case MATOP_MULT_TRANSPOSE:
1614: if (shell->managescalingshifts) {
1615: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1616: mat->ops->multtranspose = MatMultTranspose_Shell;
1617: } else {
1618: shell->ops->multtranspose = NULL;
1619: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1620: }
1621: break;
1622: default:
1623: (((void(**)(void))mat->ops)[op]) = f;
1624: break;
1625: }
1626: return(0);
1627: }
1629: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1630: {
1631: Mat_Shell *shell = (Mat_Shell*)mat->data;
1634: switch (op) {
1635: case MATOP_DESTROY:
1636: *f = (void (*)(void))shell->ops->destroy;
1637: break;
1638: case MATOP_VIEW:
1639: *f = (void (*)(void))mat->ops->view;
1640: break;
1641: case MATOP_COPY:
1642: *f = (void (*)(void))shell->ops->copy;
1643: break;
1644: case MATOP_DIAGONAL_SET:
1645: case MATOP_DIAGONAL_SCALE:
1646: case MATOP_SHIFT:
1647: case MATOP_SCALE:
1648: case MATOP_AXPY:
1649: case MATOP_ZERO_ROWS:
1650: case MATOP_ZERO_ROWS_COLUMNS:
1651: *f = (((void (**)(void))mat->ops)[op]);
1652: break;
1653: case MATOP_GET_DIAGONAL:
1654: if (shell->ops->getdiagonal)
1655: *f = (void (*)(void))shell->ops->getdiagonal;
1656: else
1657: *f = (((void (**)(void))mat->ops)[op]);
1658: break;
1659: case MATOP_MULT:
1660: if (shell->ops->mult)
1661: *f = (void (*)(void))shell->ops->mult;
1662: else
1663: *f = (((void (**)(void))mat->ops)[op]);
1664: break;
1665: case MATOP_MULT_TRANSPOSE:
1666: if (shell->ops->multtranspose)
1667: *f = (void (*)(void))shell->ops->multtranspose;
1668: else
1669: *f = (((void (**)(void))mat->ops)[op]);
1670: break;
1671: default:
1672: *f = (((void (**)(void))mat->ops)[op]);
1673: }
1674: return(0);
1675: }
1677: /*MC
1678: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
1680: Level: advanced
1682: .seealso: MatCreateShell()
1683: M*/
1685: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1686: {
1687: Mat_Shell *b;
1691: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1693: PetscNewLog(A,&b);
1694: A->data = (void*)b;
1696: b->ctx = NULL;
1697: b->vshift = 0.0;
1698: b->vscale = 1.0;
1699: b->managescalingshifts = PETSC_TRUE;
1700: A->assembled = PETSC_TRUE;
1701: A->preallocated = PETSC_FALSE;
1703: PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1704: PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1705: PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
1706: PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1707: PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1708: PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1709: PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell);
1710: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1711: return(0);
1712: }
1714: /*@C
1715: MatCreateShell - Creates a new matrix class for use with a user-defined
1716: private data storage format.
1718: Collective
1720: Input Parameters:
1721: + comm - MPI communicator
1722: . m - number of local rows (must be given)
1723: . n - number of local columns (must be given)
1724: . M - number of global rows (may be PETSC_DETERMINE)
1725: . N - number of global columns (may be PETSC_DETERMINE)
1726: - ctx - pointer to data needed by the shell matrix routines
1728: Output Parameter:
1729: . A - the matrix
1731: Level: advanced
1733: Usage:
1734: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1735: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
1736: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1737: $ [ Use matrix for operations that have been set ]
1738: $ MatDestroy(mat);
1740: Notes:
1741: The shell matrix type is intended to provide a simple class to use
1742: with KSP (such as, for use with matrix-free methods). You should not
1743: use the shell type if you plan to define a complete matrix class.
1745: Fortran Notes:
1746: To use this from Fortran with a ctx you must write an interface definition for this
1747: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
1748: in as the ctx argument.
1750: PETSc requires that matrices and vectors being used for certain
1751: operations are partitioned accordingly. For example, when
1752: creating a shell matrix, A, that supports parallel matrix-vector
1753: products using MatMult(A,x,y) the user should set the number
1754: of local matrix rows to be the number of local elements of the
1755: corresponding result vector, y. Note that this is information is
1756: required for use of the matrix interface routines, even though
1757: the shell matrix may not actually be physically partitioned.
1758: For example,
1760: $
1761: $ Vec x, y
1762: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1763: $ Mat A
1764: $
1765: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1766: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1767: $ VecGetLocalSize(y,&m);
1768: $ VecGetLocalSize(x,&n);
1769: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1770: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1771: $ MatMult(A,x,y);
1772: $ MatDestroy(&A);
1773: $ VecDestroy(&y);
1774: $ VecDestroy(&x);
1775: $
1777: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these
1778: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
1780: For rectangular matrices do all the scalings and shifts make sense?
1782: Developers Notes:
1783: Regarding shifting and scaling. The general form is
1785: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1787: The order you apply the operations is important. For example if you have a dshift then
1788: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1789: you get s*vscale*A + diag(shift)
1791: A is the user provided function.
1793: KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
1794: for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
1795: an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
1796: each time the MATSHELL matrix has changed.
1798: Matrix product operations (i.e. MatMat, MatTranposeMat etc) can be specified using MatShellSetMatProductOperation()
1800: Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
1801: with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().
1803: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
1804: @*/
1805: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1806: {
1810: MatCreate(comm,A);
1811: MatSetSizes(*A,m,n,M,N);
1812: MatSetType(*A,MATSHELL);
1813: MatShellSetContext(*A,ctx);
1814: MatSetUp(*A);
1815: return(0);
1816: }
1818: /*@
1819: MatShellSetContext - sets the context for a shell matrix
1821: Logically Collective on Mat
1823: Input Parameters:
1824: + mat - the shell matrix
1825: - ctx - the context
1827: Level: advanced
1829: Fortran Notes:
1830: To use this from Fortran you must write a Fortran interface definition for this
1831: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1833: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1834: @*/
1835: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
1836: {
1841: PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1842: return(0);
1843: }
1845: /*@C
1846: MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs()
1848: Logically collective
1850: Input Parameters:
1851: + mat - the shell matrix
1852: - vtype - type to use for creating vectors
1854: Notes:
1856: Level: advanced
1858: .seealso: MatCreateVecs()
1859: @*/
1860: PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype)
1861: {
1865: PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1866: return(0);
1867: }
1869: /*@
1870: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
1871: after MatCreateShell()
1873: Logically Collective on Mat
1875: Input Parameter:
1876: . mat - the shell matrix
1878: Level: advanced
1880: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1881: @*/
1882: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1883: {
1888: PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1889: return(0);
1890: }
1892: /*@C
1893: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
1895: Logically Collective on Mat
1897: Input Parameters:
1898: + mat - the shell matrix
1899: . f - the function
1900: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1901: - ctx - an optional context for the function
1903: Output Parameter:
1904: . flg - PETSC_TRUE if the multiply is likely correct
1906: Options Database:
1907: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1909: Level: advanced
1911: Fortran Notes:
1912: Not supported from Fortran
1914: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1915: @*/
1916: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1917: {
1919: PetscInt m,n;
1920: Mat mf,Dmf,Dmat,Ddiff;
1921: PetscReal Diffnorm,Dmfnorm;
1922: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1926: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1927: MatGetLocalSize(mat,&m,&n);
1928: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1929: MatMFFDSetFunction(mf,f,ctx);
1930: MatMFFDSetBase(mf,base,NULL);
1932: MatComputeOperator(mf,MATAIJ,&Dmf);
1933: MatComputeOperator(mat,MATAIJ,&Dmat);
1935: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1936: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1937: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1938: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1939: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1940: flag = PETSC_FALSE;
1941: if (v) {
1942: 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));
1943: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1944: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1945: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1946: }
1947: } else if (v) {
1948: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1949: }
1950: if (flg) *flg = flag;
1951: MatDestroy(&Ddiff);
1952: MatDestroy(&mf);
1953: MatDestroy(&Dmf);
1954: MatDestroy(&Dmat);
1955: return(0);
1956: }
1958: /*@C
1959: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
1961: Logically Collective on Mat
1963: Input Parameters:
1964: + mat - the shell matrix
1965: . f - the function
1966: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1967: - ctx - an optional context for the function
1969: Output Parameter:
1970: . flg - PETSC_TRUE if the multiply is likely correct
1972: Options Database:
1973: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1975: Level: advanced
1977: Fortran Notes:
1978: Not supported from Fortran
1980: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1981: @*/
1982: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1983: {
1985: Vec x,y,z;
1986: PetscInt m,n,M,N;
1987: Mat mf,Dmf,Dmat,Ddiff;
1988: PetscReal Diffnorm,Dmfnorm;
1989: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1993: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1994: MatCreateVecs(mat,&x,&y);
1995: VecDuplicate(y,&z);
1996: MatGetLocalSize(mat,&m,&n);
1997: MatGetSize(mat,&M,&N);
1998: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
1999: MatMFFDSetFunction(mf,f,ctx);
2000: MatMFFDSetBase(mf,base,NULL);
2001: MatComputeOperator(mf,MATAIJ,&Dmf);
2002: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
2003: MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);
2005: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
2006: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
2007: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
2008: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
2009: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
2010: flag = PETSC_FALSE;
2011: if (v) {
2012: 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));
2013: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2014: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2015: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2016: }
2017: } else if (v) {
2018: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
2019: }
2020: if (flg) *flg = flag;
2021: MatDestroy(&mf);
2022: MatDestroy(&Dmat);
2023: MatDestroy(&Ddiff);
2024: MatDestroy(&Dmf);
2025: VecDestroy(&x);
2026: VecDestroy(&y);
2027: VecDestroy(&z);
2028: return(0);
2029: }
2031: /*@C
2032: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
2034: Logically Collective on Mat
2036: Input Parameters:
2037: + mat - the shell matrix
2038: . op - the name of the operation
2039: - g - the function that provides the operation.
2041: Level: advanced
2043: Usage:
2044: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
2045: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
2046: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
2048: Notes:
2049: See the file include/petscmat.h for a complete list of matrix
2050: operations, which all have the form MATOP_<OPERATION>, where
2051: <OPERATION> is the name (in all capital letters) of the
2052: user interface routine (e.g., MatMult() -> MATOP_MULT).
2054: All user-provided functions (except for MATOP_DESTROY) should have the same calling
2055: sequence as the usual matrix interface routines, since they
2056: are intended to be accessed via the usual matrix interface
2057: routines, e.g.,
2058: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
2060: In particular each function MUST return an error code of 0 on success and
2061: nonzero on failure.
2063: Within each user-defined routine, the user should call
2064: MatShellGetContext() to obtain the user-defined context that was
2065: set by MatCreateShell().
2067: Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation()
2069: Fortran Notes:
2070: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
2071: generate a matrix. See src/mat/tests/ex120f.F
2073: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
2074: @*/
2075: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
2076: {
2081: PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
2082: return(0);
2083: }
2085: /*@C
2086: MatShellGetOperation - Gets a matrix function for a shell matrix.
2088: Not Collective
2090: Input Parameters:
2091: + mat - the shell matrix
2092: - op - the name of the operation
2094: Output Parameter:
2095: . g - the function that provides the operation.
2097: Level: advanced
2099: Notes:
2100: See the file include/petscmat.h for a complete list of matrix
2101: operations, which all have the form MATOP_<OPERATION>, where
2102: <OPERATION> is the name (in all capital letters) of the
2103: user interface routine (e.g., MatMult() -> MATOP_MULT).
2105: All user-provided functions have the same calling
2106: sequence as the usual matrix interface routines, since they
2107: are intended to be accessed via the usual matrix interface
2108: routines, e.g.,
2109: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
2111: Within each user-defined routine, the user should call
2112: MatShellGetContext() to obtain the user-defined context that was
2113: set by MatCreateShell().
2115: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
2116: @*/
2117: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
2118: {
2123: PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
2124: return(0);
2125: }
2127: /*@
2128: MatIsShell - Inquires if a matrix is derived from MATSHELL
2130: Input Parameter:
2131: . mat - the matrix
2133: Output Parameter:
2134: . flg - the boolean value
2136: Level: developer
2138: 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)
2140: .seealso: MatCreateShell()
2141: @*/
2142: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2143: {
2147: *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2148: return(0);
2149: }