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;
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) goto set;
860: matmat = entry;
861: entry = entry->next;
862: }
863: PetscNew(&matmat->next);
864: matmat = matmat->next;
865: }
867: set:
868: matmat->symbolic = symbolic;
869: matmat->numeric = numeric;
870: matmat->destroy = destroy;
871: matmat->ptype = ptype;
872: PetscFree(matmat->composedname);
873: PetscFree(matmat->resultname);
874: PetscStrallocpy(composedname,&matmat->composedname);
875: PetscStrallocpy(resultname,&matmat->resultname);
876: PetscInfo3(A,"Composing %s for product type %s with result %s\n",matmat->composedname,MatProductTypes[matmat->ptype],matmat->resultname ? matmat->resultname : "not specified");
877: PetscObjectComposeFunction((PetscObject)A,matmat->composedname,MatProductSetFromOptions_Shell_X);
878: return(0);
879: }
881: /*@C
882: MatShellSetMatProductOperation - Allows user to set a matrix matrix operation for a shell matrix.
884: Logically Collective on Mat
886: Input Parameters:
887: + A - the shell matrix
888: . ptype - the product type
889: . symbolic - the function for the symbolic phase (can be NULL)
890: . numeric - the function for the numerical phase
891: . destroy - the function for the destruction of the needed data generated during the symbolic phase (can be NULL)
892: . Btype - the matrix type for the matrix to be multiplied against
893: - Ctype - the matrix type for the result (can be NULL)
895: Level: advanced
897: Usage:
898: $ extern PetscErrorCode usersymbolic(Mat,Mat,Mat,void**);
899: $ extern PetscErrorCode usernumeric(Mat,Mat,Mat,void*);
900: $ extern PetscErrorCode userdestroy(void*);
901: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
902: $ MatShellSetMatProductOperation(A,MATPRODUCT_AB,usersymbolic,usernumeric,userdestroy,MATSEQAIJ,MATDENSE);
903: $ [ create B of type SEQAIJ etc..]
904: $ MatProductCreate(A,B,NULL,&C);
905: $ MatProductSetType(C,MATPRODUCT_AB);
906: $ MatProductSetFromOptions(C);
907: $ MatProductSymbolic(C); -> actually runs the user defined symbolic operation
908: $ MatProductNumeric(C); -> actually runs the user defined numeric operation
909: $ [ use C = A*B ]
911: Notes:
912: MATPRODUCT_ABC is not supported yet. Not supported in Fortran.
913: If the symbolic phase is not specified, MatSetUp() is called on the result matrix that must have its type set if Ctype is NULL.
914: Any additional data needed by the matrix product needs to be returned during the symbolic phase and destroyed with the destroy callback.
915: PETSc will take care of calling the user-defined callbacks.
916: It is allowed to specify the same callbacks for different Btype matrix types.
917: The couple (Btype,ptype) uniquely identifies the operation: the last specified callbacks takes precedence.
919: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatProductType, MatType, MatSetUp()
920: @*/
921: 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)
922: {
928: if (ptype == MATPRODUCT_ABC) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Not for product type %s",MatProductTypes[ptype]);
929: if (!numeric) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_USER,"Missing numeric routine, argument 4");
932: 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));
933: return(0);
934: }
936: 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)
937: {
938: PetscBool flg;
940: char composedname[256];
941: MatRootName Bnames = MatRootNameList, Cnames = MatRootNameList;
942: PetscMPIInt size;
946: while (Bnames) { /* user passed in the root name */
947: PetscStrcmp(Btype,Bnames->rname,&flg);
948: if (flg) break;
949: Bnames = Bnames->next;
950: }
951: while (Cnames) { /* user passed in the root name */
952: PetscStrcmp(Ctype,Cnames->rname,&flg);
953: if (flg) break;
954: Cnames = Cnames->next;
955: }
956: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
957: Btype = Bnames ? (size > 1 ? Bnames->mname : Bnames->sname) : Btype;
958: Ctype = Cnames ? (size > 1 ? Cnames->mname : Cnames->sname) : Ctype;
959: PetscSNPrintf(composedname,sizeof(composedname),"MatProductSetFromOptions_%s_%s_C",((PetscObject)A)->type_name,Btype);
960: MatShellSetMatProductOperation_Private(A,ptype,symbolic,numeric,destroy,composedname,Ctype);
961: return(0);
962: }
964: PetscErrorCode MatCopy_Shell(Mat A,Mat B,MatStructure str)
965: {
966: Mat_Shell *shellA = (Mat_Shell*)A->data,*shellB = (Mat_Shell*)B->data;
967: PetscErrorCode ierr;
968: PetscBool matflg;
969: MatShellMatFunctionList matmatA;
972: MatIsShell(B,&matflg);
973: if (!matflg) SETERRQ1(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix %s not derived from MATSHELL",((PetscObject)B)->type_name);
975: PetscMemcpy(B->ops,A->ops,sizeof(struct _MatOps));
976: PetscMemcpy(shellB->ops,shellA->ops,sizeof(struct _MatShellOps));
978: if (shellA->ops->copy) {
979: (*shellA->ops->copy)(A,B,str);
980: }
981: shellB->vscale = shellA->vscale;
982: shellB->vshift = shellA->vshift;
983: if (shellA->dshift) {
984: if (!shellB->dshift) {
985: VecDuplicate(shellA->dshift,&shellB->dshift);
986: }
987: VecCopy(shellA->dshift,shellB->dshift);
988: } else {
989: VecDestroy(&shellB->dshift);
990: }
991: if (shellA->left) {
992: if (!shellB->left) {
993: VecDuplicate(shellA->left,&shellB->left);
994: }
995: VecCopy(shellA->left,shellB->left);
996: } else {
997: VecDestroy(&shellB->left);
998: }
999: if (shellA->right) {
1000: if (!shellB->right) {
1001: VecDuplicate(shellA->right,&shellB->right);
1002: }
1003: VecCopy(shellA->right,shellB->right);
1004: } else {
1005: VecDestroy(&shellB->right);
1006: }
1007: MatDestroy(&shellB->axpy);
1008: shellB->axpy_vscale = 0.0;
1009: shellB->axpy_state = 0;
1010: if (shellA->axpy) {
1011: PetscObjectReference((PetscObject)shellA->axpy);
1012: shellB->axpy = shellA->axpy;
1013: shellB->axpy_vscale = shellA->axpy_vscale;
1014: shellB->axpy_state = shellA->axpy_state;
1015: }
1016: if (shellA->zrows) {
1017: ISDuplicate(shellA->zrows,&shellB->zrows);
1018: if (shellA->zcols) {
1019: ISDuplicate(shellA->zcols,&shellB->zcols);
1020: }
1021: VecDuplicate(shellA->zvals,&shellB->zvals);
1022: VecCopy(shellA->zvals,shellB->zvals);
1023: VecDuplicate(shellA->zvals_w,&shellB->zvals_w);
1024: PetscObjectReference((PetscObject)shellA->zvals_sct_r);
1025: PetscObjectReference((PetscObject)shellA->zvals_sct_c);
1026: shellB->zvals_sct_r = shellA->zvals_sct_r;
1027: shellB->zvals_sct_c = shellA->zvals_sct_c;
1028: }
1030: matmatA = shellA->matmat;
1031: if (matmatA) {
1032: while (matmatA->next) {
1033: MatShellSetMatProductOperation_Private(B,matmatA->ptype,matmatA->symbolic,matmatA->numeric,matmatA->destroy,matmatA->composedname,matmatA->resultname);
1034: matmatA = matmatA->next;
1035: }
1036: }
1037: return(0);
1038: }
1040: PetscErrorCode MatDuplicate_Shell(Mat mat,MatDuplicateOption op,Mat *M)
1041: {
1043: void *ctx;
1046: MatShellGetContext(mat,&ctx);
1047: MatCreateShell(PetscObjectComm((PetscObject)mat),mat->rmap->n,mat->cmap->n,mat->rmap->N,mat->cmap->N,ctx,M);
1048: PetscObjectChangeTypeName((PetscObject)(*M),((PetscObject)mat)->type_name);
1049: if (op != MAT_DO_NOT_COPY_VALUES) {
1050: MatCopy(mat,*M,SAME_NONZERO_PATTERN);
1051: }
1052: return(0);
1053: }
1055: PetscErrorCode MatMult_Shell(Mat A,Vec x,Vec y)
1056: {
1057: Mat_Shell *shell = (Mat_Shell*)A->data;
1058: PetscErrorCode ierr;
1059: Vec xx;
1060: PetscObjectState instate,outstate;
1063: if (!shell->ops->mult) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMult() for this MATSHELL");
1064: MatShellPreZeroRight(A,x,&xx);
1065: MatShellPreScaleRight(A,xx,&xx);
1066: PetscObjectStateGet((PetscObject)y, &instate);
1067: (*shell->ops->mult)(A,xx,y);
1068: PetscObjectStateGet((PetscObject)y, &outstate);
1069: if (instate == outstate) {
1070: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1071: PetscObjectStateIncrease((PetscObject)y);
1072: }
1073: MatShellShiftAndScale(A,xx,y);
1074: MatShellPostScaleLeft(A,y);
1075: MatShellPostZeroLeft(A,y);
1077: if (shell->axpy) {
1078: Mat X;
1079: PetscObjectState axpy_state;
1081: MatShellGetContext(shell->axpy,(void *)&X);
1082: PetscObjectStateGet((PetscObject)X,&axpy_state);
1083: 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,...)");
1085: MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1086: VecCopy(x,shell->axpy_right);
1087: MatMult(shell->axpy,shell->axpy_right,shell->axpy_left);
1088: VecAXPY(y,shell->axpy_vscale,shell->axpy_left);
1089: }
1090: return(0);
1091: }
1093: PetscErrorCode MatMultAdd_Shell(Mat A,Vec x,Vec y,Vec z)
1094: {
1095: Mat_Shell *shell = (Mat_Shell*)A->data;
1099: if (y == z) {
1100: if (!shell->right_add_work) {VecDuplicate(z,&shell->right_add_work);}
1101: MatMult(A,x,shell->right_add_work);
1102: VecAXPY(z,1.0,shell->right_add_work);
1103: } else {
1104: MatMult(A,x,z);
1105: VecAXPY(z,1.0,y);
1106: }
1107: return(0);
1108: }
1110: PetscErrorCode MatMultTranspose_Shell(Mat A,Vec x,Vec y)
1111: {
1112: Mat_Shell *shell = (Mat_Shell*)A->data;
1113: PetscErrorCode ierr;
1114: Vec xx;
1115: PetscObjectState instate,outstate;
1118: if (!shell->ops->multtranspose) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Have not provided a MatMultTranspose() for this MATSHELL");
1119: MatShellPreZeroLeft(A,x,&xx);
1120: MatShellPreScaleLeft(A,xx,&xx);
1121: PetscObjectStateGet((PetscObject)y, &instate);
1122: (*shell->ops->multtranspose)(A,xx,y);
1123: PetscObjectStateGet((PetscObject)y, &outstate);
1124: if (instate == outstate) {
1125: /* increase the state of the output vector since the user did not update its state themself as should have been done */
1126: PetscObjectStateIncrease((PetscObject)y);
1127: }
1128: MatShellShiftAndScale(A,xx,y);
1129: MatShellPostScaleRight(A,y);
1130: MatShellPostZeroRight(A,y);
1132: if (shell->axpy) {
1133: Mat X;
1134: PetscObjectState axpy_state;
1136: MatShellGetContext(shell->axpy,(void *)&X);
1137: PetscObjectStateGet((PetscObject)X,&axpy_state);
1138: 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,...)");
1139: MatCreateVecs(shell->axpy,shell->axpy_right ? NULL : &shell->axpy_right,shell->axpy_left ? NULL : &shell->axpy_left);
1140: VecCopy(x,shell->axpy_left);
1141: MatMultTranspose(shell->axpy,shell->axpy_left,shell->axpy_right);
1142: VecAXPY(y,shell->axpy_vscale,shell->axpy_right);
1143: }
1144: return(0);
1145: }
1147: PetscErrorCode MatMultTransposeAdd_Shell(Mat A,Vec x,Vec y,Vec z)
1148: {
1149: Mat_Shell *shell = (Mat_Shell*)A->data;
1153: if (y == z) {
1154: if (!shell->left_add_work) {VecDuplicate(z,&shell->left_add_work);}
1155: MatMultTranspose(A,x,shell->left_add_work);
1156: VecAXPY(z,1.0,shell->left_add_work);
1157: } else {
1158: MatMultTranspose(A,x,z);
1159: VecAXPY(z,1.0,y);
1160: }
1161: return(0);
1162: }
1164: /*
1165: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1166: */
1167: PetscErrorCode MatGetDiagonal_Shell(Mat A,Vec v)
1168: {
1169: Mat_Shell *shell = (Mat_Shell*)A->data;
1173: if (shell->ops->getdiagonal) {
1174: (*shell->ops->getdiagonal)(A,v);
1175: } else SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONGSTATE,"Must provide shell matrix with routine to return diagonal using\nMatShellSetOperation(S,MATOP_GET_DIAGONAL,...)");
1176: VecScale(v,shell->vscale);
1177: if (shell->dshift) {
1178: VecAXPY(v,1.0,shell->dshift);
1179: }
1180: VecShift(v,shell->vshift);
1181: if (shell->left) {VecPointwiseMult(v,v,shell->left);}
1182: if (shell->right) {VecPointwiseMult(v,v,shell->right);}
1183: if (shell->zrows) {
1184: VecScatterBegin(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
1185: VecScatterEnd(shell->zvals_sct_r,shell->zvals,v,INSERT_VALUES,SCATTER_REVERSE);
1186: }
1187: if (shell->axpy) {
1188: Mat X;
1189: PetscObjectState axpy_state;
1191: MatShellGetContext(shell->axpy,(void *)&X);
1192: PetscObjectStateGet((PetscObject)X,&axpy_state);
1193: 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,...)");
1194: MatCreateVecs(shell->axpy,NULL,shell->axpy_left ? NULL : &shell->axpy_left);
1195: MatGetDiagonal(shell->axpy,shell->axpy_left);
1196: VecAXPY(v,shell->axpy_vscale,shell->axpy_left);
1197: }
1198: return(0);
1199: }
1201: PetscErrorCode MatShift_Shell(Mat Y,PetscScalar a)
1202: {
1203: Mat_Shell *shell = (Mat_Shell*)Y->data;
1205: PetscBool flg;
1208: MatHasCongruentLayouts(Y,&flg);
1209: if (!flg) SETERRQ(PetscObjectComm((PetscObject)Y),PETSC_ERR_SUP,"Cannot shift shell matrix if it is not congruent");
1210: if (shell->left || shell->right) {
1211: if (!shell->dshift) {
1212: VecDuplicate(shell->left ? shell->left : shell->right, &shell->dshift);
1213: VecSet(shell->dshift,a);
1214: } else {
1215: if (shell->left) {VecPointwiseMult(shell->dshift,shell->dshift,shell->left);}
1216: if (shell->right) {VecPointwiseMult(shell->dshift,shell->dshift,shell->right);}
1217: VecShift(shell->dshift,a);
1218: }
1219: if (shell->left) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->left);}
1220: if (shell->right) {VecPointwiseDivide(shell->dshift,shell->dshift,shell->right);}
1221: } else shell->vshift += a;
1222: if (shell->zrows) {
1223: VecShift(shell->zvals,a);
1224: }
1225: return(0);
1226: }
1228: PetscErrorCode MatDiagonalSet_Shell_Private(Mat A,Vec D,PetscScalar s)
1229: {
1230: Mat_Shell *shell = (Mat_Shell*)A->data;
1234: if (!shell->dshift) {VecDuplicate(D,&shell->dshift);}
1235: if (shell->left || shell->right) {
1236: if (!shell->right_work) {VecDuplicate(shell->left ? shell->left : shell->right, &shell->right_work);}
1237: if (shell->left && shell->right) {
1238: VecPointwiseDivide(shell->right_work,D,shell->left);
1239: VecPointwiseDivide(shell->right_work,shell->right_work,shell->right);
1240: } else if (shell->left) {
1241: VecPointwiseDivide(shell->right_work,D,shell->left);
1242: } else {
1243: VecPointwiseDivide(shell->right_work,D,shell->right);
1244: }
1245: VecAXPY(shell->dshift,s,shell->right_work);
1246: } else {
1247: VecAXPY(shell->dshift,s,D);
1248: }
1249: return(0);
1250: }
1252: PetscErrorCode MatDiagonalSet_Shell(Mat A,Vec D,InsertMode ins)
1253: {
1254: Mat_Shell *shell = (Mat_Shell*)A->data;
1255: Vec d;
1257: PetscBool flg;
1260: MatHasCongruentLayouts(A,&flg);
1261: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Cannot diagonal set or shift shell matrix if it is not congruent");
1262: if (ins == INSERT_VALUES) {
1263: if (!A->ops->getdiagonal) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_SUP,"Operation MATOP_GETDIAGONAL must be set first");
1264: VecDuplicate(D,&d);
1265: MatGetDiagonal(A,d);
1266: MatDiagonalSet_Shell_Private(A,d,-1.);
1267: MatDiagonalSet_Shell_Private(A,D,1.);
1268: VecDestroy(&d);
1269: if (shell->zrows) {
1270: VecCopy(D,shell->zvals);
1271: }
1272: } else {
1273: MatDiagonalSet_Shell_Private(A,D,1.);
1274: if (shell->zrows) {
1275: VecAXPY(shell->zvals,1.0,D);
1276: }
1277: }
1278: return(0);
1279: }
1281: PetscErrorCode MatScale_Shell(Mat Y,PetscScalar a)
1282: {
1283: Mat_Shell *shell = (Mat_Shell*)Y->data;
1287: shell->vscale *= a;
1288: shell->vshift *= a;
1289: if (shell->dshift) {
1290: VecScale(shell->dshift,a);
1291: }
1292: shell->axpy_vscale *= a;
1293: if (shell->zrows) {
1294: VecScale(shell->zvals,a);
1295: }
1296: return(0);
1297: }
1299: static PetscErrorCode MatDiagonalScale_Shell(Mat Y,Vec left,Vec right)
1300: {
1301: Mat_Shell *shell = (Mat_Shell*)Y->data;
1305: if (left) {
1306: if (!shell->left) {
1307: VecDuplicate(left,&shell->left);
1308: VecCopy(left,shell->left);
1309: } else {
1310: VecPointwiseMult(shell->left,shell->left,left);
1311: }
1312: if (shell->zrows) {
1313: VecPointwiseMult(shell->zvals,shell->zvals,left);
1314: }
1315: }
1316: if (right) {
1317: if (!shell->right) {
1318: VecDuplicate(right,&shell->right);
1319: VecCopy(right,shell->right);
1320: } else {
1321: VecPointwiseMult(shell->right,shell->right,right);
1322: }
1323: if (shell->zrows) {
1324: if (!shell->left_work) {
1325: MatCreateVecs(Y,NULL,&shell->left_work);
1326: }
1327: VecSet(shell->zvals_w,1.0);
1328: VecScatterBegin(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
1329: VecScatterEnd(shell->zvals_sct_c,right,shell->zvals_w,INSERT_VALUES,SCATTER_FORWARD);
1330: VecPointwiseMult(shell->zvals,shell->zvals,shell->zvals_w);
1331: }
1332: }
1333: if (shell->axpy) {
1334: MatDiagonalScale(shell->axpy,left,right);
1335: }
1336: return(0);
1337: }
1339: PetscErrorCode MatAssemblyEnd_Shell(Mat Y,MatAssemblyType t)
1340: {
1341: Mat_Shell *shell = (Mat_Shell*)Y->data;
1345: if (t == MAT_FINAL_ASSEMBLY) {
1346: shell->vshift = 0.0;
1347: shell->vscale = 1.0;
1348: shell->axpy_vscale = 0.0;
1349: shell->axpy_state = 0;
1350: VecDestroy(&shell->dshift);
1351: VecDestroy(&shell->left);
1352: VecDestroy(&shell->right);
1353: MatDestroy(&shell->axpy);
1354: VecDestroy(&shell->axpy_left);
1355: VecDestroy(&shell->axpy_right);
1356: VecScatterDestroy(&shell->zvals_sct_c);
1357: VecScatterDestroy(&shell->zvals_sct_r);
1358: ISDestroy(&shell->zrows);
1359: ISDestroy(&shell->zcols);
1360: }
1361: return(0);
1362: }
1364: static PetscErrorCode MatMissingDiagonal_Shell(Mat A,PetscBool *missing,PetscInt *d)
1365: {
1367: *missing = PETSC_FALSE;
1368: return(0);
1369: }
1371: PetscErrorCode MatAXPY_Shell(Mat Y,PetscScalar a,Mat X,MatStructure str)
1372: {
1373: Mat_Shell *shell = (Mat_Shell*)Y->data;
1377: if (X == Y) {
1378: MatScale(Y,1.0 + a);
1379: return(0);
1380: }
1381: if (!shell->axpy) {
1382: MatConvertFrom_Shell(X,MATSHELL,MAT_INITIAL_MATRIX,&shell->axpy);
1383: shell->axpy_vscale = a;
1384: PetscObjectStateGet((PetscObject)X,&shell->axpy_state);
1385: } else {
1386: MatAXPY(shell->axpy,a/shell->axpy_vscale,X,str);
1387: }
1388: return(0);
1389: }
1391: static struct _MatOps MatOps_Values = {NULL,
1392: NULL,
1393: NULL,
1394: NULL,
1395: /* 4*/ MatMultAdd_Shell,
1396: NULL,
1397: MatMultTransposeAdd_Shell,
1398: NULL,
1399: NULL,
1400: NULL,
1401: /*10*/ NULL,
1402: NULL,
1403: NULL,
1404: NULL,
1405: NULL,
1406: /*15*/ NULL,
1407: NULL,
1408: NULL,
1409: MatDiagonalScale_Shell,
1410: NULL,
1411: /*20*/ NULL,
1412: MatAssemblyEnd_Shell,
1413: NULL,
1414: NULL,
1415: /*24*/ MatZeroRows_Shell,
1416: NULL,
1417: NULL,
1418: NULL,
1419: NULL,
1420: /*29*/ NULL,
1421: NULL,
1422: NULL,
1423: NULL,
1424: NULL,
1425: /*34*/ MatDuplicate_Shell,
1426: NULL,
1427: NULL,
1428: NULL,
1429: NULL,
1430: /*39*/ MatAXPY_Shell,
1431: NULL,
1432: NULL,
1433: NULL,
1434: MatCopy_Shell,
1435: /*44*/ NULL,
1436: MatScale_Shell,
1437: MatShift_Shell,
1438: MatDiagonalSet_Shell,
1439: MatZeroRowsColumns_Shell,
1440: /*49*/ NULL,
1441: NULL,
1442: NULL,
1443: NULL,
1444: NULL,
1445: /*54*/ NULL,
1446: NULL,
1447: NULL,
1448: NULL,
1449: NULL,
1450: /*59*/ NULL,
1451: MatDestroy_Shell,
1452: NULL,
1453: MatConvertFrom_Shell,
1454: NULL,
1455: /*64*/ NULL,
1456: NULL,
1457: NULL,
1458: NULL,
1459: NULL,
1460: /*69*/ NULL,
1461: NULL,
1462: MatConvert_Shell,
1463: NULL,
1464: NULL,
1465: /*74*/ NULL,
1466: NULL,
1467: NULL,
1468: NULL,
1469: NULL,
1470: /*79*/ NULL,
1471: NULL,
1472: NULL,
1473: NULL,
1474: NULL,
1475: /*84*/ NULL,
1476: NULL,
1477: NULL,
1478: NULL,
1479: NULL,
1480: /*89*/ NULL,
1481: NULL,
1482: NULL,
1483: NULL,
1484: NULL,
1485: /*94*/ NULL,
1486: NULL,
1487: NULL,
1488: NULL,
1489: NULL,
1490: /*99*/ NULL,
1491: NULL,
1492: NULL,
1493: NULL,
1494: NULL,
1495: /*104*/ NULL,
1496: NULL,
1497: NULL,
1498: NULL,
1499: NULL,
1500: /*109*/ NULL,
1501: NULL,
1502: NULL,
1503: NULL,
1504: MatMissingDiagonal_Shell,
1505: /*114*/ NULL,
1506: NULL,
1507: NULL,
1508: NULL,
1509: NULL,
1510: /*119*/ NULL,
1511: NULL,
1512: NULL,
1513: NULL,
1514: NULL,
1515: /*124*/ NULL,
1516: NULL,
1517: NULL,
1518: NULL,
1519: NULL,
1520: /*129*/ NULL,
1521: NULL,
1522: NULL,
1523: NULL,
1524: NULL,
1525: /*134*/ NULL,
1526: NULL,
1527: NULL,
1528: NULL,
1529: NULL,
1530: /*139*/ NULL,
1531: NULL,
1532: NULL
1533: };
1535: PetscErrorCode MatShellSetContext_Shell(Mat mat,void *ctx)
1536: {
1537: Mat_Shell *shell = (Mat_Shell*)mat->data;
1540: shell->ctx = ctx;
1541: return(0);
1542: }
1544: static PetscErrorCode MatShellSetVecType_Shell(Mat mat,VecType vtype)
1545: {
1549: PetscFree(mat->defaultvectype);
1550: PetscStrallocpy(vtype,(char**)&mat->defaultvectype);
1551: return(0);
1552: }
1554: PetscErrorCode MatShellSetManageScalingShifts_Shell(Mat A)
1555: {
1556: Mat_Shell *shell = (Mat_Shell*)A->data;
1559: shell->managescalingshifts = PETSC_FALSE;
1560: A->ops->diagonalset = NULL;
1561: A->ops->diagonalscale = NULL;
1562: A->ops->scale = NULL;
1563: A->ops->shift = NULL;
1564: A->ops->axpy = NULL;
1565: return(0);
1566: }
1568: PetscErrorCode MatShellSetOperation_Shell(Mat mat,MatOperation op,void (*f)(void))
1569: {
1570: Mat_Shell *shell = (Mat_Shell*)mat->data;
1573: switch (op) {
1574: case MATOP_DESTROY:
1575: shell->ops->destroy = (PetscErrorCode (*)(Mat))f;
1576: break;
1577: case MATOP_VIEW:
1578: if (!mat->ops->viewnative) {
1579: mat->ops->viewnative = mat->ops->view;
1580: }
1581: mat->ops->view = (PetscErrorCode (*)(Mat,PetscViewer))f;
1582: break;
1583: case MATOP_COPY:
1584: shell->ops->copy = (PetscErrorCode (*)(Mat,Mat,MatStructure))f;
1585: break;
1586: case MATOP_DIAGONAL_SET:
1587: case MATOP_DIAGONAL_SCALE:
1588: case MATOP_SHIFT:
1589: case MATOP_SCALE:
1590: case MATOP_AXPY:
1591: case MATOP_ZERO_ROWS:
1592: case MATOP_ZERO_ROWS_COLUMNS:
1593: if (shell->managescalingshifts) SETERRQ(PetscObjectComm((PetscObject)mat),PETSC_ERR_ARG_WRONGSTATE,"MATSHELL is managing scalings and shifts, see MatShellSetManageScalingShifts()");
1594: (((void(**)(void))mat->ops)[op]) = f;
1595: break;
1596: case MATOP_GET_DIAGONAL:
1597: if (shell->managescalingshifts) {
1598: shell->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1599: mat->ops->getdiagonal = MatGetDiagonal_Shell;
1600: } else {
1601: shell->ops->getdiagonal = NULL;
1602: mat->ops->getdiagonal = (PetscErrorCode (*)(Mat,Vec))f;
1603: }
1604: break;
1605: case MATOP_MULT:
1606: if (shell->managescalingshifts) {
1607: shell->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1608: mat->ops->mult = MatMult_Shell;
1609: } else {
1610: shell->ops->mult = NULL;
1611: mat->ops->mult = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1612: }
1613: break;
1614: case MATOP_MULT_TRANSPOSE:
1615: if (shell->managescalingshifts) {
1616: shell->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1617: mat->ops->multtranspose = MatMultTranspose_Shell;
1618: } else {
1619: shell->ops->multtranspose = NULL;
1620: mat->ops->multtranspose = (PetscErrorCode (*)(Mat,Vec,Vec))f;
1621: }
1622: break;
1623: default:
1624: (((void(**)(void))mat->ops)[op]) = f;
1625: break;
1626: }
1627: return(0);
1628: }
1630: PetscErrorCode MatShellGetOperation_Shell(Mat mat,MatOperation op,void(**f)(void))
1631: {
1632: Mat_Shell *shell = (Mat_Shell*)mat->data;
1635: switch (op) {
1636: case MATOP_DESTROY:
1637: *f = (void (*)(void))shell->ops->destroy;
1638: break;
1639: case MATOP_VIEW:
1640: *f = (void (*)(void))mat->ops->view;
1641: break;
1642: case MATOP_COPY:
1643: *f = (void (*)(void))shell->ops->copy;
1644: break;
1645: case MATOP_DIAGONAL_SET:
1646: case MATOP_DIAGONAL_SCALE:
1647: case MATOP_SHIFT:
1648: case MATOP_SCALE:
1649: case MATOP_AXPY:
1650: case MATOP_ZERO_ROWS:
1651: case MATOP_ZERO_ROWS_COLUMNS:
1652: *f = (((void (**)(void))mat->ops)[op]);
1653: break;
1654: case MATOP_GET_DIAGONAL:
1655: if (shell->ops->getdiagonal)
1656: *f = (void (*)(void))shell->ops->getdiagonal;
1657: else
1658: *f = (((void (**)(void))mat->ops)[op]);
1659: break;
1660: case MATOP_MULT:
1661: if (shell->ops->mult)
1662: *f = (void (*)(void))shell->ops->mult;
1663: else
1664: *f = (((void (**)(void))mat->ops)[op]);
1665: break;
1666: case MATOP_MULT_TRANSPOSE:
1667: if (shell->ops->multtranspose)
1668: *f = (void (*)(void))shell->ops->multtranspose;
1669: else
1670: *f = (((void (**)(void))mat->ops)[op]);
1671: break;
1672: default:
1673: *f = (((void (**)(void))mat->ops)[op]);
1674: }
1675: return(0);
1676: }
1678: /*MC
1679: MATSHELL - MATSHELL = "shell" - A matrix type to be used to define your own matrix type -- perhaps matrix free.
1681: Level: advanced
1683: .seealso: MatCreateShell()
1684: M*/
1686: PETSC_EXTERN PetscErrorCode MatCreate_Shell(Mat A)
1687: {
1688: Mat_Shell *b;
1692: PetscMemcpy(A->ops,&MatOps_Values,sizeof(struct _MatOps));
1694: PetscNewLog(A,&b);
1695: A->data = (void*)b;
1697: b->ctx = NULL;
1698: b->vshift = 0.0;
1699: b->vscale = 1.0;
1700: b->managescalingshifts = PETSC_TRUE;
1701: A->assembled = PETSC_TRUE;
1702: A->preallocated = PETSC_FALSE;
1704: PetscObjectComposeFunction((PetscObject)A,"MatShellGetContext_C",MatShellGetContext_Shell);
1705: PetscObjectComposeFunction((PetscObject)A,"MatShellSetContext_C",MatShellSetContext_Shell);
1706: PetscObjectComposeFunction((PetscObject)A,"MatShellSetVecType_C",MatShellSetVecType_Shell);
1707: PetscObjectComposeFunction((PetscObject)A,"MatShellSetManageScalingShifts_C",MatShellSetManageScalingShifts_Shell);
1708: PetscObjectComposeFunction((PetscObject)A,"MatShellSetOperation_C",MatShellSetOperation_Shell);
1709: PetscObjectComposeFunction((PetscObject)A,"MatShellGetOperation_C",MatShellGetOperation_Shell);
1710: PetscObjectComposeFunction((PetscObject)A,"MatShellSetMatProductOperation_C",MatShellSetMatProductOperation_Shell);
1711: PetscObjectChangeTypeName((PetscObject)A,MATSHELL);
1712: return(0);
1713: }
1715: /*@C
1716: MatCreateShell - Creates a new matrix class for use with a user-defined
1717: private data storage format.
1719: Collective
1721: Input Parameters:
1722: + comm - MPI communicator
1723: . m - number of local rows (must be given)
1724: . n - number of local columns (must be given)
1725: . M - number of global rows (may be PETSC_DETERMINE)
1726: . N - number of global columns (may be PETSC_DETERMINE)
1727: - ctx - pointer to data needed by the shell matrix routines
1729: Output Parameter:
1730: . A - the matrix
1732: Level: advanced
1734: Usage:
1735: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1736: $ MatCreateShell(comm,m,n,M,N,ctx,&mat);
1737: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1738: $ [ Use matrix for operations that have been set ]
1739: $ MatDestroy(mat);
1741: Notes:
1742: The shell matrix type is intended to provide a simple class to use
1743: with KSP (such as, for use with matrix-free methods). You should not
1744: use the shell type if you plan to define a complete matrix class.
1746: Fortran Notes:
1747: To use this from Fortran with a ctx you must write an interface definition for this
1748: function and for MatShellGetContext() that tells Fortran the Fortran derived data type you are passing
1749: in as the ctx argument.
1751: PETSc requires that matrices and vectors being used for certain
1752: operations are partitioned accordingly. For example, when
1753: creating a shell matrix, A, that supports parallel matrix-vector
1754: products using MatMult(A,x,y) the user should set the number
1755: of local matrix rows to be the number of local elements of the
1756: corresponding result vector, y. Note that this is information is
1757: required for use of the matrix interface routines, even though
1758: the shell matrix may not actually be physically partitioned.
1759: For example,
1761: $
1762: $ Vec x, y
1763: $ extern PetscErrorCode mult(Mat,Vec,Vec);
1764: $ Mat A
1765: $
1766: $ VecCreateMPI(comm,PETSC_DECIDE,M,&y);
1767: $ VecCreateMPI(comm,PETSC_DECIDE,N,&x);
1768: $ VecGetLocalSize(y,&m);
1769: $ VecGetLocalSize(x,&n);
1770: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
1771: $ MatShellSetOperation(mat,MATOP_MULT,(void(*)(void))mult);
1772: $ MatMult(A,x,y);
1773: $ MatDestroy(&A);
1774: $ VecDestroy(&y);
1775: $ VecDestroy(&x);
1776: $
1779: MATSHELL handles MatShift(), MatDiagonalSet(), MatDiagonalScale(), MatAXPY(), MatScale(), MatZeroRows() and MatZeroRowsColumns() internally, so these
1780: operations cannot be overwritten unless MatShellSetManageScalingShifts() is called.
1783: For rectangular matrices do all the scalings and shifts make sense?
1785: Developers Notes:
1786: Regarding shifting and scaling. The general form is
1788: diag(left)(vscale*A + diag(dshift) + vshift I)diag(right)
1790: The order you apply the operations is important. For example if you have a dshift then
1791: apply a MatScale(s) you get s*vscale*A + s*diag(shift). But if you first scale and then shift
1792: you get s*vscale*A + diag(shift)
1794: A is the user provided function.
1796: KSP/PC uses changes in the Mat's "state" to decide if preconditioners need to be rebuilt: PCSetUp() only calls the setup() for
1797: for the PC implementation if the Mat state has increased from the previous call. Thus to get changes in a MATSHELL to trigger
1798: an update in the preconditioner you must call MatAssemblyBegin()/MatAssemblyEnd() or PetscObjectStateIncrease((PetscObject)mat);
1799: each time the MATSHELL matrix has changed.
1801: Matrix product operations (i.e. MatMat, MatTranposeMat etc) can be specified using MatShellSetMatProductOperation()
1803: Calling MatAssemblyBegin()/MatAssemblyEnd() on a MATSHELL removes any previously supplied shift and scales that were provided
1804: with MatDiagonalSet(), MatShift(), MatScale(), or MatDiagonalScale().
1806: .seealso: MatShellSetOperation(), MatHasOperation(), MatShellGetContext(), MatShellSetContext(), MATSHELL, MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
1807: @*/
1808: PetscErrorCode MatCreateShell(MPI_Comm comm,PetscInt m,PetscInt n,PetscInt M,PetscInt N,void *ctx,Mat *A)
1809: {
1813: MatCreate(comm,A);
1814: MatSetSizes(*A,m,n,M,N);
1815: MatSetType(*A,MATSHELL);
1816: MatShellSetContext(*A,ctx);
1817: MatSetUp(*A);
1818: return(0);
1819: }
1821: /*@
1822: MatShellSetContext - sets the context for a shell matrix
1824: Logically Collective on Mat
1826: Input Parameters:
1827: + mat - the shell matrix
1828: - ctx - the context
1830: Level: advanced
1832: Fortran Notes:
1833: To use this from Fortran you must write a Fortran interface definition for this
1834: function that tells Fortran the Fortran derived data type that you are passing in as the ctx argument.
1836: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation()
1837: @*/
1838: PetscErrorCode MatShellSetContext(Mat mat,void *ctx)
1839: {
1844: PetscTryMethod(mat,"MatShellSetContext_C",(Mat,void*),(mat,ctx));
1845: return(0);
1846: }
1848: /*@C
1849: MatShellSetVecType - Sets the type of Vec returned by MatCreateVecs()
1851: Logically collective
1853: Input Parameters:
1854: + mat - the shell matrix
1855: - vtype - type to use for creating vectors
1857: Notes:
1859: Level: advanced
1861: .seealso: MatCreateVecs()
1862: @*/
1863: PetscErrorCode MatShellSetVecType(Mat mat,VecType vtype)
1864: {
1868: PetscTryMethod(mat,"MatShellSetVecType_C",(Mat,VecType),(mat,vtype));
1869: return(0);
1870: }
1872: /*@
1873: MatShellSetManageScalingShifts - Allows the user to control the scaling and shift operations of the MATSHELL. Must be called immediately
1874: after MatCreateShell()
1876: Logically Collective on Mat
1878: Input Parameter:
1879: . mat - the shell matrix
1881: Level: advanced
1883: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatShellSetOperation()
1884: @*/
1885: PetscErrorCode MatShellSetManageScalingShifts(Mat A)
1886: {
1891: PetscTryMethod(A,"MatShellSetManageScalingShifts_C",(Mat),(A));
1892: return(0);
1893: }
1895: /*@C
1896: MatShellTestMult - Compares the multiply routine provided to the MATSHELL with differencing on a given function.
1898: Logically Collective on Mat
1900: Input Parameters:
1901: + mat - the shell matrix
1902: . f - the function
1903: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1904: - ctx - an optional context for the function
1906: Output Parameter:
1907: . flg - PETSC_TRUE if the multiply is likely correct
1909: Options Database:
1910: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1912: Level: advanced
1914: Fortran Notes:
1915: Not supported from Fortran
1917: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMultTranspose()
1918: @*/
1919: PetscErrorCode MatShellTestMult(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1920: {
1922: PetscInt m,n;
1923: Mat mf,Dmf,Dmat,Ddiff;
1924: PetscReal Diffnorm,Dmfnorm;
1925: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1929: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_view",&v);
1930: MatGetLocalSize(mat,&m,&n);
1931: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,PETSC_DECIDE,PETSC_DECIDE,&mf);
1932: MatMFFDSetFunction(mf,f,ctx);
1933: MatMFFDSetBase(mf,base,NULL);
1935: MatComputeOperator(mf,MATAIJ,&Dmf);
1936: MatComputeOperator(mat,MATAIJ,&Dmat);
1938: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
1939: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
1940: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
1941: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
1942: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
1943: flag = PETSC_FALSE;
1944: if (v) {
1945: 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));
1946: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_view");
1947: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_view");
1948: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_view");
1949: }
1950: } else if (v) {
1951: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL and matrix free multiple appear to produce the same results\n");
1952: }
1953: if (flg) *flg = flag;
1954: MatDestroy(&Ddiff);
1955: MatDestroy(&mf);
1956: MatDestroy(&Dmf);
1957: MatDestroy(&Dmat);
1958: return(0);
1959: }
1961: /*@C
1962: MatShellTestMultTranpose - Compares the multiply transpose routine provided to the MATSHELL with differencing on a given function.
1964: Logically Collective on Mat
1966: Input Parameters:
1967: + mat - the shell matrix
1968: . f - the function
1969: . base - differences are computed around this vector, see MatMFFDSetBase(), for Jacobians this is the point at which the Jacobian is being evaluated
1970: - ctx - an optional context for the function
1972: Output Parameter:
1973: . flg - PETSC_TRUE if the multiply is likely correct
1975: Options Database:
1976: . -mat_shell_test_mult_view - print if any differences are detected between the products and print the difference
1978: Level: advanced
1980: Fortran Notes:
1981: Not supported from Fortran
1983: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellTestMult()
1984: @*/
1985: PetscErrorCode MatShellTestMultTranspose(Mat mat,PetscErrorCode (*f)(void*,Vec,Vec),Vec base,void *ctx,PetscBool *flg)
1986: {
1988: Vec x,y,z;
1989: PetscInt m,n,M,N;
1990: Mat mf,Dmf,Dmat,Ddiff;
1991: PetscReal Diffnorm,Dmfnorm;
1992: PetscBool v = PETSC_FALSE, flag = PETSC_TRUE;
1996: PetscOptionsHasName(NULL,((PetscObject)mat)->prefix,"-mat_shell_test_mult_transpose_view",&v);
1997: MatCreateVecs(mat,&x,&y);
1998: VecDuplicate(y,&z);
1999: MatGetLocalSize(mat,&m,&n);
2000: MatGetSize(mat,&M,&N);
2001: MatCreateMFFD(PetscObjectComm((PetscObject)mat),m,n,M,N,&mf);
2002: MatMFFDSetFunction(mf,f,ctx);
2003: MatMFFDSetBase(mf,base,NULL);
2004: MatComputeOperator(mf,MATAIJ,&Dmf);
2005: MatTranspose(Dmf,MAT_INPLACE_MATRIX,&Dmf);
2006: MatComputeOperatorTranspose(mat,MATAIJ,&Dmat);
2008: MatDuplicate(Dmat,MAT_COPY_VALUES,&Ddiff);
2009: MatAXPY(Ddiff,-1.0,Dmf,DIFFERENT_NONZERO_PATTERN);
2010: MatNorm(Ddiff,NORM_FROBENIUS,&Diffnorm);
2011: MatNorm(Dmf,NORM_FROBENIUS,&Dmfnorm);
2012: if (Diffnorm/Dmfnorm > 10*PETSC_SQRT_MACHINE_EPSILON) {
2013: flag = PETSC_FALSE;
2014: if (v) {
2015: 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));
2016: MatViewFromOptions(Ddiff,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2017: MatViewFromOptions(Dmf,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2018: MatViewFromOptions(Dmat,(PetscObject)mat,"-mat_shell_test_mult_transpose_view");
2019: }
2020: } else if (v) {
2021: PetscPrintf(PetscObjectComm((PetscObject)mat),"MATSHELL transpose and matrix free multiple appear to produce the same results\n");
2022: }
2023: if (flg) *flg = flag;
2024: MatDestroy(&mf);
2025: MatDestroy(&Dmat);
2026: MatDestroy(&Ddiff);
2027: MatDestroy(&Dmf);
2028: VecDestroy(&x);
2029: VecDestroy(&y);
2030: VecDestroy(&z);
2031: return(0);
2032: }
2034: /*@C
2035: MatShellSetOperation - Allows user to set a matrix operation for a shell matrix.
2037: Logically Collective on Mat
2039: Input Parameters:
2040: + mat - the shell matrix
2041: . op - the name of the operation
2042: - g - the function that provides the operation.
2044: Level: advanced
2046: Usage:
2047: $ extern PetscErrorCode usermult(Mat,Vec,Vec);
2048: $ MatCreateShell(comm,m,n,M,N,ctx,&A);
2049: $ MatShellSetOperation(A,MATOP_MULT,(void(*)(void))usermult);
2051: Notes:
2052: See the file include/petscmat.h for a complete list of matrix
2053: operations, which all have the form MATOP_<OPERATION>, where
2054: <OPERATION> is the name (in all capital letters) of the
2055: user interface routine (e.g., MatMult() -> MATOP_MULT).
2057: All user-provided functions (except for MATOP_DESTROY) should have the same calling
2058: sequence as the usual matrix interface routines, since they
2059: are intended to be accessed via the usual matrix interface
2060: routines, e.g.,
2061: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
2063: In particular each function MUST return an error code of 0 on success and
2064: nonzero on failure.
2066: Within each user-defined routine, the user should call
2067: MatShellGetContext() to obtain the user-defined context that was
2068: set by MatCreateShell().
2070: Use MatSetOperation() to set an operation for any matrix type. For matrix product operations (i.e. MatMat, MatTransposeMat etc) use MatShellSetMatProductOperation()
2072: Fortran Notes:
2073: For MatCreateVecs() the user code should check if the input left or right matrix is -1 and in that case not
2074: generate a matrix. See src/mat/tests/ex120f.F
2076: .seealso: MatCreateShell(), MatShellGetContext(), MatShellGetOperation(), MatShellSetContext(), MatSetOperation(), MatShellSetManageScalingShifts(), MatShellSetMatProductOperation()
2077: @*/
2078: PetscErrorCode MatShellSetOperation(Mat mat,MatOperation op,void (*g)(void))
2079: {
2084: PetscTryMethod(mat,"MatShellSetOperation_C",(Mat,MatOperation,void (*)(void)),(mat,op,g));
2085: return(0);
2086: }
2088: /*@C
2089: MatShellGetOperation - Gets a matrix function for a shell matrix.
2091: Not Collective
2093: Input Parameters:
2094: + mat - the shell matrix
2095: - op - the name of the operation
2097: Output Parameter:
2098: . g - the function that provides the operation.
2100: Level: advanced
2102: Notes:
2103: See the file include/petscmat.h for a complete list of matrix
2104: operations, which all have the form MATOP_<OPERATION>, where
2105: <OPERATION> is the name (in all capital letters) of the
2106: user interface routine (e.g., MatMult() -> MATOP_MULT).
2108: All user-provided functions have the same calling
2109: sequence as the usual matrix interface routines, since they
2110: are intended to be accessed via the usual matrix interface
2111: routines, e.g.,
2112: $ MatMult(Mat,Vec,Vec) -> usermult(Mat,Vec,Vec)
2114: Within each user-defined routine, the user should call
2115: MatShellGetContext() to obtain the user-defined context that was
2116: set by MatCreateShell().
2118: .seealso: MatCreateShell(), MatShellGetContext(), MatShellSetOperation(), MatShellSetContext()
2119: @*/
2120: PetscErrorCode MatShellGetOperation(Mat mat,MatOperation op,void(**g)(void))
2121: {
2126: PetscUseMethod(mat,"MatShellGetOperation_C",(Mat,MatOperation,void (**)(void)),(mat,op,g));
2127: return(0);
2128: }
2130: /*@
2131: MatIsShell - Inquires if a matrix is derived from MATSHELL
2133: Input Parameter:
2134: . mat - the matrix
2136: Output Parameter:
2137: . flg - the boolean value
2139: Level: developer
2141: 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)
2143: .seealso: MatCreateShell()
2144: @*/
2145: PetscErrorCode MatIsShell(Mat mat, PetscBool *flg)
2146: {
2150: *flg = (PetscBool)(mat->ops->destroy == MatDestroy_Shell);
2151: return(0);
2152: }