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;
133: PetscMPIInt iN;
135: PetscFunctionBegin;
136: PetscCall(MatShellGetContext(N, &Na));
137: A = Na->A;
138: PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
139: PetscCall(PetscArrayzero(work, A->cmap->N));
140: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
141: for (i = rstart; i < rend; i++) {
142: PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
143: for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
144: PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
145: }
146: PetscCall(PetscMPIIntCast(A->cmap->N, &iN));
147: PetscCallMPI(MPIU_Allreduce(work, diag, iN, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
148: rstart = N->cmap->rstart;
149: rend = N->cmap->rend;
150: PetscCall(VecGetArray(v, &values));
151: PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
152: PetscCall(VecRestoreArray(v, &values));
153: PetscCall(PetscFree2(diag, work));
154: PetscFunctionReturn(PETSC_SUCCESS);
155: }
157: static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
158: {
159: Mat_NormalHermitian *Na;
160: Mat M, A;
162: PetscFunctionBegin;
163: PetscCall(MatShellGetContext(N, &Na));
164: A = Na->A;
165: PetscCall(MatGetDiagonalBlock(A, &M));
166: PetscCall(MatCreateNormalHermitian(M, &Na->D));
167: *D = Na->D;
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: static PetscErrorCode MatNormalHermitianGetMat_NormalHermitian(Mat A, Mat *M)
172: {
173: Mat_NormalHermitian *Aa;
175: PetscFunctionBegin;
176: PetscCall(MatShellGetContext(A, &Aa));
177: *M = Aa->A;
178: PetscFunctionReturn(PETSC_SUCCESS);
179: }
181: /*@
182: MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`
184: Logically Collective
186: Input Parameter:
187: . A - the `MATNORMALHERMITIAN` matrix
189: Output Parameter:
190: . M - the matrix object stored inside `A`
192: Level: intermediate
194: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
195: @*/
196: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
197: {
198: PetscFunctionBegin;
201: PetscAssertPointer(M, 2);
202: PetscUseMethod(A, "MatNormalHermitianGetMat_C", (Mat, Mat *), (A, M));
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
207: {
208: Mat_NormalHermitian *Aa;
209: Mat B, conjugate;
210: PetscInt m, n, M, N;
212: PetscFunctionBegin;
213: PetscCall(MatShellGetContext(A, &Aa));
214: PetscCall(MatGetSize(A, &M, &N));
215: PetscCall(MatGetLocalSize(A, &m, &n));
216: if (reuse == MAT_REUSE_MATRIX) {
217: B = *newmat;
218: PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
219: } else {
220: PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
221: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
222: PetscCall(MatProductSetFromOptions(B));
223: PetscCall(MatProductSymbolic(B));
224: PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
225: }
226: if (PetscDefined(USE_COMPLEX)) {
227: PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
228: PetscCall(MatConjugate(conjugate));
229: PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
230: }
231: PetscCall(MatProductNumeric(B));
232: if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
233: if (reuse == MAT_INPLACE_MATRIX) {
234: PetscCall(MatHeaderReplace(A, &B));
235: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
236: PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
237: PetscFunctionReturn(PETSC_SUCCESS);
238: }
240: #if defined(PETSC_HAVE_HYPRE)
241: static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
242: {
243: PetscFunctionBegin;
244: if (reuse == MAT_INITIAL_MATRIX) {
245: PetscCall(MatConvert(A, MATAIJ, reuse, B));
246: PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
247: } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
248: PetscFunctionReturn(PETSC_SUCCESS);
249: }
250: #endif
252: /*MC
253: MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A
255: Level: intermediate
257: Developer Notes:
258: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
260: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
262: .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
263: M*/
265: /*@
266: MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A.
268: Collective
270: Input Parameter:
271: . A - the (possibly rectangular complex) matrix
273: Output Parameter:
274: . N - the matrix that represents (A*)'*A
276: Level: intermediate
278: Note:
279: The product (A*)'*A is NOT actually formed! Rather the new matrix
280: object performs the matrix-vector product, `MatMult()`, by first multiplying by
281: A and then (A*)'
283: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
284: @*/
285: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
286: {
287: Mat_NormalHermitian *Na;
288: VecType vtype;
290: PetscFunctionBegin;
291: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
292: PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
293: PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
294: PetscCall(MatSetType(*N, MATSHELL));
295: PetscCall(PetscNew(&Na));
296: PetscCall(MatShellSetContext(*N, Na));
297: PetscCall(PetscObjectReference((PetscObject)A));
298: Na->A = A;
299: PetscCall(MatCreateVecs(A, NULL, &Na->w));
301: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
302: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_NormalHermitian));
303: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_NormalHermitian));
304: PetscCall(MatShellSetOperation(*N, MATOP_MULT_HERMITIAN_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian));
305: #if !defined(PETSC_USE_COMPLEX)
306: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_NormalHermitian));
307: #endif
308: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_NormalHermitian));
309: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_NormalHermitian));
310: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (void (*)(void))MatGetDiagonalBlock_NormalHermitian));
311: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_NormalHermitian));
312: (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
313: (*N)->ops->permute = MatPermute_NormalHermitian;
315: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalHermitianGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
316: #if !defined(PETSC_USE_COMPLEX)
317: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalHermitianGetMat_NormalHermitian));
318: #endif
319: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ));
320: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ));
321: #if defined(PETSC_HAVE_HYPRE)
322: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE));
323: #endif
324: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
325: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
326: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
327: PetscCall(MatSetOption(*N, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
328: PetscCall(MatGetVecType(A, &vtype));
329: PetscCall(MatSetVecType(*N, vtype));
330: #if defined(PETSC_HAVE_DEVICE)
331: PetscCall(MatBindToCPU(*N, A->boundtocpu));
332: #endif
333: PetscCall(MatSetUp(*N));
334: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN));
335: PetscFunctionReturn(PETSC_SUCCESS);
336: }