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: a->scale *= scale;
16: return 0;
17: }
19: PetscErrorCode MatDiagonalScale_NormalHermitian(Mat inA, Vec left, Vec right)
20: {
21: Mat_Normal *a = (Mat_Normal *)inA->data;
23: if (left) {
24: if (!a->left) {
25: VecDuplicate(left, &a->left);
26: VecCopy(left, a->left);
27: } else {
28: VecPointwiseMult(a->left, left, a->left);
29: }
30: }
31: if (right) {
32: if (!a->right) {
33: VecDuplicate(right, &a->right);
34: VecCopy(right, a->right);
35: } else {
36: VecPointwiseMult(a->right, right, a->right);
37: }
38: }
39: return 0;
40: }
42: PetscErrorCode MatCreateSubMatrices_NormalHermitian(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
43: {
44: Mat_Normal *a = (Mat_Normal *)mat->data;
45: Mat B = a->A, *suba;
46: IS *row;
47: PetscInt M;
50: if (scall != MAT_REUSE_MATRIX) PetscCalloc1(n, submat);
51: MatGetSize(B, &M, NULL);
52: PetscMalloc1(n, &row);
53: ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]);
54: ISSetIdentity(row[0]);
55: for (M = 1; M < n; ++M) row[M] = row[0];
56: MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba);
57: for (M = 0; M < n; ++M) {
58: MatCreateNormalHermitian(suba[M], *submat + M);
59: ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale;
60: }
61: ISDestroy(&row[0]);
62: PetscFree(row);
63: MatDestroySubMatrices(n, &suba);
64: return 0;
65: }
67: PetscErrorCode MatPermute_NormalHermitian(Mat A, IS rowp, IS colp, Mat *B)
68: {
69: Mat_Normal *a = (Mat_Normal *)A->data;
70: Mat C, Aa = a->A;
71: IS row;
74: ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row);
75: ISSetIdentity(row);
76: MatPermute(Aa, row, colp, &C);
77: ISDestroy(&row);
78: MatCreateNormalHermitian(C, B);
79: MatDestroy(&C);
80: return 0;
81: }
83: PetscErrorCode MatDuplicate_NormalHermitian(Mat A, MatDuplicateOption op, Mat *B)
84: {
85: Mat_Normal *a = (Mat_Normal *)A->data;
86: Mat C;
89: MatDuplicate(a->A, op, &C);
90: MatCreateNormalHermitian(C, B);
91: MatDestroy(&C);
92: if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale;
93: return 0;
94: }
96: PetscErrorCode MatCopy_NormalHermitian(Mat A, Mat B, MatStructure str)
97: {
98: Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data;
101: MatCopy(a->A, b->A, str);
102: b->scale = a->scale;
103: VecDestroy(&b->left);
104: VecDestroy(&b->right);
105: VecDestroy(&b->leftwork);
106: VecDestroy(&b->rightwork);
107: return 0;
108: }
110: PetscErrorCode MatMult_NormalHermitian(Mat N, Vec x, Vec y)
111: {
112: Mat_Normal *Na = (Mat_Normal *)N->data;
113: Vec in;
115: in = x;
116: if (Na->right) {
117: if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
118: VecPointwiseMult(Na->rightwork, Na->right, in);
119: in = Na->rightwork;
120: }
121: MatMult(Na->A, in, Na->w);
122: MatMultHermitianTranspose(Na->A, Na->w, y);
123: if (Na->left) VecPointwiseMult(y, Na->left, y);
124: VecScale(y, Na->scale);
125: return 0;
126: }
128: PetscErrorCode MatMultHermitianAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
129: {
130: Mat_Normal *Na = (Mat_Normal *)N->data;
131: Vec in;
133: in = v1;
134: if (Na->right) {
135: if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
136: VecPointwiseMult(Na->rightwork, Na->right, in);
137: in = Na->rightwork;
138: }
139: MatMult(Na->A, in, Na->w);
140: VecScale(Na->w, Na->scale);
141: if (Na->left) {
142: MatMultHermitianTranspose(Na->A, Na->w, v3);
143: VecPointwiseMult(v3, Na->left, v3);
144: VecAXPY(v3, 1.0, v2);
145: } else {
146: MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3);
147: }
148: return 0;
149: }
151: PetscErrorCode MatMultHermitianTranspose_Normal(Mat N, Vec x, Vec y)
152: {
153: Mat_Normal *Na = (Mat_Normal *)N->data;
154: Vec in;
156: in = x;
157: if (Na->left) {
158: if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
159: VecPointwiseMult(Na->leftwork, Na->left, in);
160: in = Na->leftwork;
161: }
162: MatMult(Na->A, in, Na->w);
163: MatMultHermitianTranspose(Na->A, Na->w, y);
164: if (Na->right) VecPointwiseMult(y, Na->right, y);
165: VecScale(y, Na->scale);
166: return 0;
167: }
169: PetscErrorCode MatMultHermitianTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
170: {
171: Mat_Normal *Na = (Mat_Normal *)N->data;
172: Vec in;
174: in = v1;
175: if (Na->left) {
176: if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
177: VecPointwiseMult(Na->leftwork, Na->left, in);
178: in = Na->leftwork;
179: }
180: MatMult(Na->A, in, Na->w);
181: VecScale(Na->w, Na->scale);
182: if (Na->right) {
183: MatMultHermitianTranspose(Na->A, Na->w, v3);
184: VecPointwiseMult(v3, Na->right, v3);
185: VecAXPY(v3, 1.0, v2);
186: } else {
187: MatMultHermitianTransposeAdd(Na->A, Na->w, v2, v3);
188: }
189: return 0;
190: }
192: PetscErrorCode MatDestroy_NormalHermitian(Mat N)
193: {
194: Mat_Normal *Na = (Mat_Normal *)N->data;
196: MatDestroy(&Na->A);
197: MatDestroy(&Na->D);
198: VecDestroy(&Na->w);
199: VecDestroy(&Na->left);
200: VecDestroy(&Na->right);
201: VecDestroy(&Na->leftwork);
202: VecDestroy(&Na->rightwork);
203: PetscFree(N->data);
204: PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMatHermitian_C", NULL);
205: PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_seqaij_C", NULL);
206: PetscObjectComposeFunction((PetscObject)N, "MatConvert_normalh_mpiaij_C", NULL);
207: return 0;
208: }
210: /*
211: Slow, nonscalable version
212: */
213: PetscErrorCode MatGetDiagonal_NormalHermitian(Mat N, Vec v)
214: {
215: Mat_Normal *Na = (Mat_Normal *)N->data;
216: Mat A = Na->A;
217: PetscInt i, j, rstart, rend, nnz;
218: const PetscInt *cols;
219: PetscScalar *diag, *work, *values;
220: const PetscScalar *mvalues;
222: PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work);
223: PetscArrayzero(work, A->cmap->N);
224: MatGetOwnershipRange(A, &rstart, &rend);
225: for (i = rstart; i < rend; i++) {
226: MatGetRow(A, i, &nnz, &cols, &mvalues);
227: for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * PetscConj(mvalues[j]);
228: MatRestoreRow(A, i, &nnz, &cols, &mvalues);
229: }
230: MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N));
231: rstart = N->cmap->rstart;
232: rend = N->cmap->rend;
233: VecGetArray(v, &values);
234: PetscArraycpy(values, diag + rstart, rend - rstart);
235: VecRestoreArray(v, &values);
236: PetscFree2(diag, work);
237: VecScale(v, Na->scale);
238: return 0;
239: }
241: PetscErrorCode MatGetDiagonalBlock_NormalHermitian(Mat N, Mat *D)
242: {
243: Mat_Normal *Na = (Mat_Normal *)N->data;
244: Mat M, A = Na->A;
246: MatGetDiagonalBlock(A, &M);
247: MatCreateNormalHermitian(M, &Na->D);
248: *D = Na->D;
249: return 0;
250: }
252: PetscErrorCode MatNormalGetMat_NormalHermitian(Mat A, Mat *M)
253: {
254: Mat_Normal *Aa = (Mat_Normal *)A->data;
256: *M = Aa->A;
257: return 0;
258: }
260: /*@
261: MatNormalHermitianGetMat - Gets the `Mat` object stored inside a `MATNORMALHERMITIAN`
263: Logically collective
265: Input Parameter:
266: . A - the `MATNORMALHERMITIAN` matrix
268: Output Parameter:
269: . M - the matrix object stored inside A
271: Level: intermediate
273: .seealso: `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
274: @*/
275: PetscErrorCode MatNormalHermitianGetMat(Mat A, Mat *M)
276: {
280: PetscUseMethod(A, "MatNormalGetMatHermitian_C", (Mat, Mat *), (A, M));
281: return 0;
282: }
284: PetscErrorCode MatConvert_NormalHermitian_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
285: {
286: Mat_Normal *Aa = (Mat_Normal *)A->data;
287: Mat B, conjugate;
288: PetscInt m, n, M, N;
290: MatGetSize(A, &M, &N);
291: MatGetLocalSize(A, &m, &n);
292: if (reuse == MAT_REUSE_MATRIX) {
293: B = *newmat;
294: MatProductReplaceMats(Aa->A, Aa->A, NULL, B);
295: } else {
296: MatProductCreate(Aa->A, Aa->A, NULL, &B);
297: MatProductSetType(B, MATPRODUCT_AtB);
298: MatProductSetFromOptions(B);
299: MatProductSymbolic(B);
300: MatSetOption(B, !PetscDefined(USE_COMPLEX) ? MAT_SYMMETRIC : MAT_HERMITIAN, PETSC_TRUE);
301: }
302: if (PetscDefined(USE_COMPLEX)) {
303: MatDuplicate(Aa->A, MAT_COPY_VALUES, &conjugate);
304: MatConjugate(conjugate);
305: MatProductReplaceMats(conjugate, Aa->A, NULL, B);
306: }
307: MatProductNumeric(B);
308: if (PetscDefined(USE_COMPLEX)) MatDestroy(&conjugate);
309: if (reuse == MAT_INPLACE_MATRIX) {
310: MatHeaderReplace(A, &B);
311: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
312: MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat);
313: return 0;
314: }
316: /*@
317: MatCreateNormalHermitian - Creates a new matrix object `MATNORMALHERMITIAN` that behaves like (A*)'*A.
319: Collective
321: Input Parameter:
322: . A - the (possibly rectangular complex) matrix
324: Output Parameter:
325: . N - the matrix that represents (A*)'*A
327: Level: intermediate
329: Note:
330: The product (A*)'*A is NOT actually formed! Rather the new matrix
331: object performs the matrix-vector product, `MatMult()`, by first multiplying by
332: A and then (A*)'
334: .seealso: `MATNORMAL`, `MATNORMALHERMITIAN`, `MatNormalHermitianGetMat()`
335: @*/
336: PetscErrorCode MatCreateNormalHermitian(Mat A, Mat *N)
337: {
338: PetscInt m, n;
339: Mat_Normal *Na;
340: VecType vtype;
342: MatGetLocalSize(A, &m, &n);
343: MatCreate(PetscObjectComm((PetscObject)A), N);
344: MatSetSizes(*N, n, n, PETSC_DECIDE, PETSC_DECIDE);
345: PetscObjectChangeTypeName((PetscObject)*N, MATNORMALHERMITIAN);
346: PetscLayoutReference(A->cmap, &(*N)->rmap);
347: PetscLayoutReference(A->cmap, &(*N)->cmap);
349: PetscNew(&Na);
350: (*N)->data = (void *)Na;
351: PetscObjectReference((PetscObject)A);
352: Na->A = A;
353: Na->scale = 1.0;
355: MatCreateVecs(A, NULL, &Na->w);
357: (*N)->ops->destroy = MatDestroy_NormalHermitian;
358: (*N)->ops->mult = MatMult_NormalHermitian;
359: (*N)->ops->multtranspose = MatMultHermitianTranspose_Normal;
360: (*N)->ops->multtransposeadd = MatMultHermitianTransposeAdd_Normal;
361: (*N)->ops->multadd = MatMultHermitianAdd_Normal;
362: (*N)->ops->getdiagonal = MatGetDiagonal_NormalHermitian;
363: (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_NormalHermitian;
364: (*N)->ops->scale = MatScale_NormalHermitian;
365: (*N)->ops->diagonalscale = MatDiagonalScale_NormalHermitian;
366: (*N)->ops->createsubmatrices = MatCreateSubMatrices_NormalHermitian;
367: (*N)->ops->permute = MatPermute_NormalHermitian;
368: (*N)->ops->duplicate = MatDuplicate_NormalHermitian;
369: (*N)->ops->copy = MatCopy_NormalHermitian;
370: (*N)->assembled = PETSC_TRUE;
371: (*N)->preallocated = PETSC_TRUE;
373: PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMatHermitian_C", MatNormalGetMat_NormalHermitian);
374: PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_seqaij_C", MatConvert_NormalHermitian_AIJ);
375: PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normalh_mpiaij_C", MatConvert_NormalHermitian_AIJ);
376: MatSetOption(*N, MAT_HERMITIAN, PETSC_TRUE);
377: MatGetVecType(A, &vtype);
378: MatSetVecType(*N, vtype);
379: #if defined(PETSC_HAVE_DEVICE)
380: MatBindToCPU(*N, A->boundtocpu);
381: #endif
382: return 0;
383: }