Actual source code: kaczmarz.c
1: #include <petsc/private/pcimpl.h>
3: typedef struct {
4: PetscReal lambda; /* damping parameter */
5: PetscBool symmetric; /* apply the projections symmetrically */
6: } PC_Kaczmarz;
8: static PetscErrorCode PCDestroy_Kaczmarz(PC pc)
9: {
10: PetscFunctionBegin;
11: PetscCall(PetscFree(pc->data));
12: PetscFunctionReturn(PETSC_SUCCESS);
13: }
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;
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: }
70: static PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc, PetscOptionItems *PetscOptionsObject)
71: {
72: PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
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: }
82: static PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer)
83: {
84: PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
85: PetscBool iascii;
87: PetscFunctionBegin;
88: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
89: if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, " lambda = %g\n", (double)jac->lambda));
90: PetscFunctionReturn(PETSC_SUCCESS);
91: }
93: /*MC
94: PCKACZMARZ - Kaczmarz iteration {cite}`kaczmarz1937angenaherte`
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
100: Level: beginner
102: Note:
103: In parallel this is block-Jacobi with Kaczmarz inner solve on each processor.
105: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCBJACOBI`
106: M*/
108: PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
109: {
110: PC_Kaczmarz *jac;
112: PetscFunctionBegin;
113: PetscCall(PetscNew(&jac));
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: }