Actual source code: normm.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_Normal(Mat inA, PetscScalar scale)
12: {
13: Mat_Normal *a = (Mat_Normal *)inA->data;
15: a->scale *= scale;
16: return 0;
17: }
19: PetscErrorCode MatDiagonalScale_Normal(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 MatIncreaseOverlap_Normal(Mat A, PetscInt is_max, IS is[], PetscInt ov)
43: {
44: Mat_Normal *a = (Mat_Normal *)A->data;
45: Mat pattern;
48: MatProductCreate(a->A, a->A, NULL, &pattern);
49: MatProductSetType(pattern, MATPRODUCT_AtB);
50: MatProductSetFromOptions(pattern);
51: MatProductSymbolic(pattern);
52: MatIncreaseOverlap(pattern, is_max, is, ov);
53: MatDestroy(&pattern);
54: return 0;
55: }
57: PetscErrorCode MatCreateSubMatrices_Normal(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
58: {
59: Mat_Normal *a = (Mat_Normal *)mat->data;
60: Mat B = a->A, *suba;
61: IS *row;
62: PetscInt M;
65: if (scall != MAT_REUSE_MATRIX) PetscCalloc1(n, submat);
66: MatGetSize(B, &M, NULL);
67: PetscMalloc1(n, &row);
68: ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]);
69: ISSetIdentity(row[0]);
70: for (M = 1; M < n; ++M) row[M] = row[0];
71: MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba);
72: for (M = 0; M < n; ++M) {
73: MatCreateNormal(suba[M], *submat + M);
74: ((Mat_Normal *)(*submat)[M]->data)->scale = a->scale;
75: }
76: ISDestroy(&row[0]);
77: PetscFree(row);
78: MatDestroySubMatrices(n, &suba);
79: return 0;
80: }
82: PetscErrorCode MatPermute_Normal(Mat A, IS rowp, IS colp, Mat *B)
83: {
84: Mat_Normal *a = (Mat_Normal *)A->data;
85: Mat C, Aa = a->A;
86: IS row;
89: ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row);
90: ISSetIdentity(row);
91: MatPermute(Aa, row, colp, &C);
92: ISDestroy(&row);
93: MatCreateNormal(C, B);
94: MatDestroy(&C);
95: return 0;
96: }
98: PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
99: {
100: Mat_Normal *a = (Mat_Normal *)A->data;
101: Mat C;
104: MatDuplicate(a->A, op, &C);
105: MatCreateNormal(C, B);
106: MatDestroy(&C);
107: if (op == MAT_COPY_VALUES) ((Mat_Normal *)(*B)->data)->scale = a->scale;
108: return 0;
109: }
111: PetscErrorCode MatCopy_Normal(Mat A, Mat B, MatStructure str)
112: {
113: Mat_Normal *a = (Mat_Normal *)A->data, *b = (Mat_Normal *)B->data;
116: MatCopy(a->A, b->A, str);
117: b->scale = a->scale;
118: VecDestroy(&b->left);
119: VecDestroy(&b->right);
120: VecDestroy(&b->leftwork);
121: VecDestroy(&b->rightwork);
122: return 0;
123: }
125: PetscErrorCode MatMult_Normal(Mat N, Vec x, Vec y)
126: {
127: Mat_Normal *Na = (Mat_Normal *)N->data;
128: Vec in;
130: in = x;
131: if (Na->right) {
132: if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
133: VecPointwiseMult(Na->rightwork, Na->right, in);
134: in = Na->rightwork;
135: }
136: MatMult(Na->A, in, Na->w);
137: MatMultTranspose(Na->A, Na->w, y);
138: if (Na->left) VecPointwiseMult(y, Na->left, y);
139: VecScale(y, Na->scale);
140: return 0;
141: }
143: PetscErrorCode MatMultAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
144: {
145: Mat_Normal *Na = (Mat_Normal *)N->data;
146: Vec in;
148: in = v1;
149: if (Na->right) {
150: if (!Na->rightwork) VecDuplicate(Na->right, &Na->rightwork);
151: VecPointwiseMult(Na->rightwork, Na->right, in);
152: in = Na->rightwork;
153: }
154: MatMult(Na->A, in, Na->w);
155: VecScale(Na->w, Na->scale);
156: if (Na->left) {
157: MatMultTranspose(Na->A, Na->w, v3);
158: VecPointwiseMult(v3, Na->left, v3);
159: VecAXPY(v3, 1.0, v2);
160: } else {
161: MatMultTransposeAdd(Na->A, Na->w, v2, v3);
162: }
163: return 0;
164: }
166: PetscErrorCode MatMultTranspose_Normal(Mat N, Vec x, Vec y)
167: {
168: Mat_Normal *Na = (Mat_Normal *)N->data;
169: Vec in;
171: in = x;
172: if (Na->left) {
173: if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
174: VecPointwiseMult(Na->leftwork, Na->left, in);
175: in = Na->leftwork;
176: }
177: MatMult(Na->A, in, Na->w);
178: MatMultTranspose(Na->A, Na->w, y);
179: if (Na->right) VecPointwiseMult(y, Na->right, y);
180: VecScale(y, Na->scale);
181: return 0;
182: }
184: PetscErrorCode MatMultTransposeAdd_Normal(Mat N, Vec v1, Vec v2, Vec v3)
185: {
186: Mat_Normal *Na = (Mat_Normal *)N->data;
187: Vec in;
189: in = v1;
190: if (Na->left) {
191: if (!Na->leftwork) VecDuplicate(Na->left, &Na->leftwork);
192: VecPointwiseMult(Na->leftwork, Na->left, in);
193: in = Na->leftwork;
194: }
195: MatMult(Na->A, in, Na->w);
196: VecScale(Na->w, Na->scale);
197: if (Na->right) {
198: MatMultTranspose(Na->A, Na->w, v3);
199: VecPointwiseMult(v3, Na->right, v3);
200: VecAXPY(v3, 1.0, v2);
201: } else {
202: MatMultTransposeAdd(Na->A, Na->w, v2, v3);
203: }
204: return 0;
205: }
207: PetscErrorCode MatDestroy_Normal(Mat N)
208: {
209: Mat_Normal *Na = (Mat_Normal *)N->data;
211: MatDestroy(&Na->A);
212: MatDestroy(&Na->D);
213: VecDestroy(&Na->w);
214: VecDestroy(&Na->left);
215: VecDestroy(&Na->right);
216: VecDestroy(&Na->leftwork);
217: VecDestroy(&Na->rightwork);
218: PetscFree(N->data);
219: PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL);
220: PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_seqaij_C", NULL);
221: PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_mpiaij_C", NULL);
222: PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_seqdense_C", NULL);
223: PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_mpidense_C", NULL);
224: PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_dense_C", NULL);
225: return 0;
226: }
228: /*
229: Slow, nonscalable version
230: */
231: PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v)
232: {
233: Mat_Normal *Na = (Mat_Normal *)N->data;
234: Mat A = Na->A;
235: PetscInt i, j, rstart, rend, nnz;
236: const PetscInt *cols;
237: PetscScalar *diag, *work, *values;
238: const PetscScalar *mvalues;
240: PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work);
241: PetscArrayzero(work, A->cmap->N);
242: MatGetOwnershipRange(A, &rstart, &rend);
243: for (i = rstart; i < rend; i++) {
244: MatGetRow(A, i, &nnz, &cols, &mvalues);
245: for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j];
246: MatRestoreRow(A, i, &nnz, &cols, &mvalues);
247: }
248: MPIU_Allreduce(work, diag, A->cmap->N, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N));
249: rstart = N->cmap->rstart;
250: rend = N->cmap->rend;
251: VecGetArray(v, &values);
252: PetscArraycpy(values, diag + rstart, rend - rstart);
253: VecRestoreArray(v, &values);
254: PetscFree2(diag, work);
255: VecScale(v, Na->scale);
256: return 0;
257: }
259: PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
260: {
261: Mat_Normal *Na = (Mat_Normal *)N->data;
262: Mat M, A = Na->A;
264: MatGetDiagonalBlock(A, &M);
265: MatCreateNormal(M, &Na->D);
266: *D = Na->D;
267: return 0;
268: }
270: PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M)
271: {
272: Mat_Normal *Aa = (Mat_Normal *)A->data;
274: *M = Aa->A;
275: return 0;
276: }
278: /*@
279: MatNormalGetMat - Gets the `Mat` object stored inside a `MATNORMAL`
281: Logically collective on A
283: Input Parameter:
284: . A - the `MATNORMAL` matrix
286: Output Parameter:
287: . M - the matrix object stored inside A
289: Level: intermediate
291: .seealso: `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()`
292: @*/
293: PetscErrorCode MatNormalGetMat(Mat A, Mat *M)
294: {
298: PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M));
299: return 0;
300: }
302: PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
303: {
304: Mat_Normal *Aa = (Mat_Normal *)A->data;
305: Mat B;
306: PetscInt m, n, M, N;
308: MatGetSize(A, &M, &N);
309: MatGetLocalSize(A, &m, &n);
310: if (reuse == MAT_REUSE_MATRIX) {
311: B = *newmat;
312: MatProductReplaceMats(Aa->A, Aa->A, NULL, B);
313: } else {
314: MatProductCreate(Aa->A, Aa->A, NULL, &B);
315: MatProductSetType(B, MATPRODUCT_AtB);
316: MatProductSetFromOptions(B);
317: MatProductSymbolic(B);
318: MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE);
319: }
320: MatProductNumeric(B);
321: if (reuse == MAT_INPLACE_MATRIX) {
322: MatHeaderReplace(A, &B);
323: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
324: MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat);
325: return 0;
326: }
328: typedef struct {
329: Mat work[2];
330: } Normal_Dense;
332: PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
333: {
334: Mat A, B;
335: Normal_Dense *contents;
336: Mat_Normal *a;
337: PetscScalar *array;
339: MatCheckProduct(C, 1);
340: A = C->product->A;
341: a = (Mat_Normal *)A->data;
342: B = C->product->B;
343: contents = (Normal_Dense *)C->product->data;
345: if (a->right) {
346: MatCopy(B, C, SAME_NONZERO_PATTERN);
347: MatDiagonalScale(C, a->right, NULL);
348: }
349: MatProductNumeric(contents->work[0]);
350: MatDenseGetArrayWrite(C, &array);
351: MatDensePlaceArray(contents->work[1], array);
352: MatProductNumeric(contents->work[1]);
353: MatDenseRestoreArrayWrite(C, &array);
354: MatDenseResetArray(contents->work[1]);
355: MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
356: MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY);
357: MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY);
358: MatScale(C, a->scale);
359: return 0;
360: }
362: PetscErrorCode MatNormal_DenseDestroy(void *ctx)
363: {
364: Normal_Dense *contents = (Normal_Dense *)ctx;
366: MatDestroy(contents->work);
367: MatDestroy(contents->work + 1);
368: PetscFree(contents);
369: return 0;
370: }
372: PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
373: {
374: Mat A, B;
375: Normal_Dense *contents = NULL;
376: Mat_Normal *a;
377: PetscScalar *array;
378: PetscInt n, N, m, M;
380: MatCheckProduct(C, 1);
382: A = C->product->A;
383: a = (Mat_Normal *)A->data;
385: B = C->product->B;
386: MatGetLocalSize(C, &m, &n);
387: MatGetSize(C, &M, &N);
388: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
389: MatGetLocalSize(B, NULL, &n);
390: MatGetSize(B, NULL, &N);
391: MatGetLocalSize(A, &m, NULL);
392: MatGetSize(A, &M, NULL);
393: MatSetSizes(C, m, n, M, N);
394: }
395: MatSetType(C, ((PetscObject)B)->type_name);
396: MatSetUp(C);
397: PetscNew(&contents);
398: C->product->data = contents;
399: C->product->destroy = MatNormal_DenseDestroy;
400: if (a->right) {
401: MatProductCreate(a->A, C, NULL, contents->work);
402: } else {
403: MatProductCreate(a->A, B, NULL, contents->work);
404: }
405: MatProductSetType(contents->work[0], MATPRODUCT_AB);
406: MatProductSetFromOptions(contents->work[0]);
407: MatProductSymbolic(contents->work[0]);
408: MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1);
409: MatProductSetType(contents->work[1], MATPRODUCT_AtB);
410: MatProductSetFromOptions(contents->work[1]);
411: MatProductSymbolic(contents->work[1]);
412: MatDenseGetArrayWrite(C, &array);
413: MatSeqDenseSetPreallocation(contents->work[1], array);
414: MatMPIDenseSetPreallocation(contents->work[1], array);
415: MatDenseRestoreArrayWrite(C, &array);
416: C->ops->productnumeric = MatProductNumeric_Normal_Dense;
417: return 0;
418: }
420: PetscErrorCode MatProductSetFromOptions_Normal_Dense_AB(Mat C)
421: {
422: C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
423: return 0;
424: }
426: PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
427: {
428: Mat_Product *product = C->product;
430: if (product->type == MATPRODUCT_AB) MatProductSetFromOptions_Normal_Dense_AB(C);
431: return 0;
432: }
434: /*@
435: MatCreateNormal - Creates a new `MATNORMAL` matrix object that behaves like A'*A.
437: Collective
439: Input Parameter:
440: . A - the (possibly rectangular) matrix
442: Output Parameter:
443: . N - the matrix that represents A'*A
445: Level: intermediate
447: Notes:
448: The product A'*A is NOT actually formed! Rather the new matrix
449: object performs the matrix-vector product, `MatMult()`, by first multiplying by
450: A and then A'
452: .seealso: `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`,
453: @*/
454: PetscErrorCode MatCreateNormal(Mat A, Mat *N)
455: {
456: PetscInt n, nn;
457: Mat_Normal *Na;
458: VecType vtype;
460: MatGetSize(A, NULL, &nn);
461: MatGetLocalSize(A, NULL, &n);
462: MatCreate(PetscObjectComm((PetscObject)A), N);
463: MatSetSizes(*N, n, n, nn, nn);
464: PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL);
465: PetscLayoutReference(A->cmap, &(*N)->rmap);
466: PetscLayoutReference(A->cmap, &(*N)->cmap);
468: PetscNew(&Na);
469: (*N)->data = (void *)Na;
470: PetscObjectReference((PetscObject)A);
471: Na->A = A;
472: Na->scale = 1.0;
474: MatCreateVecs(A, NULL, &Na->w);
476: (*N)->ops->destroy = MatDestroy_Normal;
477: (*N)->ops->mult = MatMult_Normal;
478: (*N)->ops->multtranspose = MatMultTranspose_Normal;
479: (*N)->ops->multtransposeadd = MatMultTransposeAdd_Normal;
480: (*N)->ops->multadd = MatMultAdd_Normal;
481: (*N)->ops->getdiagonal = MatGetDiagonal_Normal;
482: (*N)->ops->getdiagonalblock = MatGetDiagonalBlock_Normal;
483: (*N)->ops->scale = MatScale_Normal;
484: (*N)->ops->diagonalscale = MatDiagonalScale_Normal;
485: (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal;
486: (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
487: (*N)->ops->permute = MatPermute_Normal;
488: (*N)->ops->duplicate = MatDuplicate_Normal;
489: (*N)->ops->copy = MatCopy_Normal;
490: (*N)->assembled = PETSC_TRUE;
491: (*N)->preallocated = PETSC_TRUE;
493: PetscObjectComposeFunction((PetscObject)(*N), "MatNormalGetMat_C", MatNormalGetMat_Normal);
494: PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ);
495: PetscObjectComposeFunction((PetscObject)(*N), "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ);
496: PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense);
497: PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense);
498: PetscObjectComposeFunction((PetscObject)(*N), "MatProductSetFromOptions_normal_dense_C", MatProductSetFromOptions_Normal_Dense);
499: MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE);
500: MatGetVecType(A, &vtype);
501: MatSetVecType(*N, vtype);
502: #if defined(PETSC_HAVE_DEVICE)
503: MatBindToCPU(*N, A->boundtocpu);
504: #endif
505: return 0;
506: }