Actual source code: submatfree.c

petsc-3.11.4 2019-09-28
Report Typos and Errors
  1:  #include <petsctao.h>
  2:  #include <../src/tao/matrix/submatfree.h>

  4: /*@C
  5:   MatCreateSubMatrixFree - Creates a reduced matrix by masking a
  6:   full matrix.

  8:    Collective on matrix

 10:    Input Parameters:
 11: +  mat - matrix of arbitrary type
 12: .  Rows - the rows that will be in the submatrix
 13: -  Cols - the columns that will be in the submatrix

 15:    Output Parameters:
 16: .  J - New matrix

 18:    Notes:
 19:    The caller is responsible for destroying the input objects after matrix J has been destroyed.

 21:    Level: developer

 23: .seealso: MatCreate()
 24: @*/
 25: PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J)
 26: {
 27:   MPI_Comm         comm=PetscObjectComm((PetscObject)mat);
 28:   MatSubMatFreeCtx ctx;
 29:   PetscErrorCode   ierr;
 30:   PetscInt         mloc,nloc,m,n;

 33:   PetscNew(&ctx);
 34:   ctx->A  = mat;
 35:   MatGetSize(mat,&m,&n);
 36:   MatGetLocalSize(mat,&mloc,&nloc);
 37:   MatCreateVecs(mat,NULL,&ctx->VC);
 38:   ctx->VR = ctx->VC;
 39:    PetscObjectReference((PetscObject)mat);


 42:   ctx->Rows = Rows;
 43:   ctx->Cols = Cols;
 44:   PetscObjectReference((PetscObject)Rows);
 45:   PetscObjectReference((PetscObject)Cols);
 46:   MatCreateShell(comm,mloc,nloc,m,n,ctx,J);
 47:   MatShellSetManageScalingShifts(*J);
 48:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);
 49:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);
 50:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);
 51:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);
 52:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);
 53:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);
 54:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);
 55:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);
 56:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);
 57:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);
 58:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_SMF);
 59:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);
 60:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);
 61:   MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_SMF);
 62:   MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);

 64:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
 65:   return(0);
 66: }

 68: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){
 69:   MatSubMatFreeCtx ctx;
 70:   PetscErrorCode   ierr;

 73:   MatShellGetContext(mat,(void **)&ctx);
 74:   ISDestroy(&ctx->Rows);
 75:   ISDestroy(&ctx->Cols);
 76:   PetscObjectReference((PetscObject)Rows);
 77:   PetscObjectReference((PetscObject)Cols);
 78:   ctx->Cols=Cols;
 79:   ctx->Rows=Rows;
 80:   return(0);
 81: }

 83: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
 84: {
 85:   MatSubMatFreeCtx ctx;
 86:   PetscErrorCode   ierr;

 89:   MatShellGetContext(mat,(void **)&ctx);
 90:   VecCopy(a,ctx->VR);
 91:   VecISSet(ctx->VR,ctx->Cols,0.0);
 92:   MatMult(ctx->A,ctx->VR,y);
 93:   VecISSet(y,ctx->Rows,0.0);
 94:   return(0);
 95: }

 97: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
 98: {
 99:   MatSubMatFreeCtx ctx;
100:   PetscErrorCode   ierr;

103:   MatShellGetContext(mat,(void **)&ctx);
104:   VecCopy(a,ctx->VC);
105:   VecISSet(ctx->VC,ctx->Rows,0.0);
106:   MatMultTranspose(ctx->A,ctx->VC,y);
107:   VecISSet(y,ctx->Cols,0.0);
108:   return(0);
109: }

111: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
112: {
113:   MatSubMatFreeCtx ctx;
114:   PetscErrorCode   ierr;

117:   MatShellGetContext(M,(void **)&ctx);
118:   MatDiagonalSet(ctx->A,D,is);
119:   return(0);
120: }

122: PetscErrorCode MatDestroy_SMF(Mat mat)
123: {
124:   PetscErrorCode   ierr;
125:   MatSubMatFreeCtx ctx;

128:   MatShellGetContext(mat,(void **)&ctx);
129:   MatDestroy(&ctx->A);
130:   ISDestroy(&ctx->Rows);
131:   ISDestroy(&ctx->Cols);
132:   VecDestroy(&ctx->VC);
133:   PetscFree(ctx);
134:   return(0);
135: }



