Actual source code: submatfree.c
petsc-3.6.1 2015-08-06
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: }