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: }