139: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
140: {
141:   PetscErrorCode   ierr;
142:   MatSubMatFreeCtx ctx;

145:   MatShellGetContext(mat,(void **)&ctx);
146:   MatView(ctx->A,viewer);
147:   return(0);
148: }

150: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
151: {
152:   PetscErrorCode   ierr;
153:   MatSubMatFreeCtx ctx;

156:   MatShellGetContext(Y,(void **)&ctx);
157:   MatShift(ctx->A,a);
158:   return(0);
159: }

161: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
162: {
163:   PetscErrorCode   ierr;
164:   MatSubMatFreeCtx ctx;

167:   MatShellGetContext(mat,(void **)&ctx);
168:   MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
169:   return(0);
170: }

172: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
173: {
174:   PetscErrorCode    ierr;
175:   MatSubMatFreeCtx  ctx1,ctx2;
176:   PetscBool         flg1,flg2,flg3;

179:   MatShellGetContext(A,(void **)&ctx1);
180:   MatShellGetContext(B,(void **)&ctx2);
181:   ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
182:   ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
183:   if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
184:     *flg=PETSC_FALSE;
185:   } else {
186:     MatEqual(ctx1->A,ctx2->A,&flg1);
187:     if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
188:     else { *flg=PETSC_TRUE;}
189:   }
190:   return(0);
191: }

193: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
194: {
195:   PetscErrorCode   ierr;
196:   MatSubMatFreeCtx ctx;

199:   MatShellGetContext(mat,(void **)&ctx);
200:   MatScale(ctx->A,a);
201:   return(0);
202: }

204: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
205: {
207:   PetscFunctionReturn(1);
208: }

210: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
211: {
212:   PetscErrorCode   ierr;
213:   MatSubMatFreeCtx ctx;

216:   MatShellGetContext(mat,(void **)&ctx);
217:   MatGetDiagonal(ctx->A,v);
218:   return(0);
219: }

221: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
222: {
223:   MatSubMatFreeCtx ctx;
224:   PetscErrorCode   ierr;

227:   MatShellGetContext(M,(void **)&ctx);
228:   MatGetRowMax(ctx->A,D,NULL);
229:   return(0);
230: }

232: PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
233: {
235:   PetscInt       i;

238:   if (scall == MAT_INITIAL_MATRIX) {
239:     PetscCalloc1(n+1,B );
240:   }

242:   for ( i=0; i<n; i++ ) {
243:     MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
244:   }
245:   return(0);
246: }

248: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
249:                         Mat *newmat)
250: {
251:   PetscErrorCode   ierr;
252:   MatSubMatFreeCtx ctx;

255:   MatShellGetContext(mat,(void **)&ctx);
256:   if (newmat){
257:     MatDestroy(&*newmat);
258:   }
259:   MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
260:   return(0);
261: }

263: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
264: {
265:   PetscErrorCode   ierr;
266:   MatSubMatFreeCtx ctx;

269:   MatShellGetContext(mat,(void **)&ctx);
270:   MatGetRow(ctx->A,row,ncols,cols,vals);
271:   return(0);
272: }

274: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
275: {
276:   PetscErrorCode   ierr;
277:   MatSubMatFreeCtx ctx;

280:   MatShellGetContext(mat,(void **)&ctx);
281:   MatRestoreRow(ctx->A,row,ncols,cols,vals);
282:   return(0);
283: }

285: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
286: {
287:   PetscErrorCode   ierr;
288:   MatSubMatFreeCtx ctx;

291:   MatShellGetContext(mat,(void **)&ctx);
292:   MatGetColumnVector(ctx->A,Y,col);
293:   return(0);
294: }

296: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
297: {
298:   PetscErrorCode    ierr;
299:   MatSubMatFreeCtx  ctx;

302:   MatShellGetContext(mat,(void **)&ctx);
303:   if (type == NORM_FROBENIUS) {
304:     *norm = 1.0;
305:   } else if (type == NORM_1 || type == NORM_INFINITY) {
306:     *norm = 1.0;
307:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
308:   return(0);
309: }