Actual source code: submat.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/matimpl.h> /*I "petscmat.h" I*/
4: typedef struct {
5: IS isrow,iscol; /* rows and columns in submatrix, only used to check consistency */
6: Vec left,right; /* optional scaling */
7: Vec olwork,orwork; /* work vectors outside the scatters, only touched by PreScale and only created if needed*/
8: Vec lwork,rwork; /* work vectors inside the scatters */
9: VecScatter lrestrict,rprolong;
10: Mat A;
11: PetscScalar scale;
12: } Mat_SubMatrix;
16: static PetscErrorCode PreScaleLeft(Mat N,Vec x,Vec *xx)
17: {
18: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
22: if (!Na->left) {
23: *xx = x;
24: } else {
25: if (!Na->olwork) {
26: VecDuplicate(Na->left,&Na->olwork);
27: }
28: VecPointwiseMult(Na->olwork,x,Na->left);
29: *xx = Na->olwork;
30: }
31: return(0);
32: }
36: static PetscErrorCode PreScaleRight(Mat N,Vec x,Vec *xx)
37: {
38: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
42: if (!Na->right) {
43: *xx = x;
44: } else {
45: if (!Na->orwork) {
46: VecDuplicate(Na->right,&Na->orwork);
47: }
48: VecPointwiseMult(Na->orwork,x,Na->right);
49: *xx = Na->orwork;
50: }
51: return(0);
52: }
56: static PetscErrorCode PostScaleLeft(Mat N,Vec x)
57: {
58: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
62: if (Na->left) {
63: VecPointwiseMult(x,x,Na->left);
64: }
65: return(0);
66: }
70: static PetscErrorCode PostScaleRight(Mat N,Vec x)
71: {
72: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
76: if (Na->right) {
77: VecPointwiseMult(x,x,Na->right);
78: }
79: return(0);
80: }
84: static PetscErrorCode MatScale_SubMatrix(Mat N,PetscScalar scale)
85: {
86: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
89: Na->scale *= scale;
90: return(0);
91: }
95: static PetscErrorCode MatDiagonalScale_SubMatrix(Mat N,Vec left,Vec right)
96: {
97: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
101: if (left) {
102: if (!Na->left) {
103: VecDuplicate(left,&Na->left);
104: VecCopy(left,Na->left);
105: } else {
106: VecPointwiseMult(Na->left,left,Na->left);
107: }
108: }
109: if (right) {
110: if (!Na->right) {
111: VecDuplicate(right,&Na->right);
112: VecCopy(right,Na->right);
113: } else {
114: VecPointwiseMult(Na->right,right,Na->right);
115: }
116: }
117: return(0);
118: }
122: static PetscErrorCode MatMult_SubMatrix(Mat N,Vec x,Vec y)
123: {
124: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
125: Vec xx = 0;
129: PreScaleRight(N,x,&xx);
130: VecZeroEntries(Na->rwork);
131: VecScatterBegin(Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
132: VecScatterEnd (Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
133: MatMult(Na->A,Na->rwork,Na->lwork);
134: VecScatterBegin(Na->lrestrict,Na->lwork,y,INSERT_VALUES,SCATTER_FORWARD);
135: VecScatterEnd (Na->lrestrict,Na->lwork,y,INSERT_VALUES,SCATTER_FORWARD);
136: PostScaleLeft(N,y);
137: VecScale(y,Na->scale);
138: return(0);
139: }
143: static PetscErrorCode MatMultAdd_SubMatrix(Mat N,Vec v1,Vec v2,Vec v3)
144: {
145: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
146: Vec xx = 0;
150: PreScaleRight(N,v1,&xx);
151: VecZeroEntries(Na->rwork);
152: VecScatterBegin(Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
153: VecScatterEnd (Na->rprolong,xx,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
154: MatMult(Na->A,Na->rwork,Na->lwork);
155: if (v2 == v3) {
156: if (Na->scale == (PetscScalar)1.0 && !Na->left) {
157: VecScatterBegin(Na->lrestrict,Na->lwork,v3,ADD_VALUES,SCATTER_FORWARD);
158: VecScatterEnd (Na->lrestrict,Na->lwork,v3,ADD_VALUES,SCATTER_FORWARD);
159: } else {
160: if (!Na->olwork) {VecDuplicate(v3,&Na->olwork);}
161: VecScatterBegin(Na->lrestrict,Na->lwork,Na->olwork,INSERT_VALUES,SCATTER_FORWARD);
162: VecScatterEnd (Na->lrestrict,Na->lwork,Na->olwork,INSERT_VALUES,SCATTER_FORWARD);
163: PostScaleLeft(N,Na->olwork);
164: VecAXPY(v3,Na->scale,Na->olwork);
165: }
166: } else {
167: VecScatterBegin(Na->lrestrict,Na->lwork,v3,INSERT_VALUES,SCATTER_FORWARD);
168: VecScatterEnd (Na->lrestrict,Na->lwork,v3,INSERT_VALUES,SCATTER_FORWARD);
169: PostScaleLeft(N,v3);
170: VecAYPX(v3,Na->scale,v2);
171: }
172: return(0);
173: }
177: static PetscErrorCode MatMultTranspose_SubMatrix(Mat N,Vec x,Vec y)
178: {
179: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
180: Vec xx = 0;
184: PreScaleLeft(N,x,&xx);
185: VecZeroEntries(Na->lwork);
186: VecScatterBegin(Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
187: VecScatterEnd (Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
188: MatMultTranspose(Na->A,Na->lwork,Na->rwork);
189: VecScatterBegin(Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);
190: VecScatterEnd (Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);
191: PostScaleRight(N,y);
192: VecScale(y,Na->scale);
193: return(0);
194: }
198: static PetscErrorCode MatMultTransposeAdd_SubMatrix(Mat N,Vec v1,Vec v2,Vec v3)
199: {
200: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
201: Vec xx = 0;
205: PreScaleLeft(N,v1,&xx);
206: VecZeroEntries(Na->lwork);
207: VecScatterBegin(Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
208: VecScatterEnd (Na->lrestrict,xx,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
209: MatMultTranspose(Na->A,Na->lwork,Na->rwork);
210: if (v2 == v3) {
211: if (Na->scale == (PetscScalar)1.0 && !Na->right) {
212: VecScatterBegin(Na->rprolong,Na->rwork,v3,ADD_VALUES,SCATTER_REVERSE);
213: VecScatterEnd (Na->rprolong,Na->rwork,v3,ADD_VALUES,SCATTER_REVERSE);
214: } else {
215: if (!Na->orwork) {VecDuplicate(v3,&Na->orwork);}
216: VecScatterBegin(Na->rprolong,Na->rwork,Na->orwork,INSERT_VALUES,SCATTER_REVERSE);
217: VecScatterEnd (Na->rprolong,Na->rwork,Na->orwork,INSERT_VALUES,SCATTER_REVERSE);
218: PostScaleRight(N,Na->orwork);
219: VecAXPY(v3,Na->scale,Na->orwork);
220: }
221: } else {
222: VecScatterBegin(Na->rprolong,Na->rwork,v3,INSERT_VALUES,SCATTER_REVERSE);
223: VecScatterEnd (Na->rprolong,Na->rwork,v3,INSERT_VALUES,SCATTER_REVERSE);
224: PostScaleRight(N,v3);
225: VecAYPX(v3,Na->scale,v2);
226: }
227: return(0);
228: }
232: static PetscErrorCode MatDestroy_SubMatrix(Mat N)
233: {
234: Mat_SubMatrix *Na = (Mat_SubMatrix*)N->data;
238: ISDestroy(&Na->isrow);
239: ISDestroy(&Na->iscol);
240: VecDestroy(&Na->left);
241: VecDestroy(&Na->right);
242: VecDestroy(&Na->olwork);
243: VecDestroy(&Na->orwork);
244: VecDestroy(&Na->lwork);
245: VecDestroy(&Na->rwork);
246: VecScatterDestroy(&Na->lrestrict);
247: VecScatterDestroy(&Na->rprolong);
248: MatDestroy(&Na->A);
249: PetscFree(N->data);
250: return(0);
251: }
255: /*@
256: MatCreateSubMatrix - Creates a composite matrix that acts as a submatrix
258: Collective on Mat
260: Input Parameters:
261: + A - matrix that we will extract a submatrix of
262: . isrow - rows to be present in the submatrix
263: - iscol - columns to be present in the submatrix
265: Output Parameters:
266: . newmat - new matrix
268: Level: developer
270: Notes:
271: Most will use MatGetSubMatrix which provides a more efficient representation if it is available.
273: .seealso: MatGetSubMatrix(), MatSubMatrixUpdate()
274: @*/
275: PetscErrorCode MatCreateSubMatrix(Mat A,IS isrow,IS iscol,Mat *newmat)
276: {
277: Vec left,right;
278: PetscInt m,n;
279: Mat N;
280: Mat_SubMatrix *Na;
288: *newmat = 0;
290: MatCreate(PetscObjectComm((PetscObject)A),&N);
291: ISGetLocalSize(isrow,&m);
292: ISGetLocalSize(iscol,&n);
293: MatSetSizes(N,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
294: PetscObjectChangeTypeName((PetscObject)N,MATSUBMATRIX);
296: PetscNewLog(N,&Na);
297: N->data = (void*)Na;
298: PetscObjectReference((PetscObject)A);
299: PetscObjectReference((PetscObject)isrow);
300: PetscObjectReference((PetscObject)iscol);
301: Na->A = A;
302: Na->isrow = isrow;
303: Na->iscol = iscol;
304: Na->scale = 1.0;
306: N->ops->destroy = MatDestroy_SubMatrix;
307: N->ops->mult = MatMult_SubMatrix;
308: N->ops->multadd = MatMultAdd_SubMatrix;
309: N->ops->multtranspose = MatMultTranspose_SubMatrix;
310: N->ops->multtransposeadd = MatMultTransposeAdd_SubMatrix;
311: N->ops->scale = MatScale_SubMatrix;
312: N->ops->diagonalscale = MatDiagonalScale_SubMatrix;
314: MatSetBlockSizesFromMats(N,A,A);
315: PetscLayoutSetUp(N->rmap);
316: PetscLayoutSetUp(N->cmap);
318: MatCreateVecs(A,&Na->rwork,&Na->lwork);
319: VecCreate(PetscObjectComm((PetscObject)isrow),&left);
320: VecCreate(PetscObjectComm((PetscObject)iscol),&right);
321: VecSetSizes(left,m,PETSC_DETERMINE);
322: VecSetSizes(right,n,PETSC_DETERMINE);
323: VecSetUp(left);
324: VecSetUp(right);
325: VecScatterCreate(Na->lwork,isrow,left,NULL,&Na->lrestrict);
326: VecScatterCreate(right,NULL,Na->rwork,iscol,&Na->rprolong);
327: VecDestroy(&left);
328: VecDestroy(&right);
330: N->assembled = PETSC_TRUE;
332: MatSetUp(N);
334: *newmat = N;
335: return(0);
336: }
341: /*@
342: MatSubMatrixUpdate - Updates a submatrix
344: Collective on Mat
346: Input Parameters:
347: + N - submatrix to update
348: . A - full matrix in the submatrix
349: . isrow - rows in the update (same as the first time the submatrix was created)
350: - iscol - columns in the update (same as the first time the submatrix was created)
352: Level: developer
354: Notes:
355: Most will use MatGetSubMatrix which provides a more efficient representation if it is available.
357: .seealso: MatGetSubMatrix(), MatCreateSubMatrix()
358: @*/
359: PetscErrorCode MatSubMatrixUpdate(Mat N,Mat A,IS isrow,IS iscol)
360: {
362: PetscBool flg;
363: Mat_SubMatrix *Na;
370: PetscObjectTypeCompare((PetscObject)N,MATSUBMATRIX,&flg);
371: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix has wrong type");
373: Na = (Mat_SubMatrix*)N->data;
374: ISEqual(isrow,Na->isrow,&flg);
375: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different row indices");
376: ISEqual(iscol,Na->iscol,&flg);
377: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different column indices");
379: PetscObjectReference((PetscObject)A);
380: MatDestroy(&Na->A);
381: Na->A = A;
383: Na->scale = 1.0;
384: VecDestroy(&Na->left);
385: VecDestroy(&Na->right);
386: return(0);
387: }