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
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: Level: developer
20: Note:
21: The caller is responsible for destroying the input objects after matrix J has been destroyed.
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: return 0;
62: }
64: PetscErrorCode MatSMFResetRowColumn(Mat mat, IS Rows, IS Cols)
65: {
66: MatSubMatFreeCtx ctx;
68: MatShellGetContext(mat, &ctx);
69: ISDestroy(&ctx->Rows);
70: ISDestroy(&ctx->Cols);
71: PetscObjectReference((PetscObject)Rows);
72: PetscObjectReference((PetscObject)Cols);
73: ctx->Cols = Cols;
74: ctx->Rows = Rows;
75: return 0;
76: }
78: PetscErrorCode MatMult_SMF(Mat mat, Vec a, Vec y)
79: {
80: MatSubMatFreeCtx ctx;
82: MatShellGetContext(mat, &ctx);
83: VecCopy(a, ctx->VR);
84: VecISSet(ctx->VR, ctx->Cols, 0.0);
85: MatMult(ctx->A, ctx->VR, y);
86: VecISSet(y, ctx->Rows, 0.0);
87: return 0;
88: }
90: PetscErrorCode MatMultTranspose_SMF(Mat mat, Vec a, Vec y)
91: {
92: MatSubMatFreeCtx ctx;
94: MatShellGetContext(mat, &ctx);
95: VecCopy(a, ctx->VC);
96: VecISSet(ctx->VC, ctx->Rows, 0.0);
97: MatMultTranspose(ctx->A, ctx->VC, y);
98: VecISSet(y, ctx->Cols, 0.0);
99: return 0;
100: }
102: PetscErrorCode MatDiagonalSet_SMF(Mat M, Vec D, InsertMode is)
103: {
104: MatSubMatFreeCtx ctx;
106: MatShellGetContext(M, &ctx);
107: MatDiagonalSet(ctx->A, D, is);
108: return 0;
109: }
111: PetscErrorCode MatDestroy_SMF(Mat mat)
112: {
113: MatSubMatFreeCtx ctx;
115: MatShellGetContext(mat, &ctx);
116: MatDestroy(&ctx->A);
117: ISDestroy(&ctx->Rows);
118: ISDestroy(&ctx->Cols);
119: VecDestroy(&ctx->VC);
120: PetscFree(ctx);
121: return 0;
122: }
124: PetscErrorCode MatView_SMF(Mat mat, PetscViewer viewer)
125: {
126: MatSubMatFreeCtx ctx;
128: MatShellGetContext(mat, &ctx);
129: MatView(ctx->A, viewer);
130: return 0;
131: }
133: PetscErrorCode MatShift_SMF(Mat Y, PetscReal a)
134: {
135: MatSubMatFreeCtx ctx;
137: MatShellGetContext(Y, &ctx);
138: MatShift(ctx->A, a);
139: return 0;
140: }
142: PetscErrorCode MatDuplicate_SMF(Mat mat, MatDuplicateOption op, Mat *M)
143: {
144: MatSubMatFreeCtx ctx;
146: MatShellGetContext(mat, &ctx);
147: MatCreateSubMatrixFree(ctx->A, ctx->Rows, ctx->Cols, M);
148: return 0;
149: }
151: PetscErrorCode MatEqual_SMF(Mat A, Mat B, PetscBool *flg)
152: {
153: MatSubMatFreeCtx ctx1, ctx2;
154: PetscBool flg1, flg2, flg3;
156: MatShellGetContext(A, &ctx1);
157: MatShellGetContext(B, &ctx2);
158: ISEqual(ctx1->Rows, ctx2->Rows, &flg2);
159: ISEqual(ctx1->Cols, ctx2->Cols, &flg3);
160: if (flg2 == PETSC_FALSE || flg3 == PETSC_FALSE) {
161: *flg = PETSC_FALSE;
162: } else {
163: MatEqual(ctx1->A, ctx2->A, &flg1);
164: if (flg1 == PETSC_FALSE) {
165: *flg = PETSC_FALSE;
166: } else {
167: *flg = PETSC_TRUE;
168: }
169: }
170: return 0;
171: }
173: PetscErrorCode MatScale_SMF(Mat mat, PetscReal a)
174: {
175: MatSubMatFreeCtx ctx;
177: MatShellGetContext(mat, &ctx);
178: MatScale(ctx->A, a);
179: return 0;
180: }
182: PetscErrorCode MatTranspose_SMF(Mat mat, Mat *B)
183: {
184: SETERRQ(PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "No support for transpose for MatCreateSubMatrixFree() matrix");
185: return 0;
186: }
188: PetscErrorCode MatGetDiagonal_SMF(Mat mat, Vec v)
189: {
190: MatSubMatFreeCtx ctx;
192: MatShellGetContext(mat, &ctx);
193: MatGetDiagonal(ctx->A, v);
194: return 0;
195: }
197: PetscErrorCode MatGetRowMax_SMF(Mat M, Vec D)
198: {
199: MatSubMatFreeCtx ctx;
201: MatShellGetContext(M, &ctx);
202: MatGetRowMax(ctx->A, D, NULL);
203: return 0;
204: }
206: PetscErrorCode MatCreateSubMatrices_SMF(Mat A, PetscInt n, IS *irow, IS *icol, MatReuse scall, Mat **B)
207: {
208: PetscInt i;
210: if (scall == MAT_INITIAL_MATRIX) PetscCalloc1(n + 1, B);
212: for (i = 0; i < n; i++) MatCreateSubMatrix_SMF(A, irow[i], icol[i], scall, &(*B)[i]);
213: return 0;
214: }
216: PetscErrorCode MatCreateSubMatrix_SMF(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
217: {
218: MatSubMatFreeCtx ctx;
220: MatShellGetContext(mat, &ctx);
221: if (newmat) MatDestroy(&*newmat);
222: MatCreateSubMatrixFree(ctx->A, isrow, iscol, newmat);
223: return 0;
224: }
226: PetscErrorCode MatGetRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
227: {
228: MatSubMatFreeCtx ctx;
230: MatShellGetContext(mat, &ctx);
231: MatGetRow(ctx->A, row, ncols, cols, vals);
232: return 0;
233: }
235: PetscErrorCode MatRestoreRow_SMF(Mat mat, PetscInt row, PetscInt *ncols, const PetscInt **cols, const PetscScalar **vals)
236: {
237: MatSubMatFreeCtx ctx;
239: MatShellGetContext(mat, &ctx);
240: MatRestoreRow(ctx->A, row, ncols, cols, vals);
241: return 0;
242: }
244: PetscErrorCode MatGetColumnVector_SMF(Mat mat, Vec Y, PetscInt col)
245: {
246: MatSubMatFreeCtx ctx;
248: MatShellGetContext(mat, &ctx);
249: MatGetColumnVector(ctx->A, Y, col);
250: return 0;
251: }
253: PetscErrorCode MatNorm_SMF(Mat mat, NormType type, PetscReal *norm)
254: {
255: MatSubMatFreeCtx ctx;
257: MatShellGetContext(mat, &ctx);
258: if (type == NORM_FROBENIUS) {
259: *norm = 1.0;
260: } else if (type == NORM_1 || type == NORM_INFINITY) {
261: *norm = 1.0;
262: } else SETERRQ(PETSC_COMM_SELF, PETSC_ERR_SUP, "No two norm");
263: return 0;
264: }