1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
2
3 typedef struct {
4 PetscReal lambda; /* damping parameter */
5 PetscBool symmetric; /* apply the projections symmetrically */
6 } PC_Kaczmarz;
7
PCDestroy_Kaczmarz(PC pc)8 static PetscErrorCode PCDestroy_Kaczmarz(PC pc)
9 {
10 PetscFunctionBegin;
11 PetscCall(PetscFree(pc->data));
12 PetscFunctionReturn(PETSC_SUCCESS);
13 }
14
PCApply_Kaczmarz(PC pc,Vec x,Vec y)15 static PetscErrorCode PCApply_Kaczmarz(PC pc, Vec x, Vec y)
16 {
17 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
18 PetscInt xs, xe, ys, ye, ncols, i, j;
19 const PetscInt *cols;
20 const PetscScalar *vals, *xarray;
21 PetscScalar r;
22 PetscReal anrm;
23 PetscScalar *yarray;
24 PetscReal lambda = jac->lambda;
25
26 PetscFunctionBegin;
27 PetscCall(MatGetOwnershipRange(pc->pmat, &xs, &xe));
28 PetscCall(MatGetOwnershipRangeColumn(pc->pmat, &ys, &ye));
29 PetscCall(VecSet(y, 0.));
30 PetscCall(VecGetArrayRead(x, &xarray));
31 PetscCall(VecGetArray(y, &yarray));
32 for (i = xs; i < xe; i++) {
33 /* get the maximum row width and row norms */
34 PetscCall(MatGetRow(pc->pmat, i, &ncols, &cols, &vals));
35 r = xarray[i - xs];
36 anrm = 0.;
37 for (j = 0; j < ncols; j++) {
38 if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
39 anrm += PetscRealPart(PetscSqr(vals[j]));
40 }
41 if (anrm > 0.) {
42 for (j = 0; j < ncols; j++) {
43 if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
44 }
45 }
46 PetscCall(MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals));
47 }
48 if (jac->symmetric) {
49 for (i = xe - 1; i >= xs; i--) {
50 PetscCall(MatGetRow(pc->pmat, i, &ncols, &cols, &vals));
51 r = xarray[i - xs];
52 anrm = 0.;
53 for (j = 0; j < ncols; j++) {
54 if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
55 anrm += PetscRealPart(PetscSqr(vals[j]));
56 }
57 if (anrm > 0.) {
58 for (j = 0; j < ncols; j++) {
59 if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
60 }
61 }
62 PetscCall(MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals));
63 }
64 }
65 PetscCall(VecRestoreArray(y, &yarray));
66 PetscCall(VecRestoreArrayRead(x, &xarray));
67 PetscFunctionReturn(PETSC_SUCCESS);
68 }
69
PCSetFromOptions_Kaczmarz(PC pc,PetscOptionItems PetscOptionsObject)70 static PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc, PetscOptionItems PetscOptionsObject)
71 {
72 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
73
74 PetscFunctionBegin;
75 PetscOptionsHeadBegin(PetscOptionsObject, "Kaczmarz options");
76 PetscCall(PetscOptionsReal("-pc_kaczmarz_lambda", "relaxation factor (0 < lambda)", "", jac->lambda, &jac->lambda, NULL));
77 PetscCall(PetscOptionsBool("-pc_kaczmarz_symmetric", "apply row projections symmetrically", "", jac->symmetric, &jac->symmetric, NULL));
78 PetscOptionsHeadEnd();
79 PetscFunctionReturn(PETSC_SUCCESS);
80 }
81
PCView_Kaczmarz(PC pc,PetscViewer viewer)82 static PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer)
83 {
84 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
85 PetscBool isascii;
86
87 PetscFunctionBegin;
88 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
89 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, " lambda = %g\n", (double)jac->lambda));
90 PetscFunctionReturn(PETSC_SUCCESS);
91 }
92
93 /*MC
94 PCKACZMARZ - Kaczmarz iteration {cite}`kaczmarz1937angenaherte`
95
96 Options Database Keys:
97 + -pc_kaczmarz_lambda <lambda> - Sets damping parameter defaults to 1.0
98 - -pc_kaczmarz_symmetric - Apply the row projections symmetrically
99
100 Level: beginner
101
102 Note:
103 In parallel this is block-Jacobi with Kaczmarz inner solve on each processor.
104
105 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCBJACOBI`
106 M*/
107
PCCreate_Kaczmarz(PC pc)108 PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
109 {
110 PC_Kaczmarz *jac;
111
112 PetscFunctionBegin;
113 PetscCall(PetscNew(&jac));
114
115 pc->ops->apply = PCApply_Kaczmarz;
116 pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz;
117 pc->ops->setup = NULL;
118 pc->ops->view = PCView_Kaczmarz;
119 pc->ops->destroy = PCDestroy_Kaczmarz;
120 pc->data = (void *)jac;
121 jac->lambda = 1.0;
122 jac->symmetric = PETSC_FALSE;
123 PetscFunctionReturn(PETSC_SUCCESS);
124 }
125