Actual source code: submatfree.c

  1: #include "tao_util.h"   /*I "tao_util.h" I*/
  2: #include "submatfree.h" /*I "submatfree.h" I*/



  8: /*@C
  9:   MatCreateSubMatrixFree - Creates a reduced matrix by masking a
 10:   full matrix.

 12:    Collective on matrix

 14:    Input Parameters:
 15: +  mat - matrix of arbitrary type
 16: .  Rows - the rows that will be in the submatrix
 17: -  Cols - the columns that will be in the submatrix

 19:    Output Parameters:
 20: .  J - New matrix

 22:    Notes: 
 23:    The user provides the input data and is responsible for destroying
 24:    this data after matrix J has been destroyed.  
 25:  
 26:    Level: developer

 28: .seealso: MatCreate()
 29: @*/
 30: PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J)
 31: {
 32:   MPI_Comm     comm=((PetscObject)mat)->comm;
 33:   MatSubMatFreeCtx ctx;
 35:   PetscMPIInt size;
 36:   PetscInt mloc,nloc,m,n;

 39:   PetscNew(_p_MatSubMatFreeCtx,&ctx);
 40:   ctx->A=mat;
 41:   MatGetSize(mat,&m,&n);
 42:   MatGetLocalSize(mat,&mloc,&nloc);
 43:   MPI_Comm_size(comm,&size);
 44:   if (size == 1) {
 45:     VecCreateSeq(comm,n,&ctx->VC); 
 46:   } else {
 47:     VecCreateMPI(comm,nloc,n,&ctx->VC);
 48:   }
 49:   ctx->VR=ctx->VC;
 50:    PetscObjectReference((PetscObject)mat);
 51:   
 52:   
 53:   ctx->Rows = Rows;
 54:   ctx->Cols = Cols;
 55:   PetscObjectReference((PetscObject)Rows); 
 56:   PetscObjectReference((PetscObject)Cols); 
 57:   MatCreateShell(comm,mloc,nloc,m,n,ctx,J);

 59:   MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);
 60:   MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);
 61:   MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);
 62:   MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);
 63:   MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);
 64:   MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);
 65:   MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);
 66:   MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);
 67:   MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);
 68:   MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);
 69:   MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)(void))MatGetSubMatrices_SMF);
 70:   MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);
 71:   MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);
 72:   MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)(void))MatGetSubMatrix_SMF);
 73:   MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);

 75:   PetscLogObjectParent(mat,*J); 

 77:   return(0);  
 78: }

 82: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){
 83:   MatSubMatFreeCtx ctx;

 87:   MatShellGetContext(mat,(void **)&ctx);
 88:   ISDestroy(&ctx->Rows);
 89:   ISDestroy(&ctx->Cols);
 90:   PetscObjectReference((PetscObject)Rows);
 91:   PetscObjectReference((PetscObject)Cols);
 92:   ctx->Cols=Cols;
 93:   ctx->Rows=Rows;
 94:   return(0);  
 95: }

 99: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
100: {
101:   MatSubMatFreeCtx ctx;

105:   MatShellGetContext(mat,(void **)&ctx);
106:   VecCopy(a,ctx->VR);
107:   VecISSetToConstant(ctx->Cols,0.0,ctx->VR);
108:   MatMult(ctx->A,ctx->VR,y);
109:   VecISSetToConstant(ctx->Rows,0.0,y);
110:   return(0);
111: }

115: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
116: {
117:   MatSubMatFreeCtx ctx;

121:   MatShellGetContext(mat,(void **)&ctx);
122:   VecCopy(a,ctx->VC);
123:   VecISSetToConstant(ctx->Rows,0.0,ctx->VC);
124:   MatMultTranspose(ctx->A,ctx->VC,y);
125:   VecISSetToConstant(ctx->Cols,0.0,y);
126:   return(0);
127: } 

131: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
132: {
133:   MatSubMatFreeCtx ctx;

137:   MatShellGetContext(M,(void **)&ctx);
138:   MatDiagonalSet(ctx->A,D,is);
139:   return(0);
140: } 

