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