Actual source code: submatfree.c
petsc-3.14.6 2021-03-30
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 caller is responsible for destroying the input objects after matrix J has been destroyed.
21: Level: developer
23: .seealso: MatCreate()
24: @*/
25: PetscErrorCode MatCreateSubMatrixFree(Mat mat,IS Rows, IS Cols, Mat *J)
26: {
27: MPI_Comm comm=PetscObjectComm((PetscObject)mat);
28: MatSubMatFreeCtx ctx;
29: PetscErrorCode ierr;
30: PetscInt mloc,nloc,m,n;
33: PetscNew(&ctx);
34: ctx->A = mat;
35: MatGetSize(mat,&m,&n);
36: MatGetLocalSize(mat,&mloc,&nloc);
37: MatCreateVecs(mat,NULL,&ctx->VC);
38: ctx->VR = ctx->VC;
39: PetscObjectReference((PetscObject)mat);
42: ctx->Rows = Rows;
43: ctx->Cols = Cols;
44: PetscObjectReference((PetscObject)Rows);
45: PetscObjectReference((PetscObject)Cols);
46: MatCreateShell(comm,mloc,nloc,m,n,ctx,J);
47: MatShellSetManageScalingShifts(*J);
48: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);
49: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);
50: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);
51: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);
52: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);
53: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);
54: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);
55: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);
56: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);
57: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);
58: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_SMF);
59: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);
60: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);
61: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_SMF);
62: MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);
64: PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
65: return(0);
66: }
68: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols){
69: MatSubMatFreeCtx ctx;
70: PetscErrorCode ierr;
73: MatShellGetContext(mat,(void **)&ctx);
74: ISDestroy(&ctx->Rows);
75: ISDestroy(&ctx->Cols);
76: PetscObjectReference((PetscObject)Rows);
77: PetscObjectReference((PetscObject)Cols);
78: ctx->Cols=Cols;
79: ctx->Rows=Rows;
80: return(0);
81: }
83: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
84: {
85: MatSubMatFreeCtx ctx;
86: PetscErrorCode ierr;
89: MatShellGetContext(mat,(void **)&ctx);
90: VecCopy(a,ctx->VR);
91: VecISSet(ctx->VR,ctx->Cols,0.0);
92: MatMult(ctx->A,ctx->VR,y);
93: VecISSet(y,ctx->Rows,0.0);
94: return(0);
95: }
97: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
98: {
99: MatSubMatFreeCtx ctx;
100: PetscErrorCode ierr;
103: MatShellGetContext(mat,(void **)&ctx);
104: VecCopy(a,ctx->VC);
105: VecISSet(ctx->VC,ctx->Rows,0.0);
106: MatMultTranspose(ctx->A,ctx->VC,y);
107: VecISSet(y,ctx->Cols,0.0);
108: return(0);
109: }
111: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
112: {
113: MatSubMatFreeCtx ctx;
114: PetscErrorCode ierr;
117: MatShellGetContext(M,(void **)&ctx);
118: MatDiagonalSet(ctx->A,D,is);
119: return(0);
120: }
122: PetscErrorCode MatDestroy_SMF(Mat mat)
123: {
124: PetscErrorCode ierr;
125: MatSubMatFreeCtx ctx;
128: MatShellGetContext(mat,(void **)&ctx);
129: MatDestroy(&ctx->A);
130: ISDestroy(&ctx->Rows);
131: ISDestroy(&ctx->Cols);
132: VecDestroy(&ctx->VC);
133: PetscFree(ctx);
134: return(0);
135: }
139: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
140: {
141: PetscErrorCode ierr;
142: MatSubMatFreeCtx ctx;
145: MatShellGetContext(mat,(void **)&ctx);
146: MatView(ctx->A,viewer);
147: return(0);
148: }
150: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
151: {
152: PetscErrorCode ierr;
153: MatSubMatFreeCtx ctx;
156: MatShellGetContext(Y,(void **)&ctx);
157: MatShift(ctx->A,a);
158: return(0);
159: }
161: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
162: {
163: PetscErrorCode ierr;
164: MatSubMatFreeCtx ctx;
167: MatShellGetContext(mat,(void **)&ctx);
168: MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
169: return(0);
170: }
172: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
173: {
174: PetscErrorCode ierr;
175: MatSubMatFreeCtx ctx1,ctx2;
176: PetscBool flg1,flg2,flg3;
179: MatShellGetContext(A,(void **)&ctx1);
180: MatShellGetContext(B,(void **)&ctx2);
181: ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
182: ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
183: if (flg2==PETSC_FALSE || flg3==PETSC_FALSE){
184: *flg=PETSC_FALSE;
185: } else {
186: MatEqual(ctx1->A,ctx2->A,&flg1);
187: if (flg1==PETSC_FALSE){ *flg=PETSC_FALSE;}
188: else { *flg=PETSC_TRUE;}
189: }
190: return(0);
191: }
193: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
194: {
195: PetscErrorCode ierr;
196: MatSubMatFreeCtx ctx;
199: MatShellGetContext(mat,(void **)&ctx);
200: MatScale(ctx->A,a);
201: return(0);
202: }
204: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
205: {
207: PetscFunctionReturn(1);
208: }
210: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
211: {
212: PetscErrorCode ierr;
213: MatSubMatFreeCtx ctx;
216: MatShellGetContext(mat,(void **)&ctx);
217: MatGetDiagonal(ctx->A,v);
218: return(0);
219: }
221: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
222: {
223: MatSubMatFreeCtx ctx;
224: PetscErrorCode ierr;
227: MatShellGetContext(M,(void **)&ctx);
228: MatGetRowMax(ctx->A,D,NULL);
229: return(0);
230: }
232: PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
233: {
235: PetscInt i;
238: if (scall == MAT_INITIAL_MATRIX) {
239: PetscCalloc1(n+1,B);
240: }
242: for (i=0; i<n; i++) {
243: MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
244: }
245: return(0);
246: }
248: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
249: Mat *newmat)
250: {
251: PetscErrorCode ierr;
252: MatSubMatFreeCtx ctx;
255: MatShellGetContext(mat,(void **)&ctx);
256: if (newmat){
257: MatDestroy(&*newmat);
258: }
259: MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
260: return(0);
261: }
263: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
264: {
265: PetscErrorCode ierr;
266: MatSubMatFreeCtx ctx;
269: MatShellGetContext(mat,(void **)&ctx);
270: MatGetRow(ctx->A,row,ncols,cols,vals);
271: return(0);
272: }
274: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
275: {
276: PetscErrorCode ierr;
277: MatSubMatFreeCtx ctx;
280: MatShellGetContext(mat,(void **)&ctx);
281: MatRestoreRow(ctx->A,row,ncols,cols,vals);
282: return(0);
283: }
285: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
286: {
287: PetscErrorCode ierr;
288: MatSubMatFreeCtx ctx;
291: MatShellGetContext(mat,(void **)&ctx);
292: MatGetColumnVector(ctx->A,Y,col);
293: return(0);
294: }
296: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
297: {
298: PetscErrorCode ierr;
299: MatSubMatFreeCtx ctx;
302: MatShellGetContext(mat,(void **)&ctx);
303: if (type == NORM_FROBENIUS) {
304: *norm = 1.0;
305: } else if (type == NORM_1 || type == NORM_INFINITY) {
306: *norm = 1.0;
307: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
308: return(0);
309: }