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: }