Actual source code: submatfree.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: #include <petsctao.h>   /*I "petsctao.h" I*/
  2: #include <../src/tao/matrix/submatfree.h> /*I "submatfree.h" I*/

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

 10:    Collective on matrix

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

 17:    Output Parameters:
 18: .  J - New matrix

 20:    Notes:
 21:    The user provides the input data and is responsible for destroying
 22:    this data after matrix J has been destroyed.

 24:    Level: developer

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

 37:   PetscNew(&ctx);
 38:   ctx->A=mat;
 39:   MatGetSize(mat,&m,&n);
 40:   MatGetLocalSize(mat,&mloc,&nloc);
 41:   MPI_Comm_size(comm,&size);
 42:   if (size == 1) {
 43:     VecCreateSeq(comm,n,&ctx->VC);
 44:   } else {
 45:     VecCreateMPI(comm,nloc,n,&ctx->VC);
 46:   }
 47:   ctx->VR=ctx->VC;
 48:    PetscObjectReference((PetscObject)mat);


 51:   ctx->Rows = Rows;
 52:   ctx->Cols = Cols;
 53:   PetscObjectReference((PetscObject)Rows);
 54:   PetscObjectReference((PetscObject)Cols);
 55:   MatCreateShell(comm,mloc,nloc,m,n,ctx,J);

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

 73:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
 74:   return(0);
 75: }

 79: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){
 80:   MatSubMatFreeCtx ctx;
 81:   PetscErrorCode   ierr;

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

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

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

112: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
113: {
114:   MatSubMatFreeCtx ctx;
115:   PetscErrorCode   ierr;

118:   MatShellGetContext(mat,(void **)&ctx);
119:   VecCopy(a,ctx->VC);
120:   VecISSet(ctx->VC,ctx->Rows,0.0);
121:   MatMultTranspose(ctx->A,ctx->VC,y);
122:   VecISSet(y,ctx->Cols,0.0);
123:   return(0);
124: }

128: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
129: {
130:   MatSubMatFreeCtx ctx;
131:   PetscErrorCode   ierr;

134:   MatShellGetContext(M,(void **)&ctx);
135:   MatDiagonalSet(ctx->A,D,is);
136:   return(0);
137: }

141: PetscErrorCode MatDestroy_SMF(Mat mat)
142: {
143:   PetscErrorCode   ierr;
144:   MatSubMatFreeCtx ctx;

147:   MatShellGetContext(mat,(void **)&ctx);
148:   MatDestroy(&ctx->A);
149:   ISDestroy(&ctx->Rows);
150:   ISDestroy(&ctx->Cols);
151:   VecDestroy(&ctx->VC);
152:   PetscFree(ctx);
153:   return(0);
154: }



160: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
161: {
162:   PetscErrorCode   ierr;
163:   MatSubMatFreeCtx ctx;

166:   MatShellGetContext(mat,(void **)&ctx);
167:   MatView(ctx->A,viewer);
168:   return(0);
169: }

173: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
174: {
175:   PetscErrorCode   ierr;
176:   MatSubMatFreeCtx ctx;

179:   MatShellGetContext(Y,(void **)&ctx);
180:   MatShift(ctx->A,a);
181:   return(0);
182: }

186: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
187: {
188:   PetscErrorCode   ierr;
189:   MatSubMatFreeCtx ctx;

192:   MatShellGetContext(mat,(void **)&ctx);
193:   MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
194:   return(0);
195: }

199: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
200: {
201:   PetscErrorCode    ierr;
202:   MatSubMatFreeCtx  ctx1,ctx2;
203:   PetscBool         flg1,flg2,flg3;

206:   MatShellGetContext(A,(void **)&ctx1);
207:   MatShellGetContext(B,(void **)&ctx2);
208:   ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
209:   ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
210:   if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
211:     *flg=PETSC_FALSE;
212:   } else {
213:     MatEqual(ctx1->A,ctx2->A,&flg1);
214:     if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
215:     else { *flg=PETSC_TRUE;}
216:   }
217:   return(0);
218: }

222: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
223: {
224:   PetscErrorCode   ierr;
225:   MatSubMatFreeCtx ctx;

228:   MatShellGetContext(mat,(void **)&ctx);
229:   MatScale(ctx->A,a);
230:   return(0);
231: }

235: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
236: {
238:   PetscFunctionReturn(1);
239: }

243: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
244: {
245:   PetscErrorCode   ierr;
246:   MatSubMatFreeCtx ctx;

249:   MatShellGetContext(mat,(void **)&ctx);
250:   MatGetDiagonal(ctx->A,v);
251:   return(0);
252: }

256: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
257: {
258:   MatSubMatFreeCtx ctx;
259:   PetscErrorCode   ierr;

262:   MatShellGetContext(M,(void **)&ctx);
263:   MatGetRowMax(ctx->A,D,NULL);
264:   return(0);
265: }

269: PetscErrorCode MatGetSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
270: {
272:   PetscInt       i;

275:   if (scall == MAT_INITIAL_MATRIX) {
276:     PetscMalloc1(n+1,B );
277:   }

279:   for ( i=0; i<n; i++ ) {
280:     MatGetSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
281:   }
282:   return(0);
283: }

287: PetscErrorCode MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
288:                         Mat *newmat)
289: {
290:   PetscErrorCode   ierr;
291:   MatSubMatFreeCtx ctx;

294:   MatShellGetContext(mat,(void **)&ctx);
295:   if (newmat){
296:     MatDestroy(&*newmat);
297:   }
298:   MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
299:   return(0);
300: }

304: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
305: {
306:   PetscErrorCode   ierr;
307:   MatSubMatFreeCtx ctx;

310:   MatShellGetContext(mat,(void **)&ctx);
311:   MatGetRow(ctx->A,row,ncols,cols,vals);
312:   return(0);
313: }

317: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
318: {
319:   PetscErrorCode   ierr;
320:   MatSubMatFreeCtx ctx;

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

330: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
331: {
332:   PetscErrorCode   ierr;
333:   MatSubMatFreeCtx ctx;

336:   MatShellGetContext(mat,(void **)&ctx);
337:   MatGetColumnVector(ctx->A,Y,col);
338:   return(0);
339: }

343: PetscErrorCode MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat)
344: {
345:   PetscErrorCode    ierr;
346:   PetscMPIInt       size;
347:   MatSubMatFreeCtx  ctx;

350:   MatShellGetContext(mat,(void **)&ctx);
351:   MPI_Comm_size(((PetscObject)mat)->comm,&size);
352:   PetscFunctionReturn(1);
353: }

357: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
358: {
359:   PetscErrorCode    ierr;
360:   MatSubMatFreeCtx  ctx;

363:   MatShellGetContext(mat,(void **)&ctx);
364:   if (type == NORM_FROBENIUS) {
365:     *norm = 1.0;
366:   } else if (type == NORM_1 || type == NORM_INFINITY) {
367:     *norm = 1.0;
368:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
369:   return(0);
370: }