Actual source code: submatfree.c
1: #include "submatfree.h" /*I "mat.h" I*/
3: int ISCreateComplement(IS, Vec, IS *);
4: int VecISSetToConstant(IS, PetscScalar, Vec);
9: /*@C
10: MatCreateSubMatrixFree - Creates a reduced matrix by masking a
11: full matrix.
13: Collective on matrix
15: Input Parameters:
16: + mat - matrix of arbitrary type
17: . RowMask - the rows that will be zero
18: - ColMask - the columns that will be zero
20: Output Parameters:
21: . J - New matrix
23: Notes:
24: The user provides the input data and is responsible for destroying
25: this data after matrix J has been destroyed.
26:
27: Level: developer
29: .seealso: MatCreate()
30: @*/
31: int MatCreateSubMatrixFree(Mat mat,IS RowMask, IS ColMask, Mat *J)
32: {
33: MPI_Comm comm=mat->comm;
34: MatSubMatFreeCtx ctx;
35: int info,mloc,nloc,m,n;
39: info = PetscNew(_p_MatSubMatFreeCtx,&ctx);CHKERRQ(info);
41: ctx->A=mat;
42: // ctx->Row=Row;
43: // ctx->Col=Col;
44: ctx->ColComplement=ColMask;
45: ctx->RowComplement=RowMask;
47: info = MatGetSize(mat,&m,&n);CHKERRQ(info);
48: info = MatGetLocalSize(mat,&mloc,&nloc);CHKERRQ(info);
50: info = VecCreateMPI(comm,nloc,n,&ctx->VC);CHKERRQ(info);
51: // info = ISCreateComplement(Col, ctx->VC, &ctx->ColComplement);CHKERRQ(info);
52: // ctx->RowComplement=ctx->ColComplement;
54: ctx->VR=ctx->VC;
55: info = PetscObjectReference((PetscObject)mat);CHKERRQ(info);
56: info = PetscObjectReference((PetscObject)RowMask);CHKERRQ(info);
57: info = PetscObjectReference((PetscObject)ColMask);CHKERRQ(info);
59: info = MatCreateShell(comm,mloc,nloc,m,n,ctx,J);CHKERRQ(info);
61: // info = MatShellSetOperation(*J,MATOP_GET_ROW,(void(*)())MatGetRow_SMF);CHKERRQ(info);
62: // info = MatShellSetOperation(*J,MATOP_RESTORE_ROW,(void(*)())MatRestoreRow_SMF);CHKERRQ(info);
63: info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_SMF);CHKERRQ(info);
64: info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_SMF);CHKERRQ(info);
65: info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_SMF);CHKERRQ(info);
66: info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_SMF);CHKERRQ(info);
67: info = MatShellSetOperation(*J,MATOP_DIAGONAL_SHIFT,(void(*)())MatDiagonalSet_SMF);CHKERRQ(info);
68: info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_SMF);CHKERRQ(info);
69: info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_SMF);CHKERRQ(info);
70: info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_SMF);CHKERRQ(info);
71: info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_SMF);CHKERRQ(info);
72: info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_SMF);CHKERRQ(info);
73: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_SMF);CHKERRQ(info);
74: info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_SMF);CHKERRQ(info);
75: info = MatShellSetOperation(*J,MATOP_DUPLCIATE,(void(*)())MatDuplicate_SMF);CHKERRQ(info);
76: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_SMF);CHKERRQ(info);
77: info = MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)())MatDuplicate_SMF);CHKERRQ(info);
79: info = PetscLogObjectParent(mat,*J); CHKERRQ(info);
81: return(0);
82: }
86: int MatSMFResetRowColumn(Mat mat,IS RowMask,IS ColMask){
87: MatSubMatFreeCtx ctx;
88: int info;
91: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
92: info = PetscObjectReference((PetscObject)RowMask);CHKERRQ(info);
93: info = PetscObjectReference((PetscObject)ColMask);CHKERRQ(info);
94: info = ISDestroy(ctx->RowComplement);CHKERRQ(info);
95: info = ISDestroy(ctx->ColComplement);CHKERRQ(info);
96: ctx->ColComplement=ColMask;
97: ctx->RowComplement=RowMask;
98: return(0);
99: }
103: int MatMult_SMF(Mat mat,Vec a,Vec y)
104: {
105: MatSubMatFreeCtx ctx;
106: int info;
109: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
110: info = VecCopy(a,ctx->VR);CHKERRQ(info);
111: info = VecISSetToConstant(ctx->ColComplement,0.0,ctx->VR);CHKERRQ(info);
112: info = MatMult(ctx->A,ctx->VR,y);CHKERRQ(info);
113: info = VecISSetToConstant(ctx->RowComplement,0.0,y);CHKERRQ(info);
114: return(0);
115: }
119: int MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
120: {
121: MatSubMatFreeCtx ctx;
122: int info;
125: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
126: info = VecCopy(a,ctx->VC);CHKERRQ(info);
127: info = VecISSetToConstant(ctx->RowComplement,0.0,ctx->VC);CHKERRQ(info);
128: info = MatMultTranspose(ctx->A,ctx->VC,y);CHKERRQ(info);
129: info = VecISSetToConstant(ctx->ColComplement,0.0,y);CHKERRQ(info);
130: return(0);
131: }
135: int MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
136: {
137: MatSubMatFreeCtx ctx;
138: int info;
141: info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
142: info = MatDiagonalSet(ctx->A,D,is);
143: return(0);
144: }
148: int MatDestroy_SMF(Mat mat)
149: {
150: int info;
151: MatSubMatFreeCtx ctx;
154: info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
155: // info=ISDestroy(ctx->Row);CHKERRQ(info);
156: // info=ISDestroy(ctx->Col);CHKERRQ(info);
157: info =MatDestroy(ctx->A);CHKERRQ(info);
158: info=ISDestroy(ctx->RowComplement);CHKERRQ(info);
159: info=ISDestroy(ctx->ColComplement);CHKERRQ(info);
160: info=VecDestroy(ctx->VC);CHKERRQ(info);
161: PetscFree(ctx);
162: return(0);
163: }
169: int MatView_SMF(Mat mat,PetscViewer viewer)
170: {
171: int info;
172: MatSubMatFreeCtx ctx;
175: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
176: info = MatView(ctx->A,viewer);CHKERRQ(info);
177: // info = ISView(ctx->Row,viewer);CHKERRQ(info);
178: // info = ISView(ctx->Col,viewer);CHKERRQ(info);
179: return(0);
180: }
184: int MatShift_SMF(Mat Y, PetscScalar a)
185: {
186: int info;
187: MatSubMatFreeCtx ctx;
190: info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
191: info = MatShift(ctx->A,a);CHKERRQ(info);
192: return(0);
193: }
197: int MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
198: {
199: int info;
200: MatSubMatFreeCtx ctx;
203: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
204: info = MatCreateSubMatrixFree(ctx->A,ctx->RowComplement,ctx->ColComplement,M);CHKERRQ(info);
205: return(0);
206: }
210: int MatEqual_SMF(Mat A,Mat B,PetscTruth *flg)
211: {
212: int info;
213: MatSubMatFreeCtx ctx1,ctx2;
214: PetscTruth flg1,flg2,flg3;
217: info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
218: info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
219: info = ISEqual(ctx1->RowComplement,ctx2->RowComplement,&flg2);CHKERRQ(info);
220: info = ISEqual(ctx1->ColComplement,ctx2->ColComplement,&flg3);CHKERRQ(info);
221: if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
222: *flg=PETSC_FALSE;
223: } else {
224: info = MatEqual(ctx1->A,ctx2->A,&flg1);CHKERRQ(info);
225: if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
226: else { *flg=PETSC_TRUE;}
227: }
228: return(0);
229: }
233: int MatScale_SMF(Mat mat, PetscScalar a)
234: {
235: int info;
236: MatSubMatFreeCtx ctx;
239: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
240: info = MatScale(ctx->A,a);CHKERRQ(info);
241: return(0);
242: }
246: int MatTranspose_SMF(Mat mat,Mat *B)
247: {
249: PetscFunctionReturn(1);
250: }
254: int MatGetDiagonal_SMF(Mat mat,Vec v)
255: {
256: int info;
257: MatSubMatFreeCtx ctx;
260: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
261: info = MatGetDiagonal(ctx->A,v);CHKERRQ(info);
262: // info = VecISSetToConstant(ctx->RowComplement,0.0,v);CHKERRQ(info);
263: return(0);
264: }
268: int MatGetRowMax_SMF(Mat M, Vec D)
269: {
270: MatSubMatFreeCtx ctx;
271: int info;
274: info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
275: info = MatGetRowMax(ctx->A,D);
276: // info = VecISSetToConstant(ctx->RowComplement,0.0,D);CHKERRQ(info);
277: return(0);
278: }
283: int MatGetSubMatrices_SMF(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
284: {
285: int info,i;
288: if (scall == MAT_INITIAL_MATRIX) {
289: info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
290: }
292: for ( i=0; i<n; i++ ) {
293: info = MatGetSubMatrix_SMF(A,irow[i],icol[i],PETSC_DECIDE,scall,&(*B)[i]);CHKERRQ(info);
294: }
295: return(0);
296: }
300: int MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,int csize,MatReuse cll,
301: Mat *newmat)
302: {
303: int info;
304: MatSubMatFreeCtx ctx;
307: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
308: if (newmat){
309: info=MatDestroy(*newmat);CHKERRQ(info);
310: }
311: info = MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);CHKERRQ(info);
312: return(0);
313: }
317: int MatGetRow_SMF(Mat mat,int row,int *ncols,const int **cols,const PetscScalar **vals)
318: {
319: int info;
320: MatSubMatFreeCtx ctx;
323: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
324: info = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(info);
325: return(0);
326: }
330: int MatRestoreRow_SMF(Mat mat,int row,int *ncols,const int **cols,const PetscScalar **vals)
331: {
332: int info;
333: MatSubMatFreeCtx ctx;
336: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
337: info = MatRestoreRow(ctx->A,row,ncols,cols,vals);CHKERRQ(info);
338: return(0);
339: }
343: int MatGetColumnVector_SMF(Mat mat,Vec Y, int col)
344: {
345: int info;
346: MatSubMatFreeCtx ctx;
349: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
350: info = MatGetColumnVector(ctx->A,Y,col);CHKERRQ(info);
351: // info = VecISSetToConstant(ctx->RowComplement,0.0,Y);CHKERRQ(info);
352: return(0);
353: }
357: int MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat)
358: {
359: int info,size;
360: MatSubMatFreeCtx ctx;
363: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
364: MPI_Comm_size(mat->comm,&size);
365: PetscFunctionReturn(1);
366: }
370: int MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
371: {
372: int info;
373: MatSubMatFreeCtx ctx;
376: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
378: if (type == NORM_FROBENIUS) {
379: *norm = 1.0;
380: } else if (type == NORM_1 || type == NORM_INFINITY) {
381: *norm = 1.0;
382: } else {
383: SETERRQ(PETSC_ERR_SUP,"No two norm");
384: }
385: return(0);
386: }