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