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: PetscCall(MatShellGetContext(N, &Na));
166: A = Na->A;
167: PetscCall(MatGetDiagonalBlock(A, &M));
168: PetscCall(MatCreateNormalHermitian(M, &Na->D));
169: *D = Na->D;
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M)
174: {
175: Mat_NormalHermitian *Aa;
177: PetscFunctionBegin;
178: PetscCall(MatShellGetContext(A, &Aa));
179: *M = Aa->A;
180: PetscFunctionReturn(PETSC_SUCCESS);
181: }
183: /*@
184: MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`
186: Logically Collective
188: Input Parameter:
189: . A - the `MATNORMALHERMITIAN` matrix
191: Output Parameter:
192: . M - the matrix object stored inside `A`
194: Level: intermediate
196: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
197: @*/
198: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
199: {
200: PetscFunctionBegin;
203: PetscAssertPointer(M, 2);
204: PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M));
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
209: {
210: Mat_NormalHermitian *Aa;
211: Mat B, conjugate;
212: PetscInt m, n, M, N;
214: PetscFunctionBegin;
215: PetscCall(MatShellGetContext(A, &Aa));
216: PetscCall(MatGetSize(A, &M, &N));
217: PetscCall(MatGetLocalSize(A, &m, &n));
218: if (reuse == MAT_REUSE_MATRIX) {
219: B = *newmat;
220: PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
221: } else {
222: PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
223: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
224: PetscCall(MatProductSetFromOptions(B));
225: PetscCall(MatProductSymbolic(B));
226: PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
227: }
228: if (PetscDefined(USE_COMPLEX)) {
229: PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
230: PetscCall(MatConjugate(conjugate));
231: PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
232: }
233: PetscCall(MatProductNumeric(B));
234: if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
235: if (reuse == MAT_INPLACE_MATRIX) {
236: PetscCall(MatHeaderReplace(A, &B));
237: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
238: PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: #if defined(PETSC_HAVE_HYPRE)
243: static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
244: {
245: PetscFunctionBegin;
246: if (reuse == MAT_INITIAL_MATRIX) {
247: PetscCall(MatConvert(A, MATAIJ, reuse, B));
248: PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
249: } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
250: PetscFunctionReturn(PETSC_SUCCESS);
251: }
252: #endif
254: /*MC
255: MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A
257: Level: intermediate
259: Developer Notes:
260: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
262: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
264: .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
265: M*/
267: /*@
268: MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A.
270: Collective
272: Input Parameter:
273: . A - the (possibly rectangular complex) matrix
275: Output Parameter:
276: . N - the matrix that represents (A*)'*A
278: Level: intermediate
280: Note:
281: The product (A*)'*A is NOT actually formed! Rather the new matrix
282: object performs the matrix-vector product, `MatMult()`, by first multiplying by
283: A and then (A*)'
285: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
286: @*/
287: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
288: {
289: Mat_NormalHermitian *Na;
290: VecType vtype;
292: PetscFunctionBegin;
293: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
294: PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
295: PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
296: PetscCall(MatSetType(*N, MATSHELL));
297: PetscCall(PetscNew(&Na));
298: PetscCall(MatShellSetContext(*N, Na));
299: PetscCall(PetscObjectReference((PetscObject)A));
300: Na->A = A;
301: PetscCall(MatCreateVecs(A, NULL, &Na->w));
303: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
304: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_NormalHermitian));
305: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_NormalHermitian));
306: PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian));
307: #if !defined(PETSC_USE_COMPLEX)
308: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian));
309: #endif
310: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_NormalHermitian));
311: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_NormalHermitian));
312: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (void (*)(void))MatGetDiagonalBlock_NormalHermitian));
313: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_NormalHermitian));
314: (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
315: (*N)->ops->permute = MatPermute_NormalHermitian;
317: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalHermitianGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
318: #if !defined(PETSC_USE_COMPLEX)
319: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
320: #endif
321: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ));
322: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ));
323: #if defined(PETSC_HAVE_HYPRE)
324: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE));
325: #endif
326: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
327: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
328: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
329: PetscCall(MatSetOption(*N, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
330: PetscCall(MatGetVecType(A, &vtype));
331: PetscCall(MatSetVecType(*N, vtype));
332: #if defined(PETSC_HAVE_DEVICE)
333: PetscCall(MatBindToCPU(*N, A->boundtocpu));
334: #endif
335: PetscCall(MatSetUp(*N));
336: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN));
337: PetscFunctionReturn(PETSC_SUCCESS);
338: }