144: PetscErrorCode MatDestroy_SMF(Mat mat)
145: {
147:   MatSubMatFreeCtx ctx;

150:   ierr =MatShellGetContext(mat,(void **)&ctx);
151:   if (ctx->A) {
152:     ierr =MatDestroy(&ctx->A);
153:   }
154:   if (ctx->Rows) {
155:     ierr =ISDestroy(&ctx->Rows);
156:   }
157:   if (ctx->Cols) {
158:     ierr =ISDestroy(&ctx->Cols);
159:   }
160:   if (ctx->VC) {
161:     ierr =VecDestroy(&ctx->VC);
162:   }
163:   PetscFree(ctx); 
164:   return(0);
165: }



171: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
172: {
174:   MatSubMatFreeCtx ctx;

177:   MatShellGetContext(mat,(void **)&ctx);
178:   MatView(ctx->A,viewer);
179:   return(0);
180: }

184: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
185: {
187:   MatSubMatFreeCtx ctx;

190:   MatShellGetContext(Y,(void **)&ctx);
191:   MatShift(ctx->A,a);
192:   return(0);
193: }

197: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
198: {
200:   MatSubMatFreeCtx ctx;

203:   MatShellGetContext(mat,(void **)&ctx);
204:   MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
205:   return(0);
206: }

210: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
211: {
213:   MatSubMatFreeCtx  ctx1,ctx2;
214:   PetscBool flg1,flg2,flg3;

217:   MatShellGetContext(A,(void **)&ctx1);
218:   MatShellGetContext(B,(void **)&ctx2);
219:   ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
220:   ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
221:   if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
222:     *flg=PETSC_FALSE;
223:   } else {
224:     MatEqual(ctx1->A,ctx2->A,&flg1);
225:     if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;} 
226:     else { *flg=PETSC_TRUE;} 
227:   }
228:   return(0);
229: }

233: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
234: {
236:   MatSubMatFreeCtx ctx;

239:   MatShellGetContext(mat,(void **)&ctx);
240:   MatScale(ctx->A,a);
241:   return(0);
242: }

246: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
247: {
249:   PetscFunctionReturn(1);
250: }

254: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
255: {
257:   MatSubMatFreeCtx ctx;

260:   MatShellGetContext(mat,(void **)&ctx);
261:   MatGetDiagonal(ctx->A,v);
262:   return(0);
263: }

267: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
268: {
269:   MatSubMatFreeCtx ctx;

273:   MatShellGetContext(M,(void **)&ctx);
274:   MatGetRowMax(ctx->A,D,PETSC_NULL);
275:   return(0);
276: } 


281: PetscErrorCode MatGetSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
282: {
284:   PetscInt i;

287:   if (scall == MAT_INITIAL_MATRIX) {
288:     PetscMalloc( (n+1)*sizeof(Mat),B );
289:   }

291:   for ( i=0; i<n; i++ ) {
292:     MatGetSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
293:   }
294:   return(0);
295: }

299: PetscErrorCode MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
300:                         Mat *newmat)
301: {
303:   MatSubMatFreeCtx ctx;

306:   MatShellGetContext(mat,(void **)&ctx);
307:   if (newmat){
308:     ierr =MatDestroy(&*newmat);
309:   }
310:   MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
311:   return(0);
312: }

316: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscReal **vals)
317: {
319:   MatSubMatFreeCtx ctx;

322:   MatShellGetContext(mat,(void **)&ctx);
323:   MatGetRow(ctx->A,row,ncols,cols,vals);
324:   return(0);
325: }

329: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscReal **vals)
330: {
332:   MatSubMatFreeCtx ctx;

335:   MatShellGetContext(mat,(void **)&ctx);
336:   MatRestoreRow(ctx->A,row,ncols,cols,vals);
337:   return(0);
338: }

342: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
343: {
345:   MatSubMatFreeCtx  ctx;

348:   MatShellGetContext(mat,(void **)&ctx);
349:   MatGetColumnVector(ctx->A,Y,col);
350:   return(0);
351: }

355: PetscErrorCode MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat)
356: {
358:   PetscMPIInt size;
359:   MatSubMatFreeCtx  ctx;

362:   MatShellGetContext(mat,(void **)&ctx);
363:   MPI_Comm_size(((PetscObject)mat)->comm,&size);
364:   PetscFunctionReturn(1);
365: }

369: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
370: {
372:   MatSubMatFreeCtx  ctx;

375:   MatShellGetContext(mat,(void **)&ctx);

377:   if (type == NORM_FROBENIUS) {
378:     *norm = 1.0;
379:   } else if (type == NORM_1 || type == NORM_INFINITY) {
380:     *norm = 1.0;
381:   } else {
382:     SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
383:   }
384:   return(0);
385: }