Actual source code: submatfree.c

  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:   PetscInt         mloc,nloc,m,n;

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

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

 61:   PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
 62:   return 0;
 63: }

 65: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols)
 66: {
 67:   MatSubMatFreeCtx ctx;

 69:   MatShellGetContext(mat,&ctx);
 70:   ISDestroy(&ctx->Rows);
 71:   ISDestroy(&ctx->Cols);
 72:   PetscObjectReference((PetscObject)Rows);
 73:   PetscObjectReference((PetscObject)Cols);
 74:   ctx->Cols=Cols;
 75:   ctx->Rows=Rows;
 76:   return 0;
 77: }

 79: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
 80: {
 81:   MatSubMatFreeCtx ctx;

 83:   MatShellGetContext(mat,&ctx);
 84:   VecCopy(a,ctx->VR);
 85:   VecISSet(ctx->VR,ctx->Cols,0.0);
 86:   MatMult(ctx->A,ctx->VR,y);
 87:   VecISSet(y,ctx->Rows,0.0);
 88:   return 0;
 89: }

 91: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
 92: {
 93:   MatSubMatFreeCtx ctx;

 95:   MatShellGetContext(mat,&ctx);
 96:   VecCopy(a,ctx->VC);
 97:   VecISSet(ctx->VC,ctx->Rows,0.0);
 98:   MatMultTranspose(ctx->A,ctx->VC,y);
 99:   VecISSet(y,ctx->Cols,0.0);
100:   return 0;
101: }

103: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
104: {
105:   MatSubMatFreeCtx ctx;

107:   MatShellGetContext(M,&ctx);
108:   MatDiagonalSet(ctx->A,D,is);
109:   return 0;
110: }

112: PetscErrorCode MatDestroy_SMF(Mat mat)
113: {
114:   MatSubMatFreeCtx ctx;

116:   MatShellGetContext(mat,&ctx);
117:   MatDestroy(&ctx->A);
118:   ISDestroy(&ctx->Rows);
119:   ISDestroy(&ctx->Cols);
120:   VecDestroy(&ctx->VC);
121:   PetscFree(ctx);
122:   return 0;
123: }

125: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
126: {
127:   MatSubMatFreeCtx ctx;

129:   MatShellGetContext(mat,&ctx);
130:   MatView(ctx->A,viewer);
131:   return 0;
132: }

134: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
135: {
136:   MatSubMatFreeCtx ctx;

138:   MatShellGetContext(Y,&ctx);
139:   MatShift(ctx->A,a);
140:   return 0;
141: }

143: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
144: {
145:   MatSubMatFreeCtx ctx;

147:   MatShellGetContext(mat,&ctx);
148:   MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
149:   return 0;
150: }

152: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
153: {
154:   MatSubMatFreeCtx  ctx1,ctx2;
155:   PetscBool         flg1,flg2,flg3;

157:   MatShellGetContext(A,&ctx1);
158:   MatShellGetContext(B,&ctx2);
159:   ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
160:   ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
161:   if (flg2==PETSC_FALSE || flg3==PETSC_FALSE) {
162:     *flg=PETSC_FALSE;
163:   } else {
164:     MatEqual(ctx1->A,ctx2->A,&flg1);
165:     if (flg1==PETSC_FALSE) { *flg=PETSC_FALSE;}
166:     else { *flg=PETSC_TRUE;}
167:   }
168:   return 0;
169: }

171: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
172: {
173:   MatSubMatFreeCtx ctx;

175:   MatShellGetContext(mat,&ctx);
176:   MatScale(ctx->A,a);
177:   return 0;
178: }

180: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
181: {
182:   return 1;
183: }

185: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
186: {
187:   MatSubMatFreeCtx ctx;

189:   MatShellGetContext(mat,&ctx);
190:   MatGetDiagonal(ctx->A,v);
191:   return 0;
192: }

194: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
195: {
196:   MatSubMatFreeCtx ctx;

198:   MatShellGetContext(M,&ctx);
199:   MatGetRowMax(ctx->A,D,NULL);
200:   return 0;
201: }

203: PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
204: {
205:   PetscInt       i;

207:   if (scall == MAT_INITIAL_MATRIX) {
208:     PetscCalloc1(n+1,B);
209:   }

211:   for (i=0; i<n; i++) {
212:     MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
213:   }
214:   return 0;
215: }

217: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
218:                         Mat *newmat)
219: {
220:   MatSubMatFreeCtx ctx;

222:   MatShellGetContext(mat,&ctx);
223:   if (newmat) {
224:     MatDestroy(&*newmat);
225:   }
226:   MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
227:   return 0;
228: }

230: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
231: {
232:   MatSubMatFreeCtx ctx;

234:   MatShellGetContext(mat,&ctx);
235:   MatGetRow(ctx->A,row,ncols,cols,vals);
236:   return 0;
237: }

239: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
240: {
241:   MatSubMatFreeCtx ctx;

243:   MatShellGetContext(mat,&ctx);
244:   MatRestoreRow(ctx->A,row,ncols,cols,vals);
245:   return 0;
246: }

248: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
249: {
250:   MatSubMatFreeCtx ctx;

252:   MatShellGetContext(mat,&ctx);
253:   MatGetColumnVector(ctx->A,Y,col);
254:   return 0;
255: }

257: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
258: {
259:   MatSubMatFreeCtx  ctx;

261:   MatShellGetContext(mat,&ctx);
262:   if (type == NORM_FROBENIUS) {
263:     *norm = 1.0;
264:   } else if (type == NORM_1 || type == NORM_INFINITY) {
265:     *norm = 1.0;
266:   } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
267:   return 0;
268: }