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: }