Actual source code: submat.c
petsc-3.4.5 2014-06-29
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,Mat_SubMatrix,&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: PetscLayoutSetBlockSize(N->rmap,A->rmap->bs);
315: PetscLayoutSetBlockSize(N->cmap,A->cmap->bs);
316: PetscLayoutSetUp(N->rmap);
317: PetscLayoutSetUp(N->cmap);
319: MatGetVecs(A,&Na->rwork,&Na->lwork);
320: VecCreate(PetscObjectComm((PetscObject)isrow),&left);
321: VecCreate(PetscObjectComm((PetscObject)iscol),&right);
322: VecSetSizes(left,m,PETSC_DETERMINE);
323: VecSetSizes(right,n,PETSC_DETERMINE);
324: VecSetUp(left);
325: VecSetUp(right);
326: VecScatterCreate(Na->lwork,isrow,left,NULL,&Na->lrestrict);
327: VecScatterCreate(right,NULL,Na->rwork,iscol,&Na->rprolong);
328: VecDestroy(&left);
329: VecDestroy(&right);
331: N->assembled = PETSC_TRUE;
333: MatSetUp(N);
335: *newmat = N;
336: return(0);
337: }
342: /*@
343: MatSubMatrixUpdate - Updates a submatrix
345: Collective on Mat
347: Input Parameters:
348: + N - submatrix to update
349: . A - full matrix in the submatrix
350: . isrow - rows in the update (same as the first time the submatrix was created)
351: - iscol - columns in the update (same as the first time the submatrix was created)
353: Level: developer
355: Notes:
356: Most will use MatGetSubMatrix which provides a more efficient representation if it is available.
358: .seealso: MatGetSubMatrix(), MatCreateSubMatrix()
359: @*/
360: PetscErrorCode MatSubMatrixUpdate(Mat N,Mat A,IS isrow,IS iscol)
361: {
363: PetscBool flg;
364: Mat_SubMatrix *Na;
371: PetscObjectTypeCompare((PetscObject)N,MATSUBMATRIX,&flg);
372: if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix has wrong type");
374: Na = (Mat_SubMatrix*)N->data;
375: ISEqual(isrow,Na->isrow,&flg);
376: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different row indices");
377: ISEqual(iscol,Na->iscol,&flg);
378: if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different column indices");
380: PetscObjectReference((PetscObject)A);
381: MatDestroy(&Na->A);
382: Na->A = A;
384: Na->scale = 1.0;
385: VecDestroy(&Na->left);
386: VecDestroy(&Na->right);
387: return(0);
388: }