Actual source code: submat.c

petsc-master 2019-12-10
Report Typos and Errors

  2:  #include <petsc/private/matimpl.h>

  4: typedef struct {
  5:   IS          isrow,iscol;      /* rows and columns in submatrix, only used to check consistency */
  6:   Vec         lwork,rwork;      /* work vectors inside the scatters */
  7:   Vec         lwork2,rwork2;    /* work vectors inside the scatters */
  8:   VecScatter  lrestrict,rprolong;
  9:   Mat         A;
 10: } Mat_SubVirtual;

 12: static PetscErrorCode MatScale_SubMatrix(Mat N,PetscScalar a)
 13: {
 14:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

 18:   MatScale(Na->A,a);
 19:   return(0);
 20: }

 22: static PetscErrorCode MatShift_SubMatrix(Mat N,PetscScalar a)
 23: {
 24:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

 28:   MatShift(Na->A,a);
 29:   return(0);
 30: }

 32: static PetscErrorCode MatDiagonalScale_SubMatrix(Mat N,Vec left,Vec right)
 33: {
 34:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

 38:   if (right) {
 39:     VecZeroEntries(Na->rwork);
 40:     VecScatterBegin(Na->rprolong,right,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
 41:     VecScatterEnd(Na->rprolong,right,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
 42:   }
 43:   if (left) {
 44:     VecZeroEntries(Na->lwork);
 45:     VecScatterBegin(Na->lrestrict,left,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
 46:     VecScatterEnd(Na->lrestrict,left,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
 47:   }
 48:   MatDiagonalScale(Na->A,left ? Na->lwork : NULL,right ? Na->rwork : NULL);
 49:   return(0);
 50: }

 52: static PetscErrorCode MatGetDiagonal_SubMatrix(Mat N,Vec d)
 53: {
 54:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

 58:   MatGetDiagonal(Na->A,Na->rwork);
 59:   VecScatterBegin(Na->rprolong,Na->rwork,d,INSERT_VALUES,SCATTER_REVERSE);
 60:   VecScatterEnd(Na->rprolong,Na->rwork,d,INSERT_VALUES,SCATTER_REVERSE);
 61:   return(0);
 62: }

 64: static PetscErrorCode MatMult_SubMatrix(Mat N,Vec x,Vec y)
 65: {
 66:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

 70:   VecZeroEntries(Na->rwork);
 71:   VecScatterBegin(Na->rprolong,x,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
 72:   VecScatterEnd(Na->rprolong,x,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
 73:   MatMult(Na->A,Na->rwork,Na->lwork);
 74:   VecScatterBegin(Na->lrestrict,Na->lwork,y,INSERT_VALUES,SCATTER_FORWARD);
 75:   VecScatterEnd(Na->lrestrict,Na->lwork,y,INSERT_VALUES,SCATTER_FORWARD);
 76:   return(0);
 77: }

 79: static PetscErrorCode MatMultAdd_SubMatrix(Mat N,Vec v1,Vec v2,Vec v3)
 80: {
 81:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

 85:   VecZeroEntries(Na->rwork);
 86:   VecScatterBegin(Na->rprolong,v1,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
 87:   VecScatterEnd(Na->rprolong,v1,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
 88:   if (v1 == v2) {
 89:     MatMultAdd(Na->A,Na->rwork,Na->rwork,Na->lwork);
 90:   } else if (v2 == v3) {
 91:     VecZeroEntries(Na->lwork);
 92:     VecScatterBegin(Na->lrestrict,v2,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
 93:     VecScatterEnd(Na->lrestrict,v2,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
 94:     MatMultAdd(Na->A,Na->rwork,Na->lwork,Na->lwork);
 95:   } else {
 96:     if (!Na->lwork2) {
 97:       VecDuplicate(Na->lwork,&Na->lwork2);
 98:     } else {
 99:       VecZeroEntries(Na->lwork2);
100:     }
101:     VecScatterBegin(Na->lrestrict,v2,Na->lwork2,INSERT_VALUES,SCATTER_REVERSE);
102:     VecScatterEnd(Na->lrestrict,v2,Na->lwork2,INSERT_VALUES,SCATTER_REVERSE);
103:     MatMultAdd(Na->A,Na->rwork,Na->lwork2,Na->lwork);
104:   }
105:   VecScatterBegin(Na->lrestrict,Na->lwork,v3,INSERT_VALUES,SCATTER_FORWARD);
106:   VecScatterEnd(Na->lrestrict,Na->lwork,v3,INSERT_VALUES,SCATTER_FORWARD);
107:   return(0);
108: }

110: static PetscErrorCode MatMultTranspose_SubMatrix(Mat N,Vec x,Vec y)
111: {
112:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

116:   VecZeroEntries(Na->lwork);
117:   VecScatterBegin(Na->lrestrict,x,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
118:   VecScatterEnd(Na->lrestrict,x,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
119:   MatMultTranspose(Na->A,Na->lwork,Na->rwork);
120:   VecScatterBegin(Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);
121:   VecScatterEnd(Na->rprolong,Na->rwork,y,INSERT_VALUES,SCATTER_REVERSE);
122:   return(0);
123: }

125: static PetscErrorCode MatMultTransposeAdd_SubMatrix(Mat N,Vec v1,Vec v2,Vec v3)
126: {
127:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

131:   VecZeroEntries(Na->lwork);
132:   VecScatterBegin(Na->lrestrict,v1,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
133:   VecScatterEnd(Na->lrestrict,v1,Na->lwork,INSERT_VALUES,SCATTER_REVERSE);
134:   if (v1 == v2) {
135:     MatMultTransposeAdd(Na->A,Na->lwork,Na->lwork,Na->rwork);
136:   } else if (v2 == v3) {
137:     VecZeroEntries(Na->rwork);
138:     VecScatterBegin(Na->rprolong,v2,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
139:     VecScatterEnd(Na->rprolong,v2,Na->rwork,INSERT_VALUES,SCATTER_FORWARD);
140:     MatMultTransposeAdd(Na->A,Na->lwork,Na->rwork,Na->rwork);
141:   } else {
142:     if (!Na->rwork2) {
143:       VecDuplicate(Na->rwork,&Na->rwork2);
144:     } else {
145:       VecZeroEntries(Na->rwork2);
146:     }
147:     VecScatterBegin(Na->rprolong,v2,Na->rwork2,INSERT_VALUES,SCATTER_FORWARD);
148:     VecScatterEnd(Na->rprolong,v2,Na->rwork2,INSERT_VALUES,SCATTER_FORWARD);
149:     MatMultTransposeAdd(Na->A,Na->lwork,Na->rwork2,Na->rwork);
150:   }
151:   VecScatterBegin(Na->rprolong,Na->rwork,v3,INSERT_VALUES,SCATTER_REVERSE);
152:   VecScatterEnd(Na->rprolong,Na->rwork,v3,INSERT_VALUES,SCATTER_REVERSE);
153:   return(0);
154: }

156: static PetscErrorCode MatDestroy_SubMatrix(Mat N)
157: {
158:   Mat_SubVirtual *Na = (Mat_SubVirtual*)N->data;

162:   ISDestroy(&Na->isrow);
163:   ISDestroy(&Na->iscol);
164:   VecDestroy(&Na->lwork);
165:   VecDestroy(&Na->rwork);
166:   VecDestroy(&Na->lwork2);
167:   VecDestroy(&Na->rwork2);
168:   VecScatterDestroy(&Na->lrestrict);
169:   VecScatterDestroy(&Na->rprolong);
170:   MatDestroy(&Na->A);
171:   PetscFree(N->data);
172:   return(0);
173: }

175: /*@
176:    MatCreateSubMatrixVirtual - Creates a virtual matrix that acts as a submatrix

178:    Collective on Mat

180:    Input Parameters:
181: +  A - matrix that we will extract a submatrix of
182: .  isrow - rows to be present in the submatrix
183: -  iscol - columns to be present in the submatrix

185:    Output Parameters:
186: .  newmat - new matrix

188:    Level: developer

190:    Notes:
191:    Most will use MatCreateSubMatrix which provides a more efficient representation if it is available.

193: .seealso: MatCreateSubMatrix(), MatSubMatrixVirtualUpdate()
194: @*/
195: PetscErrorCode MatCreateSubMatrixVirtual(Mat A,IS isrow,IS iscol,Mat *newmat)
196: {
197:   Vec            left,right;
198:   PetscInt       m,n;
199:   Mat            N;
200:   Mat_SubVirtual *Na;

208:   *newmat = 0;

210:   MatCreate(PetscObjectComm((PetscObject)A),&N);
211:   ISGetLocalSize(isrow,&m);
212:   ISGetLocalSize(iscol,&n);
213:   MatSetSizes(N,m,n,PETSC_DETERMINE,PETSC_DETERMINE);
214:   PetscObjectChangeTypeName((PetscObject)N,MATSUBMATRIX);

216:   PetscNewLog(N,&Na);
217:   N->data   = (void*)Na;

219:   PetscObjectReference((PetscObject)isrow);
220:   PetscObjectReference((PetscObject)iscol);
221:   Na->isrow = isrow;
222:   Na->iscol = iscol;

224:   PetscFree(N->defaultvectype);
225:   PetscStrallocpy(A->defaultvectype,&N->defaultvectype);
226:   /* Do not use MatConvert directly since MatShell has a duplicate operation which does not increase
227:      the reference count of the context. This is a problem if A is already of type MATSHELL */
228:   MatConvertFrom_Shell(A,MATSHELL,MAT_INITIAL_MATRIX,&Na->A);

230:   N->ops->destroy          = MatDestroy_SubMatrix;
231:   N->ops->mult             = MatMult_SubMatrix;
232:   N->ops->multadd          = MatMultAdd_SubMatrix;
233:   N->ops->multtranspose    = MatMultTranspose_SubMatrix;
234:   N->ops->multtransposeadd = MatMultTransposeAdd_SubMatrix;
235:   N->ops->scale            = MatScale_SubMatrix;
236:   N->ops->diagonalscale    = MatDiagonalScale_SubMatrix;
237:   N->ops->shift            = MatShift_SubMatrix;
238:   N->ops->convert          = MatConvert_Shell;
239:   N->ops->getdiagonal      = MatGetDiagonal_SubMatrix;

241:   MatSetBlockSizesFromMats(N,A,A);
242:   PetscLayoutSetUp(N->rmap);
243:   PetscLayoutSetUp(N->cmap);

245:   MatCreateVecs(A,&Na->rwork,&Na->lwork);
246:   MatCreateVecs(N,&right,&left);
247:   VecScatterCreate(Na->lwork,isrow,left,NULL,&Na->lrestrict);
248:   VecScatterCreate(right,NULL,Na->rwork,iscol,&Na->rprolong);
249:   VecDestroy(&left);
250:   VecDestroy(&right);
251:   MatSetUp(N);

253:   N->assembled = PETSC_TRUE;
254:   *newmat      = N;
255:   return(0);
256: }


259: /*@
260:    MatSubMatrixVirtualUpdate - Updates a submatrix

262:    Collective on Mat

264:    Input Parameters:
265: +  N - submatrix to update
266: .  A - full matrix in the submatrix
267: .  isrow - rows in the update (same as the first time the submatrix was created)
268: -  iscol - columns in the update (same as the first time the submatrix was created)

270:    Level: developer

272:    Notes:
273:    Most will use MatCreateSubMatrix which provides a more efficient representation if it is available.

275: .seealso: MatCreateSubMatrixVirtual()
276: @*/
277: PetscErrorCode  MatSubMatrixVirtualUpdate(Mat N,Mat A,IS isrow,IS iscol)
278: {
280:   PetscBool      flg;
281:   Mat_SubVirtual *Na;

288:   PetscObjectTypeCompare((PetscObject)N,MATSUBMATRIX,&flg);
289:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)A),PETSC_ERR_ARG_WRONG,"Matrix has wrong type");

291:   Na   = (Mat_SubVirtual*)N->data;
292:   ISEqual(isrow,Na->isrow,&flg);
293:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different row indices");
294:   ISEqual(iscol,Na->iscol,&flg);
295:   if (!flg) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_INCOMP,"Cannot update submatrix with different column indices");

297:   PetscFree(N->defaultvectype);
298:   PetscStrallocpy(A->defaultvectype,&N->defaultvectype);
299:   MatDestroy(&Na->A);
300:   /* Do not use MatConvert directly since MatShell has a duplicate operation which does not increase
301:      the reference count of the context. This is a problem if A is already of type MATSHELL */
302:   MatConvertFrom_Shell(A,MATSHELL,MAT_INITIAL_MATRIX,&Na->A);
303:   return(0);
304: }