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: PetscScalar shift, scale;
15: PetscInt M;
17: PetscFunctionBegin;
18: PetscCheck(irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
19: PetscCall(MatShellGetScalingShifts(mat, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
20: PetscCall(MatShellGetContext(mat, &a));
21: B = a->A;
22: if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
23: PetscCall(MatGetSize(B, &M, NULL));
24: PetscCall(PetscMalloc1(n, &row));
25: PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
26: PetscCall(ISSetIdentity(row[0]));
27: for (M = 1; M < n; ++M) row[M] = row[0];
28: PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
29: for (M = 0; M < n; ++M) {
30: PetscCall(MatCreateNormalHermitian(suba[M], *submat + M));
31: PetscCall(MatShift((*submat)[M], shift));
32: PetscCall(MatScale((*submat)[M], scale));
33: }
34: PetscCall(ISDestroy(&row[0]));
35: PetscCall(PetscFree(row));
36: PetscCall(MatDestroySubMatrices(n, &suba));
37: PetscFunctionReturn(PETSC_SUCCESS);
38: }
40: static PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B)
41: {
42: Mat_NormalHermitian *a;
43: Mat C, Aa;
44: IS row;
45: PetscScalar shift, scale;
47: PetscFunctionBegin;
48: PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
49: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
50: PetscCall(MatShellGetContext(A, &a));
51: Aa = a->A;
52: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
53: PetscCall(ISSetIdentity(row));
54: PetscCall(MatPermute(Aa, row, colp, &C));
55: PetscCall(ISDestroy(&row));
56: PetscCall(MatCreateNormalHermitian(C, B));
57: PetscCall(MatDestroy(&C));
58: PetscCall(MatShift(*B, shift));
59: PetscCall(MatScale(*B, scale));
60: PetscFunctionReturn(PETSC_SUCCESS);
61: }
63: static PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
64: {
65: Mat_NormalHermitian *a;
66: Mat C;
68: PetscFunctionBegin;
69: PetscCall(MatShellGetContext(A, &a));
70: PetscCall(MatDuplicate(a->A, op, &C));
71: PetscCall(MatCreateNormalHermitian(C, B));
72: PetscCall(MatDestroy(&C));
73: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
74: PetscFunctionReturn(PETSC_SUCCESS);
75: }
77: static PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
78: {
79: Mat_NormalHermitian *a, *b;
81: PetscFunctionBegin;
82: PetscCall(MatShellGetContext(A, &a));
83: PetscCall(MatShellGetContext(B, &b));
84: PetscCall(MatCopy(a->A, b->A, str));
85: PetscFunctionReturn(PETSC_SUCCESS);
86: }
88: static PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
89: {
90: Mat_NormalHermitian *Na;
92: PetscFunctionBegin;
93: PetscCall(MatShellGetContext(N, &Na));
94: PetscCall(MatMult(Na->A, x, Na->w));
95: PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
96: PetscFunctionReturn(PETSC_SUCCESS);
97: }
99: static PetscErrorCode MatDestroy_NormalHermitian(Mat N)
100: {
101: Mat_NormalHermitian *Na;
103: PetscFunctionBegin;
104: PetscCall(MatShellGetContext(N, &Na));
105: PetscCall(MatDestroy(&Na->A));
106: PetscCall(MatDestroy(&Na->D));
107: PetscCall(VecDestroy(&Na->w));
108: PetscCall(PetscFree(Na));
109: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalHermitianGetMat_C", NULL));
110: #if !defined(PETSC_USE_COMPLEX)
111: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL));
112: #endif
113: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL));
114: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL));
115: #if defined(PETSC_HAVE_HYPRE)
116: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_hypre_C", NULL));
117: #endif
118: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
119: PetscFunctionReturn(PETSC_SUCCESS);
120: }
122: /*
123: Slow, nonscalable version
124: */
125: static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
126: {
127: Mat_NormalHermitian *Na;
128: Mat A;
129: PetscInt i, j, rstart, rend, nnz;
130: const PetscInt *cols;
131: PetscScalar *diag, *work, *values;
132: const PetscScalar *mvalues;
134: PetscFunctionBegin;
135: PetscCall(MatShellGetContext(N, &Na));
136: A = Na->A;
137: PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
138: PetscCall(PetscArrayzero(work, A->cmap->N));
139: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
140: for (i = rstart; i < rend; i++) {
141: PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
142: for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
143: PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
144: }
145: PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
146: rstart = N->cmap->rstart;
147: rend = N->cmap->rend;
148: PetscCall(VecGetArray(v, &values));
149: PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
150: PetscCall(VecRestoreArray(v, &values));
151: PetscCall(PetscFree2(diag, work));
152: PetscFunctionReturn(PETSC_SUCCESS);
153: }
155: static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
156: {
157: Mat_NormalHermitian *Na;
158: Mat M, A;
160: PetscFunctionBegin;
161: 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
162: 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
163: 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
164: 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
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_COPY, (void (*)(void))MatCopy_NormalHermitian));
313: (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_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: }