Actual source code: submatfree.c
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: PetscInt mloc,nloc,m,n;
31: PetscNew(&ctx);
32: ctx->A = mat;
33: MatGetSize(mat,&m,&n);
34: MatGetLocalSize(mat,&mloc,&nloc);
35: MatCreateVecs(mat,NULL,&ctx->VC);
36: ctx->VR = ctx->VC;
37: PetscObjectReference((PetscObject)mat);
39: ctx->Rows = Rows;
40: ctx->Cols = Cols;
41: PetscObjectReference((PetscObject)Rows);
42: PetscObjectReference((PetscObject)Cols);
43: MatCreateShell(comm,mloc,nloc,m,n,ctx,J);
44: MatShellSetManageScalingShifts(*J);
45: MatShellSetOperation(*J,MATOP_MULT,(void(*)(void))MatMult_SMF);
46: MatShellSetOperation(*J,MATOP_DESTROY,(void(*)(void))MatDestroy_SMF);
47: MatShellSetOperation(*J,MATOP_VIEW,(void(*)(void))MatView_SMF);
48: MatShellSetOperation(*J,MATOP_MULT_TRANSPOSE,(void(*)(void))MatMultTranspose_SMF);
49: MatShellSetOperation(*J,MATOP_DIAGONAL_SET,(void(*)(void))MatDiagonalSet_SMF);
50: MatShellSetOperation(*J,MATOP_SHIFT,(void(*)(void))MatShift_SMF);
51: MatShellSetOperation(*J,MATOP_EQUAL,(void(*)(void))MatEqual_SMF);
52: MatShellSetOperation(*J,MATOP_SCALE,(void(*)(void))MatScale_SMF);
53: MatShellSetOperation(*J,MATOP_TRANSPOSE,(void(*)(void))MatTranspose_SMF);
54: MatShellSetOperation(*J,MATOP_GET_DIAGONAL,(void(*)(void))MatGetDiagonal_SMF);
55: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRICES,(void(*)(void))MatCreateSubMatrices_SMF);
56: MatShellSetOperation(*J,MATOP_NORM,(void(*)(void))MatNorm_SMF);
57: MatShellSetOperation(*J,MATOP_DUPLICATE,(void(*)(void))MatDuplicate_SMF);
58: MatShellSetOperation(*J,MATOP_CREATE_SUBMATRIX,(void(*)(void))MatCreateSubMatrix_SMF);
59: MatShellSetOperation(*J,MATOP_GET_ROW_MAX,(void(*)(void))MatDuplicate_SMF);
61: PetscLogObjectParent((PetscObject)mat,(PetscObject)(*J));
62: return 0;
63: }
65: PetscErrorCode MatSMFResetRowColumn(Mat mat,IS Rows,IS Cols)
66: {
67: MatSubMatFreeCtx ctx;
69: MatShellGetContext(mat,&ctx);
70: ISDestroy(&ctx->Rows);
71: ISDestroy(&ctx->Cols);
72: PetscObjectReference((PetscObject)Rows);
73: PetscObjectReference((PetscObject)Cols);
74: ctx->Cols=Cols;
75: ctx->Rows=Rows;
76: return 0;
77: }
79: PetscErrorCode MatMult_SMF(Mat mat,Vec a,Vec y)
80: {
81: MatSubMatFreeCtx ctx;
83: MatShellGetContext(mat,&ctx);
84: VecCopy(a,ctx->VR);
85: VecISSet(ctx->VR,ctx->Cols,0.0);
86: MatMult(ctx->A,ctx->VR,y);
87: VecISSet(y,ctx->Rows,0.0);
88: return 0;
89: }
91: PetscErrorCode MatMultTranspose_SMF(Mat mat,Vec a,Vec y)
92: {
93: MatSubMatFreeCtx ctx;
95: MatShellGetContext(mat,&ctx);
96: VecCopy(a,ctx->VC);
97: VecISSet(ctx->VC,ctx->Rows,0.0);
98: MatMultTranspose(ctx->A,ctx->VC,y);
99: VecISSet(y,ctx->Cols,0.0);
100: return 0;
101: }
103: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D,InsertMode is)
104: {
105: MatSubMatFreeCtx ctx;
107: MatShellGetContext(M,&ctx);
108: MatDiagonalSet(ctx->A,D,is);
109: return 0;
110: }
112: PetscErrorCode MatDestroy_SMF(Mat mat)
113: {
114: MatSubMatFreeCtx ctx;
116: MatShellGetContext(mat,&ctx);
117: MatDestroy(&ctx->A);
118: ISDestroy(&ctx->Rows);
119: ISDestroy(&ctx->Cols);
120: VecDestroy(&ctx->VC);
121: PetscFree(ctx);
122: return 0;
123: }
125: PetscErrorCode MatView_SMF(Mat mat,PetscViewer viewer)
126: {
127: MatSubMatFreeCtx ctx;
129: MatShellGetContext(mat,&ctx);
130: MatView(ctx->A,viewer);
131: return 0;
132: }
134: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
135: {
136: MatSubMatFreeCtx ctx;
138: MatShellGetContext(Y,&ctx);
139: MatShift(ctx->A,a);
140: return 0;
141: }
143: PetscErrorCode MatDuplicate_SMF(Mat mat,MatDuplicateOption op,Mat *M)
144: {
145: MatSubMatFreeCtx ctx;
147: MatShellGetContext(mat,&ctx);
148: MatCreateSubMatrixFree(ctx->A,ctx->Rows,ctx->Cols,M);
149: return 0;
150: }
152: PetscErrorCode MatEqual_SMF(Mat A,Mat B,PetscBool *flg)
153: {
154: MatSubMatFreeCtx ctx1,ctx2;
155: PetscBool flg1,flg2,flg3;
157: MatShellGetContext(A,&ctx1);
158: MatShellGetContext(B,&ctx2);
159: ISEqual(ctx1->Rows,ctx2->Rows,&flg2);
160: ISEqual(ctx1->Cols,ctx2->Cols,&flg3);
161: if (flg2==PETSC_FALSE || flg3==PETSC_FALSE) {
162: *flg=PETSC_FALSE;
163: } else {
164: MatEqual(ctx1->A,ctx2->A,&flg1);
165: if (flg1==PETSC_FALSE) { *flg=PETSC_FALSE;}
166: else { *flg=PETSC_TRUE;}
167: }
168: return 0;
169: }
171: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
172: {
173: MatSubMatFreeCtx ctx;
175: MatShellGetContext(mat,&ctx);
176: MatScale(ctx->A,a);
177: return 0;
178: }
180: PetscErrorCode MatTranspose_SMF(Mat mat,Mat *B)
181: {
182: return 1;
183: }
185: PetscErrorCode MatGetDiagonal_SMF(Mat mat,Vec v)
186: {
187: MatSubMatFreeCtx ctx;
189: MatShellGetContext(mat,&ctx);
190: MatGetDiagonal(ctx->A,v);
191: return 0;
192: }
194: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
195: {
196: MatSubMatFreeCtx ctx;
198: MatShellGetContext(M,&ctx);
199: MatGetRowMax(ctx->A,D,NULL);
200: return 0;
201: }
203: PetscErrorCode MatCreateSubMatrices_SMF(Mat A,PetscInt n, IS *irow,IS *icol,MatReuse scall,Mat **B)
204: {
205: PetscInt i;
207: if (scall == MAT_INITIAL_MATRIX) {
208: PetscCalloc1(n+1,B);
209: }
211: for (i=0; i<n; i++) {
212: MatCreateSubMatrix_SMF(A,irow[i],icol[i],scall,&(*B)[i]);
213: }
214: return 0;
215: }
217: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat,IS isrow,IS iscol,MatReuse cll,
218: Mat *newmat)
219: {
220: MatSubMatFreeCtx ctx;
222: MatShellGetContext(mat,&ctx);
223: if (newmat) {
224: MatDestroy(&*newmat);
225: }
226: MatCreateSubMatrixFree(ctx->A,isrow,iscol, newmat);
227: return 0;
228: }
230: PetscErrorCode MatGetRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
231: {
232: MatSubMatFreeCtx ctx;
234: MatShellGetContext(mat,&ctx);
235: MatGetRow(ctx->A,row,ncols,cols,vals);
236: return 0;
237: }
239: PetscErrorCode MatRestoreRow_SMF(Mat mat,PetscInt row,PetscInt *ncols,const PetscInt **cols,const PetscScalar **vals)
240: {
241: MatSubMatFreeCtx ctx;
243: MatShellGetContext(mat,&ctx);
244: MatRestoreRow(ctx->A,row,ncols,cols,vals);
245: return 0;
246: }
248: PetscErrorCode MatGetColumnVector_SMF(Mat mat,Vec Y, PetscInt col)
249: {
250: MatSubMatFreeCtx ctx;
252: MatShellGetContext(mat,&ctx);
253: MatGetColumnVector(ctx->A,Y,col);
254: return 0;
255: }
257: PetscErrorCode MatNorm_SMF(Mat mat,NormType type,PetscReal *norm)
258: {
259: MatSubMatFreeCtx ctx;
261: MatShellGetContext(mat,&ctx);
262: if (type == NORM_FROBENIUS) {
263: *norm = 1.0;
264: } else if (type == NORM_1 || type == NORM_INFINITY) {
265: *norm = 1.0;
266: } else SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"No two norm");
267: return 0;
268: }