Actual source code: normmh.c

  1: #include <../src/mat/impls/shell/shell.h>

  3: typedef struct {
  4:   Mat A;
  5:   Mat D; /* local submatrix for diagonal part */
  6:   Vec w;
  7: } Mat_NormalHermitian;

  9: static PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
 10: {
 11:   Mat_NormalHermitian *a;
 12:   Mat                  B, *suba;
 13:   IS                  *row;
 14:   PetscInt             M;

 16:   PetscFunctionBegin;
 17:   PetscCheck(!((Mat_Shell *)mat->data)->zrows && !((Mat_Shell *)mat->data)->zcols, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatZeroRows()/MatZeroRowsColumns() after the SubMatrices creation
 18:   PetscCheck(!((Mat_Shell *)mat->data)->axpy, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatAXPY() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatAXPY() after the SubMatrices creation
 19:   PetscCheck(!((Mat_Shell *)mat->data)->left && !((Mat_Shell *)mat->data)->right, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalScale() after the SubMatrices creation with a SubVector
 20:   PetscCheck(!((Mat_Shell *)mat->data)->dshift, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalSet() after the SubMatrices creation with a SubVector
 21:   PetscCheck(irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
 22:   PetscCall(MatShellGetContext(mat, &a));
 23:   B = a->A;
 24:   if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
 25:   PetscCall(MatGetSize(B, &M, NULL));
 26:   PetscCall(PetscMalloc1(n, &row));
 27:   PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
 28:   PetscCall(ISSetIdentity(row[0]));
 29:   for (M = 1; M < n; ++M) row[M] = row[0];
 30:   PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
 31:   for (M = 0; M < n; ++M) {
 32:     PetscCall(MatCreateNormalHermitian(suba[M], *submat + M));
 33:     ((Mat_Shell *)(*submat)[M]->data)->vscale = ((Mat_Shell *)mat->data)->vscale;
 34:     ((Mat_Shell *)(*submat)[M]->data)->vshift = ((Mat_Shell *)mat->data)->vshift;
 35:   }
 36:   PetscCall(ISDestroy(&row[0]));
 37:   PetscCall(PetscFree(row));
 38:   PetscCall(MatDestroySubMatrices(n, &suba));
 39:   PetscFunctionReturn(PETSC_SUCCESS);
 40: }

 42: static PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B)
 43: {
 44:   Mat_NormalHermitian *a;
 45:   Mat                  C, Aa;
 46:   IS                   row;

 48:   PetscFunctionBegin;
 49:   PetscCheck(!((Mat_Shell *)A->data)->zrows && !((Mat_Shell *)A->data)->zcols, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatZeroRows()/MatZeroRowsColumns() after the permutation with a permuted zrows and zcols
 50:   PetscCheck(!((Mat_Shell *)A->data)->axpy, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatAXPY() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatAXPY() after the permutation with a permuted axpy
 51:   PetscCheck(!((Mat_Shell *)A->data)->left && !((Mat_Shell *)A->data)->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalScale() after the permutation with a permuted left and right
 52:   PetscCheck(!((Mat_Shell *)A->data)->dshift, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Cannot call MatCreateSubMatrices() if MatDiagonalSet() has been called on the input Mat"); // TODO FIXME: lift this limitation by calling MatDiagonalSet() after the permutation with a permuted dshift
 53:   PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
 54:   PetscCall(MatShellGetContext(A, &a));
 55:   Aa = a->A;
 56:   PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
 57:   PetscCall(ISSetIdentity(row));
 58:   PetscCall(MatPermute(Aa, row, colp, &C));
 59:   PetscCall(ISDestroy(&row));
 60:   PetscCall(MatCreateNormalHermitian(C, B));
 61:   PetscCall(MatDestroy(&C));
 62:   ((Mat_Shell *)(*B)->data)->vscale = ((Mat_Shell *)A->data)->vscale;
 63:   ((Mat_Shell *)(*B)->data)->vshift = ((Mat_Shell *)A->data)->vshift;
 64:   PetscFunctionReturn(PETSC_SUCCESS);
 65: }

 67: static PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
 68: {
 69:   Mat_NormalHermitian *a;
 70:   Mat                  C;

 72:   PetscFunctionBegin;
 73:   PetscCall(MatShellGetContext(A, &a));
 74:   PetscCall(MatDuplicate(a->A, op, &C));
 75:   PetscCall(MatCreateNormalHermitian(C, B));
 76:   PetscCall(MatDestroy(&C));
 77:   if (op == MAT_COPY_VALUES) PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
 78:   PetscFunctionReturn(PETSC_SUCCESS);
 79: }

 81: static PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
 82: {
 83:   Mat_NormalHermitian *a, *b;

 85:   PetscFunctionBegin;
 86:   PetscCall(MatShellGetContext(A, &a));
 87:   PetscCall(MatShellGetContext(B, &b));
 88:   PetscCall(MatCopy(a->A, b->A, str));
 89:   PetscFunctionReturn(PETSC_SUCCESS);
 90: }

 92: static PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
 93: {
 94:   Mat_NormalHermitian *Na;

 96:   PetscFunctionBegin;
 97:   PetscCall(MatShellGetContext(N, &Na));
 98:   PetscCall(MatMult(Na->A, x, Na->w));
 99:   PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
100:   PetscFunctionReturn(PETSC_SUCCESS);
101: }

103: static PetscErrorCode MatDestroy_NormalHermitian(Mat N)
104: {
105:   Mat_NormalHermitian *Na;

107:   PetscFunctionBegin;
108:   PetscCall(MatShellGetContext(N, &Na));
109:   PetscCall(MatDestroy(&Na->A));
110:   PetscCall(MatDestroy(&Na->D));
111:   PetscCall(VecDestroy(&Na->w));
112:   PetscCall(PetscFree(Na));
113:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalHermitianGetMat_C", NULL));
114: #if !defined(PETSC_USE_COMPLEX)
115:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL));
116: #endif
117:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL));
118:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL));
119: #if defined(PETSC_HAVE_HYPRE)
120:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_hypre_C", NULL));
121: #endif
122:   PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
123:   PetscFunctionReturn(PETSC_SUCCESS);
124: }

126: /*
127:       Slow, nonscalable version
128: */
129: static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
130: {
131:   Mat_NormalHermitian *Na;
132:   Mat                  A;
133:   PetscInt             i, j, rstart, rend, nnz;
134:   const PetscInt      *cols;
135:   PetscScalar         *diag, *work, *values;
136:   const PetscScalar   *mvalues;

138:   PetscFunctionBegin;
139:   PetscCall(MatShellGetContext(N, &Na));
140:   A = Na->A;
141:   PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
142:   PetscCall(PetscArrayzero(work, A->cmap->N));
143:   PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
144:   for (i = rstart; i < rend; i++) {
145:     PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
146:     for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
147:     PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
148:   }
149:   PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
150:   rstart = N->cmap->rstart;
151:   rend   = N->cmap->rend;
152:   PetscCall(VecGetArray(v, &values));
153:   PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
154:   PetscCall(VecRestoreArray(v, &values));
155:   PetscCall(PetscFree2(diag, work));
156:   PetscFunctionReturn(PETSC_SUCCESS);
157: }

159: static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
160: {
161:   Mat_NormalHermitian *Na;
162:   Mat                  M, A;

164:   PetscFunctionBegin;
165:   PetscCheck(!((Mat_Shell *)N->data)->zrows && !((Mat_Shell *)N->data)->zcols, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatZeroRows() or MatZeroRowsColumns() has been called on the input Mat"); // TODO FIXME
166:   PetscCheck(!((Mat_Shell *)N->data)->axpy, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatAXPY() has been called on the input Mat");                                            // TODO FIXME
167:   PetscCheck(!((Mat_Shell *)N->data)->left && !((Mat_Shell *)N->data)->right, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalScale() has been called on the input Mat"); // TODO FIXME
168:   PetscCheck(!((Mat_Shell *)N->data)->dshift, PetscObjectComm((PetscObject)N), PETSC_ERR_SUP, "Cannot call MatGetDiagonalBlock() if MatDiagonalSet() has been called on the input Mat");                                   // TODO FIXME
169:   PetscCall(MatShellGetContext(N, &Na));
170:   A = Na->A;
171:   PetscCall(MatGetDiagonalBlock(A, &M));
172:   PetscCall(MatCreateNormalHermitian(M, &Na->D));
173:   *D = Na->D;
174:   PetscFunctionReturn(PETSC_SUCCESS);
175: }

177: static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M)
178: {
179:   Mat_NormalHermitian *Aa;

181:   PetscFunctionBegin;
182:   PetscCall(MatShellGetContext(A, &Aa));
183:   *M = Aa->A;
184:   PetscFunctionReturn(PETSC_SUCCESS);
185: }

187: /*@
188:   MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`

190:   Logically Collective

192:   Input Parameter:
193: . A - the `MATNORMALHERMITIAN` matrix

195:   Output Parameter:
196: . M - the matrix object stored inside `A`

198:   Level: intermediate

200: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
201: @*/
202: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
203: {
204:   PetscFunctionBegin;
207:   PetscAssertPointer(M, 2);
208:   PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M));
209:   PetscFunctionReturn(PETSC_SUCCESS);
210: }

212: static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
213: {
214:   Mat_NormalHermitian *Aa;
215:   Mat                  B, conjugate;
216:   PetscInt             m, n, M, N;

218:   PetscFunctionBegin;
219:   PetscCall(MatShellGetContext(A, &Aa));
220:   PetscCall(MatGetSize(A, &M, &N));
221:   PetscCall(MatGetLocalSize(A, &m, &n));
222:   if (reuse == MAT_REUSE_MATRIX) {
223:     B = *newmat;
224:     PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
225:   } else {
226:     PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
227:     PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
228:     PetscCall(MatProductSetFromOptions(B));
229:     PetscCall(MatProductSymbolic(B));
230:     PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
231:   }
232:   if (PetscDefined(USE_COMPLEX)) {
233:     PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
234:     PetscCall(MatConjugate(conjugate));
235:     PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
236:   }
237:   PetscCall(MatProductNumeric(B));
238:   if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
239:   if (reuse == MAT_INPLACE_MATRIX) {
240:     PetscCall(MatHeaderReplace(A, &B));
241:   } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
242:   PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
243:   PetscFunctionReturn(PETSC_SUCCESS);
244: }

246: #if defined(PETSC_HAVE_HYPRE)
247: static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
248: {
249:   PetscFunctionBegin;
250:   if (reuse == MAT_INITIAL_MATRIX) {
251:     PetscCall(MatConvert(A, MATAIJ, reuse, B));
252:     PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
253:   } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }
256: #endif

258: /*MC
259:   MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A

261:   Level: intermediate

263:   Developer Notes:
264:   This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code

266:   Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage

268: .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
269: M*/

271: /*@
272:   MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A.

274:   Collective

276:   Input Parameter:
277: . A - the (possibly rectangular complex) matrix

279:   Output Parameter:
280: . N - the matrix that represents (A*)'*A

282:   Level: intermediate

284:   Note:
285:   The product (A*)'*A is NOT actually formed! Rather the new matrix
286:   object performs the matrix-vector product, `MatMult()`, by first multiplying by
287:   A and then (A*)'

289: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
290: @*/
291: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
292: {
293:   Mat_NormalHermitian *Na;
294:   VecType              vtype;

296:   PetscFunctionBegin;
297:   PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
298:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
299:   PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
300:   PetscCall(MatSetType(*N, MATSHELL));
301:   PetscCall(PetscNew(&Na));
302:   PetscCall(MatShellSetContext(*N, Na));
303:   PetscCall(PetscObjectReference((PetscObject)A));
304:   Na->A = A;
305:   PetscCall(MatCreateVecs(A, NULL, &Na->w));

307:   PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
308:   PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_NormalHermitian));
309:   PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_NormalHermitian));
310:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian));
311: #if !defined(PETSC_USE_COMPLEX)
312:   PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian));
313: #endif
314:   PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_NormalHermitian));
315:   PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_NormalHermitian));
316:   PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_NormalHermitian));
317:   (*N)->ops->getdiagonalblock  = MatGetDiagonalBlock_NormalHermitian;
318:   (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
319:   (*N)->ops->permute           = MatPermute_NormalHermitian;

321:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalHermitianGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
322: #if !defined(PETSC_USE_COMPLEX)
323:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
324: #endif
325:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ));
326:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ));
327: #if defined(PETSC_HAVE_HYPRE)
328:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE));
329: #endif
330:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
331:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
332:   PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
333:   PetscCall(MatSetOption(*N, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
334:   PetscCall(MatGetVecType(A, &vtype));
335:   PetscCall(MatSetVecType(*N, vtype));
336: #if defined(PETSC_HAVE_DEVICE)
337:   PetscCall(MatBindToCPU(*N, A->boundtocpu));
338: #endif
339:   PetscCall(MatSetUp(*N));
340:   PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN));
341:   PetscFunctionReturn(PETSC_SUCCESS);
342: }