Actual source code: normmh.c
1: #include <petsc/private/matimpl.h>
3: typedef struct {
4: Mat A;
5: Mat D; /* local submatrix for diagonal part */
6: Vec w, left, right, leftwork, rightwork;
7: PetscScalar scale;
8: } Mat_NormalHermitian;
10: static PetscErrorCode MatScale_NormalHermitian(Mat inA, PetscScalar scale)
11: {
12: Mat_NormalHermitian *a = (Mat_NormalHermitian *)inA->data;
14: PetscFunctionBegin;
15: a->scale *= scale;
16: PetscFunctionReturn(PETSC_SUCCESS);
17: }
19: static PetscErrorCode MatDiagonalScale_NormalHermitian(Mat inA, Vec left, Vec right)
20: {
21: Mat_NormalHermitian *a = (Mat_NormalHermitian *)inA->data;
23: PetscFunctionBegin;
24: if (left) {
25: if (!a->left) {
26: PetscCall(VecDuplicate(left, &a->left));
27: PetscCall(VecCopy(left, a->left));
28: } else {
29: PetscCall(VecPointwiseMult(a->left, left, a->left));
30: }
31: }
32: if (right) {
33: if (!a->right) {
34: PetscCall(VecDuplicate(right, &a->right));
35: PetscCall(VecCopy(right, a->right));
36: } else {
37: PetscCall(VecPointwiseMult(a->right, right, a->right));
38: }
39: }
40: PetscFunctionReturn(PETSC_SUCCESS);
41: }
43: static PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
44: {
45: Mat_NormalHermitian *a = (Mat_NormalHermitian *)mat->data;
46: Mat B = a->A, *suba;
47: IS *row;
48: PetscInt M;
50: PetscFunctionBegin;
51: PetscCheck(!a->left && !a->right && irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
52: if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
53: PetscCall(MatGetSize(B, &M, NULL));
54: PetscCall(PetscMalloc1(n, &row));
55: PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
56: PetscCall(ISSetIdentity(row[0]));
57: for (M = 1; M < n; ++M) row[M] = row[0];
58: PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
59: for (M = 0; M < n; ++M) {
60: PetscCall(MatCreateNormalHermitian(suba[M], *submat + M));
61: ((Mat_NormalHermitian *)(*submat)[M]->data)->scale = a->scale;
62: }
63: PetscCall(ISDestroy(&row[0]));
64: PetscCall(PetscFree(row));
65: PetscCall(MatDestroySubMatrices(n, &suba));
66: PetscFunctionReturn(PETSC_SUCCESS);
67: }
69: static PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B)
70: {
71: Mat_NormalHermitian *a = (Mat_NormalHermitian *)A->data;
72: Mat C, Aa = a->A;
73: IS row;
75: PetscFunctionBegin;
76: PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
77: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
78: PetscCall(ISSetIdentity(row));
79: PetscCall(MatPermute(Aa, row, colp, &C));
80: PetscCall(ISDestroy(&row));
81: PetscCall(MatCreateNormalHermitian(C, B));
82: PetscCall(MatDestroy(&C));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: static PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
87: {
88: Mat_NormalHermitian *a = (Mat_NormalHermitian *)A->data;
89: Mat C;
91: PetscFunctionBegin;
92: PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented");
93: PetscCall(MatDuplicate(a->A, op, &C));
94: PetscCall(MatCreateNormalHermitian(C, B));
95: PetscCall(MatDestroy(&C));
96: if (op == MAT_COPY_VALUES) ((Mat_NormalHermitian *)(*B)->data)->scale = a->scale;
97: PetscFunctionReturn(PETSC_SUCCESS);
98: }
100: static PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
101: {
102: Mat_NormalHermitian *a = (Mat_NormalHermitian *)A->data, *b = (Mat_NormalHermitian *)B->data;
104: PetscFunctionBegin;
105: PetscCheck(!a->left && !a->right, PetscObjectComm((PetscObject)A), PETSC_ERR_SUP, "Not implemented");
106: PetscCall(MatCopy(a->A, b->A, str));
107: b->scale = a->scale;
108: PetscCall(VecDestroy(&b->left));
109: PetscCall(VecDestroy(&b->right));
110: PetscCall(VecDestroy(&b->leftwork));
111: PetscCall(VecDestroy(&b->rightwork));
112: PetscFunctionReturn(PETSC_SUCCESS);
113: }
115: static PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
116: {
117: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
118: Vec in;
120: PetscFunctionBegin;
121: in = x;
122: if (Na->right) {
123: if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork));
124: PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in));
125: in = Na->rightwork;
126: }
127: PetscCall(MatMult(Na->A, in, Na->w));
128: PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
129: if (Na->left) PetscCall(VecPointwiseMult(y, Na->left, y));
130: PetscCall(VecScale(y, Na->scale));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: static PetscErrorCode MatMultHermitianAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
135: {
136: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
137: Vec in;
139: PetscFunctionBegin;
140: in = v1;
141: if (Na->right) {
142: if (!Na->rightwork) PetscCall(VecDuplicate(Na->right, &Na->rightwork));
143: PetscCall(VecPointwiseMult(Na->rightwork, Na->right, in));
144: in = Na->rightwork;
145: }
146: PetscCall(MatMult(Na->A, in, Na->w));
147: PetscCall(VecScale(Na->w, Na->scale));
148: if (Na->left) {
149: PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3));
150: PetscCall(VecPointwiseMult(v3, Na->left, v3));
151: PetscCall(VecAXPY(v3, 1.0, v2));
152: } else {
153: PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3));
154: }
155: PetscFunctionReturn(PETSC_SUCCESS);
156: }
158: static PetscErrorCode MatMultHermitianTranspose_Normal(Mat N, Vec x, Vec y)
159: {
160: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
161: Vec in;
163: PetscFunctionBegin;
164: in = x;
165: if (Na->left) {
166: if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork));
167: PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in));
168: in = Na->leftwork;
169: }
170: PetscCall(MatMult(Na->A, in, Na->w));
171: PetscCall(MatMultHermitianTranspose(Na->A, Na->w, y));
172: if (Na->right) PetscCall(VecPointwiseMult(y, Na->right, y));
173: PetscCall(VecScale(y, Na->scale));
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: static PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
178: {
179: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
180: Vec in;
182: PetscFunctionBegin;
183: in = v1;
184: if (Na->left) {
185: if (!Na->leftwork) PetscCall(VecDuplicate(Na->left, &Na->leftwork));
186: PetscCall(VecPointwiseMult(Na->leftwork, Na->left, in));
187: in = Na->leftwork;
188: }
189: PetscCall(MatMult(Na->A, in, Na->w));
190: PetscCall(VecScale(Na->w, Na->scale));
191: if (Na->right) {
192: PetscCall(MatMultHermitianTranspose(Na->A, Na->w, v3));
193: PetscCall(VecPointwiseMult(v3, Na->right, v3));
194: PetscCall(VecAXPY(v3, 1.0, v2));
195: } else {
196: PetscCall(MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3));
197: }
198: PetscFunctionReturn(PETSC_SUCCESS);
199: }
201: static PetscErrorCode MatDestroy_NormalHermitian(Mat N)
202: {
203: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
205: PetscFunctionBegin;
206: PetscCall(MatDestroy(&Na->A));
207: PetscCall(MatDestroy(&Na->D));
208: PetscCall(VecDestroy(&Na->w));
209: PetscCall(VecDestroy(&Na->left));
210: PetscCall(VecDestroy(&Na->right));
211: PetscCall(VecDestroy(&Na->leftwork));
212: PetscCall(VecDestroy(&Na->rightwork));
213: PetscCall(PetscFree(N->data));
214: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMatHermitian_C", NULL));
215: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL));
216: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL));
217: #if defined(PETSC_HAVE_HYPRE)
218: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_hypre_C", NULL));
219: #endif
220: PetscFunctionReturn(PETSC_SUCCESS);
221: }
223: /*
224: Slow, nonscalable version
225: */
226: static PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
227: {
228: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
229: Mat A = Na->A;
230: PetscInt i, j, rstart, rend, nnz;
231: const PetscInt *cols;
232: PetscScalar *diag, *work, *values;
233: const PetscScalar *mvalues;
235: PetscFunctionBegin;
236: PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
237: PetscCall(PetscArrayzero(work, A->cmap->N));
238: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
239: for (i = rstart; i < rend; i++) {
240: PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
241: for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
242: PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
243: }
244: PetscCall(MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
245: rstart = N->cmap->rstart;
246: rend = N->cmap->rend;
247: PetscCall(VecGetArray(v, &values));
248: PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
249: PetscCall(VecRestoreArray(v, &values));
250: PetscCall(PetscFree2(diag, work));
251: PetscCall(VecScale(v, Na->scale));
252: PetscFunctionReturn(PETSC_SUCCESS);
253: }
255: static PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
256: {
257: Mat_NormalHermitian *Na = (Mat_NormalHermitian *)N->data;
258: Mat M, A = Na->A;
260: PetscFunctionBegin;
261: PetscCall(MatGetDiagonalBlock(A, &M));
262: PetscCall(MatCreateNormalHermitian(M, &Na->D));
263: *D = Na->D;
264: PetscFunctionReturn(PETSC_SUCCESS);
265: }
267: static PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A, Mat *M)
268: {
269: Mat_NormalHermitian *Aa = (Mat_NormalHermitian *)A->data;
271: PetscFunctionBegin;
272: *M = Aa->A;
273: PetscFunctionReturn(PETSC_SUCCESS);
274: }
276: /*@
277: MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`
279: Logically Collective
281: Input Parameter:
282: . A - the `MATNORMALHERMITIAN` matrix
284: Output Parameter:
285: . M - the matrix object stored inside A
287: Level: intermediate
289: .seealso: [](ch_matrices), `Mat`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
290: @*/
291: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
292: {
293: PetscFunctionBegin;
296: PetscAssertPointer(M, 2);
297: PetscUseMethod(A, "MatNormalGetMatHermitian_C", (Mat, Mat *), (A, M));
298: PetscFunctionReturn(PETSC_SUCCESS);
299: }
301: static PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
302: {
303: Mat_NormalHermitian *Aa = (Mat_NormalHermitian *)A->data;
304: Mat B, conjugate;
305: PetscInt m, n, M, N;
307: PetscFunctionBegin;
308: PetscCall(MatGetSize(A, &M, &N));
309: PetscCall(MatGetLocalSize(A, &m, &n));
310: if (reuse == MAT_REUSE_MATRIX) {
311: B = *newmat;
312: PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
313: } else {
314: PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
315: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
316: PetscCall(MatProductSetFromOptions(B));
317: PetscCall(MatProductSymbolic(B));
318: PetscCall(MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE));
319: }
320: if (PetscDefined(USE_COMPLEX)) {
321: PetscCall(MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate));
322: PetscCall(MatConjugate(conjugate));
323: PetscCall(MatProductReplaceMats(conjugate, Aa->A, NULL, B));
324: }
325: PetscCall(MatProductNumeric(B));
326: if (PetscDefined(USE_COMPLEX)) PetscCall(MatDestroy(&conjugate));
327: if (reuse == MAT_INPLACE_MATRIX) {
328: PetscCall(MatHeaderReplace(A, &B));
329: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
330: PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
331: PetscFunctionReturn(PETSC_SUCCESS);
332: }
334: #if defined(PETSC_HAVE_HYPRE)
335: static PetscErrorCode MatConvert_NormalHermitian_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
336: {
337: PetscFunctionBegin;
338: if (reuse == MAT_INITIAL_MATRIX) {
339: PetscCall(MatConvert(A, MATAIJ, reuse, B));
340: PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
341: } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
342: PetscFunctionReturn(PETSC_SUCCESS);
343: }
344: #endif
346: /*MC
347: MATNORMALHERMITIAN - a matrix that behaves like (A*)'*A for `MatMult()` while only containing A
349: Level: intermediate
351: .seealso: [](ch_matrices), `Mat`, `MatCreateNormalHermitian()`, `MatMult()`, `MatNormalHermitianGetMat()`, `MATNORMAL`, `MatCreateNormal()`
352: M*/
354: /*@
355: MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A.
357: Collective
359: Input Parameter:
360: . A - the (possibly rectangular complex) matrix
362: Output Parameter:
363: . N - the matrix that represents (A*)'*A
365: Level: intermediate
367: Note:
368: The product (A*)'*A is NOT actually formed! Rather the new matrix
369: object performs the matrix-vector product, `MatMult()`, by first multiplying by
370: A and then (A*)'
372: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
373: @*/
374: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
375: {
376: PetscInt m, n;
377: Mat_NormalHermitian *Na;
378: VecType vtype;
380: PetscFunctionBegin;
381: PetscCall(MatGetLocalSize(A, &m, &n));
382: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
383: PetscCall(MatSetSizes(*N, n, n, PETSC_DECIDE, PETSC_DECIDE));
384: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN));
385: PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
386: PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
388: PetscCall(PetscNew(&Na));
389: (*N)->data = (void *)Na;
390: PetscCall(PetscObjectReference((PetscObject)A));
391: Na->A = A;
392: Na->scale = 1.0;
394: PetscCall(MatCreateVecs(A, NULL, &Na->w));
396: (*N)->ops->destroy = MatDestroy_NormalHermitian;
397: (*N)->ops->mult = MatMult_NormalHermitian;
398: (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal;
399: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
400: (*N)->ops->multadd = MatMultHermitianAdd_Normal;
401: (*N)->ops->getdiagonal = MatGetDiagonal_NormalHermitian;
402: (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_NormalHermitian;
403: (*N)->ops->scale = MatScale_NormalHermitian;
404: (*N)->ops->diagonalscale = MatDiagonalScale_NormalHermitian;
405: (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
406: (*N)->ops->permute = MatPermute_NormalHermitian;
407: (*N)->ops->duplicate = MatDuplicate_NormalHermitian;
408: (*N)->ops->copy = MatCopy_NormalHermitian;
409: (*N)->assembled = PETSC_TRUE;
410: (*N)->preallocated = PETSC_TRUE;
412: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMatHermitian_C", MatNormalGetMat_NormalHermitian));
413: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ));
414: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ));
415: #if defined(PETSC_HAVE_HYPRE)
416: PetscCall(PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_hypre_C", MatConvert_NormalHermitian_HYPRE));
417: #endif
418: PetscCall(MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE));
419: PetscCall(MatGetVecType(A, &vtype));
420: PetscCall(MatSetVecType(*N, vtype));
421: #if defined(PETSC_HAVE_DEVICE)
422: PetscCall(MatBindToCPU(*N, A->boundtocpu));
423: #endif
424: PetscFunctionReturn(PETSC_SUCCESS);
425: }