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=((PetscObject)mat)->comm;
34: MatSubMatFreeCtx ctx;
35: int info;
36: PetscInt mloc,nloc,m,n;
39: info = PetscNew(_p_MatSubMatFreeCtx,&ctx);CHKERRQ(info);
41: ctx->A=mat;
42: // ctx->Row=Row;
43: // ctx->Col=Col;
45: info = MatGetSize(mat,&m,&n);CHKERRQ(info);
46: info = MatGetLocalSize(mat,&mloc,&nloc);CHKERRQ(info);
48: info = VecCreateMPI(comm,nloc,n,&ctx->VC);CHKERRQ(info);
49: // info = ISCreateComplement(Col, ctx->VC, &ctx->ColComplement);CHKERRQ(info);
50: // ctx->RowComplement=ctx->ColComplement;
51: ctx->VR=ctx->VC;
52: info = PetscObjectReference((PetscObject)mat);CHKERRQ(info);
54: info=ISCreateComplement(RowMask,ctx->VC,&ctx->RowComplement);CHKERRQ(info);
55: info=ISCreateComplement(ColMask,ctx->VC,&ctx->ColComplement);CHKERRQ(info);
56: /*
57: info = PetscObjectReference((PetscObject)ctx->RowComplement);CHKERRQ(info);
58: info = PetscObjectReference((PetscObject)ctx->ColComplement);CHKERRQ(info);
59: */
60: info = MatCreateShell(comm,mloc,nloc,m,n,ctx,J);CHKERRQ(info);
62: // info = MatShellSetOperation(*J,MATOP_GET_ROW,(void(*)())MatGetRow_SMF);CHKERRQ(info);
63: // info = MatShellSetOperation(*J,MATOP_RESTORE_ROW,(void(*)())MatRestoreRow_SMF);CHKERRQ(info);
64: info = MatShellSetOperation(*J,MATOP_MULT,(void(*)())MatMult_SMF);CHKERRQ(info);
65: info = MatShellSetOperation(*J,MATOP_DESTROY,(void(*)())MatDestroy_SMF);CHKERRQ(info);
66: info = MatShellSetOperation(*J,MATOP_VIEW,(void(*)())MatView_SMF);CHKERRQ(info);
67: info = MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)())MatMultTranspose_SMF);CHKERRQ(info);
68: info = MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)())MatDiagonalSet_SMF);CHKERRQ(info);
69: info = MatShellSetOperation(*J,MATOP_SHIFT,(void(*)())MatShift_SMF);CHKERRQ(info);
70: info = MatShellSetOperation(*J,MATOP_EQUAL,(void(*)())MatEqual_SMF);CHKERRQ(info);
71: info = MatShellSetOperation(*J,MATOP_SCALE,(void(*)())MatScale_SMF);CHKERRQ(info);
72: info = MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)())MatTranspose_SMF);CHKERRQ(info);
73: info = MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)())MatGetDiagonal_SMF);CHKERRQ(info);
74: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRICES,(void(*)())MatGetSubMatrices_SMF);CHKERRQ(info);
75: info = MatShellSetOperation(*J,MATOP_NORM,(void(*)())MatNorm_SMF);CHKERRQ(info);
76: info = MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)())MatDuplicate_SMF);CHKERRQ(info);
77: info = MatShellSetOperation(*J,MATOP_GET_SUBMATRIX,(void(*)())MatGetSubMatrix_SMF);CHKERRQ(info);
78: info = MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)())MatDuplicate_SMF);CHKERRQ(info);
80: info = PetscLogObjectParent(mat,*J); CHKERRQ(info);
82: return(0);
83: }
87: int MatSMFResetRowColumn(Mat mat,IS RowMask,IS ColMask){
88: MatSubMatFreeCtx ctx;
89: int info;
92: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
93: info = PetscObjectReference((PetscObject)RowMask);CHKERRQ(info);
94: info = PetscObjectReference((PetscObject)ColMask);CHKERRQ(info);
95: info = ISDestroy(ctx->RowComplement);CHKERRQ(info);
96: info = ISDestroy(ctx->ColComplement);CHKERRQ(info);
97: ctx->ColComplement=ColMask;
98: ctx->RowComplement=RowMask;
99: return(0);
100: }
104: int MatMult_SMF(Mat mat,Vec a,Vec y)
105: {
106: MatSubMatFreeCtx ctx;
107: int info;
110: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
111: info = VecCopy(a,ctx->VR);CHKERRQ(info);
112: info = VecISSetToConstant(ctx->ColComplement,0.0,ctx->VR);CHKERRQ(info);
113: info = MatMult(ctx->A,ctx->VR,y);CHKERRQ(info);
114: info = VecISSetToConstant(ctx->RowComplement,0.0,y);CHKERRQ(info);
115: return(0);
116: }
120: int MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
121: {
122: MatSubMatFreeCtx ctx;
123: int info;
126: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
127: info = VecCopy(a,ctx->VC);CHKERRQ(info);
128: info = VecISSetToConstant(ctx->RowComplement,0.0,ctx->VC);CHKERRQ(info);
129: info = MatMultTranspose(ctx->A,ctx->VC,y);CHKERRQ(info);
130: info = VecISSetToConstant(ctx->ColComplement,0.0,y);CHKERRQ(info);
131: return(0);
132: }
136: int MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
137: {
138: MatSubMatFreeCtx ctx;
139: int info;
142: info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
143: info = MatDiagonalSet(ctx->A,D,is);
144: return(0);
145: }
149: int MatDestroy_SMF(Mat mat)
150: {
151: int info;
152: MatSubMatFreeCtx ctx;
155: info=MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
156: // info=ISDestroy(ctx->Row);CHKERRQ(info);
157: // info=ISDestroy(ctx->Col);CHKERRQ(info);
158: info =MatDestroy(ctx->A);CHKERRQ(info);
159: info=ISDestroy(ctx->RowComplement);CHKERRQ(info);
160: info=ISDestroy(ctx->ColComplement);CHKERRQ(info);
161: info=VecDestroy(ctx->VC);CHKERRQ(info);
162: info = PetscFree(ctx); CHKERRQ(info);
163: return(0);
164: }
170: int MatView_SMF(Mat mat,PetscViewer viewer)
171: {
172: int info;
173: MatSubMatFreeCtx ctx;
176: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
177: info = MatView(ctx->A,viewer);CHKERRQ(info);
178: // info = ISView(ctx->Row,viewer);CHKERRQ(info);
179: // info = ISView(ctx->Col,viewer);CHKERRQ(info);
180: return(0);
181: }
185: int MatShift_SMF(Mat Y, PetscScalar a)
186: {
187: int info;
188: MatSubMatFreeCtx ctx;
191: info = MatShellGetContext(Y,(void **)&ctx);CHKERRQ(info);
192: info = MatShift(ctx->A,a);CHKERRQ(info);
193: return(0);
194: }
198: int MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
199: {
200: int info;
201: MatSubMatFreeCtx ctx;
204: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
205: info = MatCreateSubMatrixFree(ctx->A,ctx->RowComplement,ctx->ColComplement,M);CHKERRQ(info);
206: return(0);
207: }
211: int MatEqual_SMF(Mat A,Mat B,PetscTruth *flg)
212: {
213: int info;
214: MatSubMatFreeCtx ctx1,ctx2;
215: PetscTruth flg1,flg2,flg3;
218: info = MatShellGetContext(A,(void **)&ctx1);CHKERRQ(info);
219: info = MatShellGetContext(B,(void **)&ctx2);CHKERRQ(info);
220: info = ISEqual(ctx1->RowComplement,ctx2->RowComplement,&flg2);CHKERRQ(info);
221: info = ISEqual(ctx1->ColComplement,ctx2->ColComplement,&flg3);CHKERRQ(info);
222: if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
223: *flg=PETSC_FALSE;
224: } else {
225: info = MatEqual(ctx1->A,ctx2->A,&flg1);CHKERRQ(info);
226: if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
227: else { *flg=PETSC_TRUE;}
228: }
229: return(0);
230: }
234: int MatScale_SMF(Mat mat, PetscScalar a)
235: {
236: int info;
237: MatSubMatFreeCtx ctx;
240: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
241: info = MatScale(ctx->A,a);CHKERRQ(info);
242: return(0);
243: }
247: int MatTranspose_SMF(Mat mat,Mat *B)
248: {
250: PetscFunctionReturn(1);
251: }
255: int MatGetDiagonal_SMF(Mat mat,Vec v)
256: {
257: int info;
258: MatSubMatFreeCtx ctx;
261: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
262: info = MatGetDiagonal(ctx->A,v);CHKERRQ(info);
263: // info = VecISSetToConstant(ctx->RowComplement,0.0,v);CHKERRQ(info);
264: return(0);
265: }
269: int MatGetRowMax_SMF(Mat M, Vec D)
270: {
271: MatSubMatFreeCtx ctx;
272: int info;
275: info = MatShellGetContext(M,(void **)&ctx);CHKERRQ(info);
276: info = MatGetRowMax(ctx->A,D,PETSC_NULL);
277: // info = VecISSetToConstant(ctx->RowComplement,0.0,D);CHKERRQ(info);
278: return(0);
279: }
284: int MatGetSubMatrices_SMF(Mat A,int n, IS *irow,IS *icol,MatReuse scall,Mat **B)
285: {
286: int info,i;
289: if (scall == MAT_INITIAL_MATRIX) {
290: info = PetscMalloc( (n+1)*sizeof(Mat),B );CHKERRQ(info);
291: }
293: for ( i=0; i<n; i++ ) {
294: info = MatGetSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);CHKERRQ(info);
295: }
296: return(0);
297: }
301: int MatGetSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
302: Mat *newmat)
303: {
304: int info;
305: MatSubMatFreeCtx ctx;
308: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
309: if (newmat){
310: info=MatDestroy(*newmat);CHKERRQ(info);
311: }
312: info = MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);CHKERRQ(info);
313: return(0);
314: }
318: int MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
319: {
320: int info;
321: MatSubMatFreeCtx ctx;
324: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
325: info = MatGetRow(ctx->A,row,ncols,cols,vals);CHKERRQ(info);
326: return(0);
327: }
331: int MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
332: {
333: int info;
334: MatSubMatFreeCtx ctx;
337: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
338: info = MatRestoreRow(ctx->A,row,ncols,cols,vals);CHKERRQ(info);
339: return(0);
340: }
344: int MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
345: {
346: int info;
347: MatSubMatFreeCtx ctx;
350: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
351: info = MatGetColumnVector(ctx->A,Y,col);CHKERRQ(info);
352: // info = VecISSetToConstant(ctx->RowComplement,0.0,Y);CHKERRQ(info);
353: return(0);
354: }
358: int MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat)
359: {
360: int info;
361: int size;
362: MatSubMatFreeCtx ctx;
365: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
366: MPI_Comm_size(((PetscObject)mat)->comm,&size);
367: PetscFunctionReturn(1);
368: }
372: int MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
373: {
374: int info;
375: MatSubMatFreeCtx ctx;
378: info = MatShellGetContext(mat,(void **)&ctx);CHKERRQ(info);
380: if (type == NORM_FROBENIUS) {
381: *norm = 1.0;
382: } else if (type == NORM_1 || type == NORM_INFINITY) {
383: *norm = 1.0;
384: } else {
385: SETERRQ(PETSC_ERR_SUP,"No two norm");
386: }
387: return(0);
388: }