Actual source code: normm.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_Normal;
9: static PetscErrorCode MatIncreaseOverlap_Normal(Mat A, PetscInt is_max, IS is[], PetscInt ov)
10: {
11: Mat_Normal *a;
12: Mat pattern;
14: PetscFunctionBegin;
15: PetscCheck(ov >= 0, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap specified");
16: PetscCall(MatShellGetContext(A, &a));
17: PetscCall(MatProductCreate(a->A, a->A, NULL, &pattern));
18: PetscCall(MatProductSetType(pattern, MATPRODUCT_AtB));
19: PetscCall(MatProductSetFromOptions(pattern));
20: PetscCall(MatProductSymbolic(pattern));
21: PetscCall(MatIncreaseOverlap(pattern, is_max, is, ov));
22: PetscCall(MatDestroy(&pattern));
23: PetscFunctionReturn(PETSC_SUCCESS);
24: }
26: static PetscErrorCode MatCreateSubMatrices_Normal(Mat mat, PetscInt n, const IS irow[], const IS icol[], MatReuse scall, Mat *submat[])
27: {
28: Mat_Normal *a;
29: Mat B, *suba;
30: IS *row;
31: PetscScalar shift, scale;
32: PetscInt M;
34: PetscFunctionBegin;
35: PetscCheck(irow == icol, PetscObjectComm((PetscObject)mat), PETSC_ERR_SUP, "Not implemented");
36: PetscCall(MatShellGetScalingShifts(mat, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
37: PetscCall(MatShellGetContext(mat, &a));
38: B = a->A;
39: if (scall != MAT_REUSE_MATRIX) PetscCall(PetscCalloc1(n, submat));
40: PetscCall(MatGetSize(B, &M, NULL));
41: PetscCall(PetscMalloc1(n, &row));
42: PetscCall(ISCreateStride(PETSC_COMM_SELF, M, 0, 1, &row[0]));
43: PetscCall(ISSetIdentity(row[0]));
44: for (M = 1; M < n; ++M) row[M] = row[0];
45: PetscCall(MatCreateSubMatrices(B, n, row, icol, MAT_INITIAL_MATRIX, &suba));
46: for (M = 0; M < n; ++M) {
47: PetscCall(MatCreateNormal(suba[M], *submat + M));
48: PetscCall(MatShift((*submat)[M], shift));
49: PetscCall(MatScale((*submat)[M], scale));
50: }
51: PetscCall(ISDestroy(&row[0]));
52: PetscCall(PetscFree(row));
53: PetscCall(MatDestroySubMatrices(n, &suba));
54: PetscFunctionReturn(PETSC_SUCCESS);
55: }
57: static PetscErrorCode MatPermute_Normal(Mat A, IS rowp, IS colp, Mat *B)
58: {
59: Mat_Normal *a;
60: Mat C, Aa;
61: IS row;
62: PetscScalar shift, scale;
64: PetscFunctionBegin;
65: PetscCheck(rowp == colp, PetscObjectComm((PetscObject)A), PETSC_ERR_ARG_INCOMP, "Row permutation and column permutation must be the same");
66: PetscCall(MatShellGetScalingShifts(A, &shift, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
67: PetscCall(MatShellGetContext(A, &a));
68: Aa = a->A;
69: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)Aa), Aa->rmap->n, Aa->rmap->rstart, 1, &row));
70: PetscCall(ISSetIdentity(row));
71: PetscCall(MatPermute(Aa, row, colp, &C));
72: PetscCall(ISDestroy(&row));
73: PetscCall(MatCreateNormal(C, B));
74: PetscCall(MatDestroy(&C));
75: PetscCall(MatShift(*B, shift));
76: PetscCall(MatScale(*B, scale));
77: PetscFunctionReturn(PETSC_SUCCESS);
78: }
80: static PetscErrorCode MatDuplicate_Normal(Mat A, MatDuplicateOption op, Mat *B)
81: {
82: Mat_Normal *a;
83: Mat C;
85: PetscFunctionBegin;
86: PetscCall(MatShellGetContext(A, &a));
87: PetscCall(MatDuplicate(a->A, op, &C));
88: PetscCall(MatCreateNormal(C, B));
89: PetscCall(MatDestroy(&C));
90: if (op == MAT_COPY_VALUES) PetscCall(MatCopy(A, *B, SAME_NONZERO_PATTERN));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: static PetscErrorCode MatCopy_Normal(Mat A, Mat B, MatStructure str)
95: {
96: Mat_Normal *a, *b;
98: PetscFunctionBegin;
99: PetscCall(MatShellGetContext(A, &a));
100: PetscCall(MatShellGetContext(B, &b));
101: PetscCall(MatCopy(a->A, b->A, str));
102: PetscFunctionReturn(PETSC_SUCCESS);
103: }
105: static PetscErrorCode MatMult_Normal(Mat N, Vec x, Vec y)
106: {
107: Mat_Normal *Na;
109: PetscFunctionBegin;
110: PetscCall(MatShellGetContext(N, &Na));
111: PetscCall(MatMult(Na->A, x, Na->w));
112: PetscCall(MatMultTranspose(Na->A, Na->w, y));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode MatDestroy_Normal(Mat N)
117: {
118: Mat_Normal *Na;
120: PetscFunctionBegin;
121: PetscCall(MatShellGetContext(N, &Na));
122: PetscCall(MatDestroy(&Na->A));
123: PetscCall(MatDestroy(&Na->D));
124: PetscCall(VecDestroy(&Na->w));
125: PetscCall(PetscFree(Na));
126: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatNormalGetMat_C", NULL));
127: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_seqaij_C", NULL));
128: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_mpiaij_C", NULL));
129: #if defined(PETSC_HAVE_HYPRE)
130: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatConvert_normal_hypre_C", NULL));
131: #endif
132: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_seqdense_C", NULL));
133: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatProductSetFromOptions_normal_mpidense_C", NULL));
134: PetscCall(PetscObjectComposeFunction((PetscObject)N, "MatShellSetContext_C", NULL));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: /*
139: Slow, nonscalable version
140: */
141: static PetscErrorCode MatGetDiagonal_Normal(Mat N, Vec v)
142: {
143: Mat_Normal *Na;
144: Mat A;
145: PetscInt i, j, rstart, rend, nnz;
146: const PetscInt *cols;
147: PetscScalar *diag, *work, *values;
148: const PetscScalar *mvalues;
149: PetscMPIInt iN;
151: PetscFunctionBegin;
152: PetscCall(MatShellGetContext(N, &Na));
153: A = Na->A;
154: PetscCall(PetscMalloc2(A->cmap->N, &diag, A->cmap->N, &work));
155: PetscCall(PetscArrayzero(work, A->cmap->N));
156: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
157: for (i = rstart; i < rend; i++) {
158: PetscCall(MatGetRow(A, i, &nnz, &cols, &mvalues));
159: for (j = 0; j < nnz; j++) work[cols[j]] += mvalues[j] * mvalues[j];
160: PetscCall(MatRestoreRow(A, i, &nnz, &cols, &mvalues));
161: }
162: PetscCall(PetscMPIIntCast(A->cmap->N, &iN));
163: PetscCallMPI(MPIU_Allreduce(work, diag, iN, MPIU_SCALAR, MPIU_SUM, PetscObjectComm((PetscObject)N)));
164: rstart = N->cmap->rstart;
165: rend = N->cmap->rend;
166: PetscCall(VecGetArray(v, &values));
167: PetscCall(PetscArraycpy(values, diag + rstart, rend - rstart));
168: PetscCall(VecRestoreArray(v, &values));
169: PetscCall(PetscFree2(diag, work));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: static PetscErrorCode MatGetDiagonalBlock_Normal(Mat N, Mat *D)
174: {
175: Mat_Normal *Na;
176: Mat M, A;
178: PetscFunctionBegin;
179: PetscCall(MatShellGetContext(N, &Na));
180: A = Na->A;
181: PetscCall(MatGetDiagonalBlock(A, &M));
182: PetscCall(MatCreateNormal(M, &Na->D));
183: *D = Na->D;
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: static PetscErrorCode MatNormalGetMat_Normal(Mat A, Mat *M)
188: {
189: Mat_Normal *Aa;
191: PetscFunctionBegin;
192: PetscCall(MatShellGetContext(A, &Aa));
193: *M = Aa->A;
194: PetscFunctionReturn(PETSC_SUCCESS);
195: }
197: /*@
198: MatNormalGetMat - Gets the `Mat` object stored inside a `MATNORMAL`
200: Logically Collective
202: Input Parameter:
203: . A - the `MATNORMAL` matrix
205: Output Parameter:
206: . M - the matrix object stored inside `A`
208: Level: intermediate
210: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MATNORMALHERMITIAN`, `MatCreateNormal()`
211: @*/
212: PetscErrorCode MatNormalGetMat(Mat A, Mat *M)
213: {
214: PetscFunctionBegin;
217: PetscAssertPointer(M, 2);
218: PetscUseMethod(A, "MatNormalGetMat_C", (Mat, Mat *), (A, M));
219: PetscFunctionReturn(PETSC_SUCCESS);
220: }
222: static PetscErrorCode MatConvert_Normal_AIJ(Mat A, MatType newtype, MatReuse reuse, Mat *newmat)
223: {
224: Mat_Normal *Aa;
225: Mat B;
226: PetscInt m, n, M, N;
228: PetscFunctionBegin;
229: PetscCall(MatShellGetContext(A, &Aa));
230: PetscCall(MatGetSize(A, &M, &N));
231: PetscCall(MatGetLocalSize(A, &m, &n));
232: if (reuse == MAT_REUSE_MATRIX) {
233: B = *newmat;
234: PetscCall(MatProductReplaceMats(Aa->A, Aa->A, NULL, B));
235: } else {
236: PetscCall(MatProductCreate(Aa->A, Aa->A, NULL, &B));
237: PetscCall(MatProductSetType(B, MATPRODUCT_AtB));
238: PetscCall(MatProductSetFromOptions(B));
239: PetscCall(MatProductSymbolic(B));
240: PetscCall(MatSetOption(B, MAT_SYMMETRIC, PETSC_TRUE));
241: }
242: PetscCall(MatProductNumeric(B));
243: if (reuse == MAT_INPLACE_MATRIX) {
244: PetscCall(MatHeaderReplace(A, &B));
245: } else if (reuse == MAT_INITIAL_MATRIX) *newmat = B;
246: PetscCall(MatConvert(*newmat, MATAIJ, MAT_INPLACE_MATRIX, newmat));
247: PetscFunctionReturn(PETSC_SUCCESS);
248: }
250: #if defined(PETSC_HAVE_HYPRE)
251: static PetscErrorCode MatConvert_Normal_HYPRE(Mat A, MatType type, MatReuse reuse, Mat *B)
252: {
253: PetscFunctionBegin;
254: if (reuse == MAT_INITIAL_MATRIX) {
255: PetscCall(MatConvert(A, MATAIJ, reuse, B));
256: PetscCall(MatConvert(*B, type, MAT_INPLACE_MATRIX, B));
257: } else PetscCall(MatConvert_Basic(A, type, reuse, B)); /* fall back to basic convert */
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
260: #endif
262: typedef struct {
263: Mat work[2];
264: } Normal_Dense;
266: static PetscErrorCode MatProductNumeric_Normal_Dense(Mat C)
267: {
268: Mat A, B;
269: Normal_Dense *contents;
270: Mat_Normal *a;
271: Vec right;
272: PetscScalar *array, scale;
274: PetscFunctionBegin;
275: MatCheckProduct(C, 1);
276: A = C->product->A;
277: B = C->product->B;
278: PetscCall(MatShellGetContext(A, &a));
279: contents = (Normal_Dense *)C->product->data;
280: PetscCheck(contents, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data empty");
281: PetscCall(MatShellGetScalingShifts(A, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, &right, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
282: if (right) {
283: PetscCall(MatCopy(B, C, SAME_NONZERO_PATTERN));
284: PetscCall(MatDiagonalScale(C, right, NULL));
285: }
286: PetscCall(MatProductNumeric(contents->work[0]));
287: PetscCall(MatDenseGetArrayWrite(C, &array));
288: PetscCall(MatDensePlaceArray(contents->work[1], array));
289: PetscCall(MatProductNumeric(contents->work[1]));
290: PetscCall(MatDenseRestoreArrayWrite(C, &array));
291: PetscCall(MatDenseResetArray(contents->work[1]));
292: PetscCall(MatSetOption(C, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE));
293: PetscCall(MatAssemblyBegin(C, MAT_FINAL_ASSEMBLY));
294: PetscCall(MatAssemblyEnd(C, MAT_FINAL_ASSEMBLY));
295: PetscCall(MatScale(C, scale));
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode MatNormal_DenseDestroy(void *ctx)
300: {
301: Normal_Dense *contents = (Normal_Dense *)ctx;
303: PetscFunctionBegin;
304: PetscCall(MatDestroy(contents->work));
305: PetscCall(MatDestroy(contents->work + 1));
306: PetscCall(PetscFree(contents));
307: PetscFunctionReturn(PETSC_SUCCESS);
308: }
310: static PetscErrorCode MatProductSymbolic_Normal_Dense(Mat C)
311: {
312: Mat A, B;
313: Normal_Dense *contents = NULL;
314: Mat_Normal *a;
315: Vec right;
316: PetscScalar *array, scale;
317: PetscInt n, N, m, M;
319: PetscFunctionBegin;
320: MatCheckProduct(C, 1);
321: PetscCheck(!C->product->data, PetscObjectComm((PetscObject)C), PETSC_ERR_PLIB, "Product data not empty");
322: A = C->product->A;
323: B = C->product->B;
324: PetscCall(MatShellGetScalingShifts(A, (PetscScalar *)MAT_SHELL_NOT_ALLOWED, &scale, (Vec *)MAT_SHELL_NOT_ALLOWED, (Vec *)MAT_SHELL_NOT_ALLOWED, &right, (Mat *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED, (IS *)MAT_SHELL_NOT_ALLOWED));
325: PetscCall(MatShellGetContext(A, &a));
326: PetscCall(MatGetLocalSize(C, &m, &n));
327: PetscCall(MatGetSize(C, &M, &N));
328: if (m == PETSC_DECIDE || n == PETSC_DECIDE || M == PETSC_DECIDE || N == PETSC_DECIDE) {
329: PetscCall(MatGetLocalSize(B, NULL, &n));
330: PetscCall(MatGetSize(B, NULL, &N));
331: PetscCall(MatGetLocalSize(A, &m, NULL));
332: PetscCall(MatGetSize(A, &M, NULL));
333: PetscCall(MatSetSizes(C, m, n, M, N));
334: }
335: PetscCall(MatSetType(C, ((PetscObject)B)->type_name));
336: PetscCall(MatSetUp(C));
337: PetscCall(PetscNew(&contents));
338: C->product->data = contents;
339: C->product->destroy = MatNormal_DenseDestroy;
340: if (right) PetscCall(MatProductCreate(a->A, C, NULL, contents->work));
341: else PetscCall(MatProductCreate(a->A, B, NULL, contents->work));
342: PetscCall(MatProductSetType(contents->work[0], MATPRODUCT_AB));
343: PetscCall(MatProductSetFromOptions(contents->work[0]));
344: PetscCall(MatProductSymbolic(contents->work[0]));
345: PetscCall(MatProductCreate(a->A, contents->work[0], NULL, contents->work + 1));
346: PetscCall(MatProductSetType(contents->work[1], MATPRODUCT_AtB));
347: PetscCall(MatProductSetFromOptions(contents->work[1]));
348: PetscCall(MatProductSymbolic(contents->work[1]));
349: PetscCall(MatDenseGetArrayWrite(C, &array));
350: PetscCall(MatSeqDenseSetPreallocation(contents->work[1], array));
351: PetscCall(MatMPIDenseSetPreallocation(contents->work[1], array));
352: PetscCall(MatDenseRestoreArrayWrite(C, &array));
353: C->ops->productnumeric = MatProductNumeric_Normal_Dense;
354: PetscFunctionReturn(PETSC_SUCCESS);
355: }
357: static PetscErrorCode MatProductSetFromOptions_Normal_Dense(Mat C)
358: {
359: Mat_Product *product = C->product;
361: PetscFunctionBegin;
362: if (product->type == MATPRODUCT_AB) C->ops->productsymbolic = MatProductSymbolic_Normal_Dense;
363: PetscFunctionReturn(PETSC_SUCCESS);
364: }
366: /*MC
367: MATNORMAL - a matrix that behaves like A'*A for `MatMult()` while only containing A
369: Level: intermediate
371: Developer Notes:
372: This is implemented on top of `MATSHELL` to get support for scaling and shifting without requiring duplicate code
374: Users can not call `MatShellSetOperation()` operations on this class, there is some error checking for that incorrect usage
376: .seealso: [](ch_matrices), `Mat`, `MatCreateNormal()`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
377: M*/
379: /*@
380: MatCreateNormal - Creates a new `MATNORMAL` matrix object that behaves like A'*A.
382: Collective
384: Input Parameter:
385: . A - the (possibly rectangular) matrix
387: Output Parameter:
388: . N - the matrix that represents A'*A
390: Level: intermediate
392: Notes:
393: The product A'*A is NOT actually formed! Rather the new matrix
394: object performs the matrix-vector product, `MatMult()`, by first multiplying by
395: A and then A'
397: .seealso: [](ch_matrices), `Mat`, `MATNORMAL`, `MatMult()`, `MatNormalGetMat()`, `MATNORMALHERMITIAN`, `MatCreateNormalHermitian()`
398: @*/
399: PetscErrorCode MatCreateNormal(Mat A, Mat *N)
400: {
401: Mat_Normal *Na;
402: VecType vtype;
404: PetscFunctionBegin;
405: PetscCall(MatCreate(PetscObjectComm((PetscObject)A), N));
406: PetscCall(PetscLayoutReference(A->cmap, &(*N)->rmap));
407: PetscCall(PetscLayoutReference(A->cmap, &(*N)->cmap));
408: PetscCall(MatSetType(*N, MATSHELL));
409: PetscCall(PetscNew(&Na));
410: PetscCall(MatShellSetContext(*N, Na));
411: PetscCall(PetscObjectReference((PetscObject)A));
412: Na->A = A;
413: PetscCall(MatCreateVecs(A, NULL, &Na->w));
415: PetscCall(MatSetBlockSizes(*N, PetscAbs(A->cmap->bs), PetscAbs(A->rmap->bs)));
416: PetscCall(MatShellSetOperation(*N, MATOP_DESTROY, (void (*)(void))MatDestroy_Normal));
417: PetscCall(MatShellSetOperation(*N, MATOP_MULT, (void (*)(void))MatMult_Normal));
418: PetscCall(MatShellSetOperation(*N, MATOP_MULT_TRANSPOSE, (void (*)(void))MatMult_Normal));
419: PetscCall(MatShellSetOperation(*N, MATOP_DUPLICATE, (void (*)(void))MatDuplicate_Normal));
420: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL, (void (*)(void))MatGetDiagonal_Normal));
421: PetscCall(MatShellSetOperation(*N, MATOP_GET_DIAGONAL_BLOCK, (void (*)(void))MatGetDiagonalBlock_Normal));
422: PetscCall(MatShellSetOperation(*N, MATOP_COPY, (void (*)(void))MatCopy_Normal));
423: (*N)->ops->increaseoverlap = MatIncreaseOverlap_Normal;
424: (*N)->ops->createsubmatrices = MatCreateSubMatrices_Normal;
425: (*N)->ops->permute = MatPermute_Normal;
427: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatNormalGetMat_C", MatNormalGetMat_Normal));
428: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_seqaij_C", MatConvert_Normal_AIJ));
429: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_mpiaij_C", MatConvert_Normal_AIJ));
430: #if defined(PETSC_HAVE_HYPRE)
431: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatConvert_normal_hypre_C", MatConvert_Normal_HYPRE));
432: #endif
433: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_normal_seqdense_C", MatProductSetFromOptions_Normal_Dense));
434: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatProductSetFromOptions_normal_mpidense_C", MatProductSetFromOptions_Normal_Dense));
435: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContext_C", MatShellSetContext_Immutable));
436: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetContextDestroy_C", MatShellSetContextDestroy_Immutable));
437: PetscCall(PetscObjectComposeFunction((PetscObject)*N, "MatShellSetManageScalingShifts_C", MatShellSetManageScalingShifts_Immutable));
438: PetscCall(MatSetOption(*N, MAT_SYMMETRIC, PETSC_TRUE));
439: PetscCall(MatGetVecType(A, &vtype));
440: PetscCall(MatSetVecType(*N, vtype));
441: #if defined(PETSC_HAVE_DEVICE)
442: PetscCall(MatBindToCPU(*N, A->boundtocpu));
443: #endif
444: PetscCall(MatSetUp(*N));
445: PetscCall(PetscObjectChangeTypeName((PetscObject)*N, MATNORMAL));
446: PetscFunctionReturn(PETSC_SUCCESS);
447: }