Actual source code: submatfree.c
petsc-3.8.4 2018-03-24
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 user provides the input data and is responsible for destroying
20: this data after matrix J has been destroyed.
22: Level: developer
24: .seealso: MatCreate()
25: @*/
26: PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J)
27: {
28: MPI_Comm comm=PetscObjectComm((PetscObject)mat);
29: MatSubMatFreeCtx ctx;
30: PetscErrorCode ierr;
31: PetscMPIInt size;
32: PetscInt mloc,nloc,m,n;
35: PetscNew(&ctx);
36: ctx->A=mat;
37: MatGetSize(mat,&m,&n);
38: MatGetLocalSize(mat,&mloc,&nloc);
39: MPI_Comm_size(comm,&size);
40: if (size == 1) {
41: VecCreateSeq(comm,n,&ctx->VC);
42: } else {
43: VecCreateMPI(comm,nloc,n,&ctx->VC);
44: }
45: ctx->VR=ctx->VC;
46: PetscObjectReference((PetscObject)mat);
49: ctx->Rows = Rows;
50: ctx->Cols = Cols;
51: PetscObjectReference((PetscObject)Rows);
52: PetscObjectReference((PetscObject)Cols);
53: MatCreateShell(comm,mloc,nloc,m,n,ctx,J);
55: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);
56: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);
57: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);
58: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);
59: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);
60: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);
61: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);
62: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);
63: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);
64: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);
65: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_SMF);
66: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);
67: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);
68: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_SMF);
69: MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);
71: PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
72: return(0);
73: }
75: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){
76: MatSubMatFreeCtx ctx;
77: PetscErrorCode ierr;
80: MatShellGetContext(mat,(void **)&ctx);
81: ISDestroy(&ctx->Rows);
82: ISDestroy(&ctx->Cols);
83: PetscObjectReference((PetscObject)Rows);
84: PetscObjectReference((PetscObject)Cols);
85: ctx->Cols=Cols;
86: ctx->Rows=Rows;
87: return(0);
88: }
90: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
91: {
92: MatSubMatFreeCtx ctx;
93: PetscErrorCode ierr;
96: MatShellGetContext(mat,(void **)&ctx);
97: VecCopy(a,ctx->VR);
98: VecISSet(ctx->VR,ctx->Cols,0.0);
99: MatMult(ctx->A,ctx->VR,y);
100: VecISSet(y,ctx->Rows,0.0);
101: return(0);
102: }
104: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
105: {
106: MatSubMatFreeCtx ctx;
107: PetscErrorCode ierr;
110: MatShellGetContext(mat,(void **)&ctx);
111: VecCopy(a,ctx->VC);
112: VecISSet(ctx->VC,ctx->Rows,0.0);
113: MatMultTranspose(ctx->A,ctx->VC,y);
114: VecISSet(y,ctx->Cols,0.0);
115: return(0);
116: }
118: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
119: {
120: MatSubMatFreeCtx ctx;
121: PetscErrorCode ierr;
124: MatShellGetContext(M,(void **)&ctx);
125: MatDiagonalSet(ctx->A,D,is);
126: return(0);
127: }
129: PetscErrorCode MatDestroy_SMF(Mat mat)
130: {
131: PetscErrorCode ierr;
132: MatSubMatFreeCtx ctx;
135: MatShellGetContext(mat,(void **)&ctx);
136: MatDestroy(&ctx->A);
137: ISDestroy(&ctx->Rows);
138: ISDestroy(&ctx->Cols);
139: VecDestroy(&ctx->VC);
140: PetscFree(ctx);
141: return(0);
142: }
146: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
147: {
148: PetscErrorCode ierr;
149: MatSubMatFreeCtx ctx;
152: MatShellGetContext(mat,(void **)&ctx);
153: MatView(ctx->A,viewer);
154: return(0);
155: }
157: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
158: {
159: PetscErrorCode ierr;
160: MatSubMatFreeCtx ctx;
163: MatShellGetContext(Y,(void **)&ctx);
164: MatShift(ctx->A,a);
165: return(0);
166: }
168: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
169: {
170: PetscErrorCode ierr;
171: MatSubMatFreeCtx ctx;
174: MatShellGetContext(mat,(void **)&ctx);
175: MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
176: return(0);
177: }
179: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
180: {
181: PetscErrorCode ierr;
182: MatSubMatFreeCtx ctx1,ctx2;
183: PetscBool flg1,flg2,flg3;
186: MatShellGetContext(A,(void **)&ctx1);
187: MatShellGetContext(B,(void **)&ctx2);
188: ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
189: ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
190: if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
191: *flg=PETSC_FALSE;
192: } else {
193: MatEqual(ctx1->A,ctx2->A,&flg1);
194: if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
195: else { *flg=PETSC_TRUE;}
196: }
197: return(0);
198: }
200: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
201: {
202: PetscErrorCode ierr;
203: MatSubMatFreeCtx ctx;
206: MatShellGetContext(mat,(void **)&ctx);
207: MatScale(ctx->A,a);
208: return(0);
209: }
211: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
212: {
214: PetscFunctionReturn(1);
215: }
217: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
218: {
219: PetscErrorCode ierr;
220: MatSubMatFreeCtx ctx;
223: MatShellGetContext(mat,(void **)&ctx);
224: MatGetDiagonal(ctx->A,v);
225: return(0);
226: }
228: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
229: {
230: MatSubMatFreeCtx ctx;
231: PetscErrorCode ierr;
234: MatShellGetContext(M,(void **)&ctx);
235: MatGetRowMax(ctx->A,D,NULL);
236: return(0);
237: }
239: PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
240: {
242: PetscInt i;
245: if (scall == MAT_INITIAL_MATRIX) {
246: PetscCalloc1(n+1,B );
247: }
249: for ( i=0; i<n; i++ ) {
250: MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
251: }
252: return(0);
253: }
255: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
256: Mat *newmat)
257: {
258: PetscErrorCode ierr;
259: MatSubMatFreeCtx ctx;
262: MatShellGetContext(mat,(void **)&ctx);
263: if (newmat){
264: MatDestroy(&*newmat);
265: }
266: MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
267: return(0);
268: }
270: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
271: {
272: PetscErrorCode ierr;
273: MatSubMatFreeCtx ctx;
276: MatShellGetContext(mat,(void **)&ctx);
277: MatGetRow(ctx->A,row,ncols,cols,vals);
278: return(0);
279: }
281: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
282: {
283: PetscErrorCode ierr;
284: MatSubMatFreeCtx ctx;
287: MatShellGetContext(mat,(void **)&ctx);
288: MatRestoreRow(ctx->A,row,ncols,cols,vals);
289: return(0);
290: }
292: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
293: {
294: PetscErrorCode ierr;
295: MatSubMatFreeCtx ctx;
298: MatShellGetContext(mat,(void **)&ctx);
299: MatGetColumnVector(ctx->A,Y,col);
300: return(0);
301: }
303: PetscErrorCode MatConvert_SMF(Mat mat,MatType newtype,Mat *NewMat)
304: {
305: PetscErrorCode ierr;
306: PetscMPIInt size;
307: MatSubMatFreeCtx ctx;
310: MatShellGetContext(mat,(void **)&ctx);
311: MPI_Comm_size(PetscObjectComm((PetscObject)mat),&size);
312: PetscFunctionReturn(1);
313: }
315: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
316: {
317: PetscErrorCode ierr;
318: MatSubMatFreeCtx ctx;
321: MatShellGetContext(mat,(void **)&ctx);
322: if (type == NORM_FROBENIUS) {
323: *norm = 1.0;
324: } else if (type == NORM_1 || type == NORM_INFINITY) {
325: *norm = 1.0;
326: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
327: return(0);
328: }