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: PetscFree(pc->data);
11: return 0;
12: }
14: static PetscErrorCode PCApply_Kaczmarz(PC pc, Vec x, Vec y)
15: {
16: PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
17: PetscInt xs, xe, ys, ye, ncols, i, j;
18: const PetscInt *cols;
19: const PetscScalar *vals, *xarray;
20: PetscScalar r;
21: PetscReal anrm;
22: PetscScalar *yarray;
23: PetscReal lambda = jac->lambda;
25: MatGetOwnershipRange(pc->pmat, &xs, &xe);
26: MatGetOwnershipRangeColumn(pc->pmat, &ys, &ye);
27: VecSet(y, 0.);
28: VecGetArrayRead(x, &xarray);
29: VecGetArray(y, &yarray);
30: for (i = xs; i < xe; i++) {
31: /* get the maximum row width and row norms */
32: MatGetRow(pc->pmat, i, &ncols, &cols, &vals);
33: r = xarray[i - xs];
34: anrm = 0.;
35: for (j = 0; j < ncols; j++) {
36: if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
37: anrm += PetscRealPart(PetscSqr(vals[j]));
38: }
39: if (anrm > 0.) {
40: for (j = 0; j < ncols; j++) {
41: if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
42: }
43: }
44: MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals);
45: }
46: if (jac->symmetric) {
47: for (i = xe - 1; i >= xs; i--) {
48: MatGetRow(pc->pmat, i, &ncols, &cols, &vals);
49: r = xarray[i - xs];
50: anrm = 0.;
51: for (j = 0; j < ncols; j++) {
52: if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
53: anrm += PetscRealPart(PetscSqr(vals[j]));
54: }
55: if (anrm > 0.) {
56: for (j = 0; j < ncols; j++) {
57: if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
58: }
59: }
60: MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals);
61: }
62: }
63: VecRestoreArray(y, &yarray);
64: VecRestoreArrayRead(x, &xarray);
65: return 0;
66: }
68: PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc, PetscOptionItems *PetscOptionsObject)
69: {
70: PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
72: PetscOptionsHeadBegin(PetscOptionsObject, "Kaczmarz options");
73: PetscOptionsReal("-pc_kaczmarz_lambda", "relaxation factor (0 < lambda)", "", jac->lambda, &jac->lambda, NULL);
74: PetscOptionsBool("-pc_kaczmarz_symmetric", "apply row projections symmetrically", "", jac->symmetric, &jac->symmetric, NULL);
75: PetscOptionsHeadEnd();
76: return 0;
77: }
79: PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer)
80: {
81: PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
82: PetscBool iascii;
84: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
85: if (iascii) PetscViewerASCIIPrintf(viewer, " lambda = %g\n", (double)jac->lambda);
86: return 0;
87: }
89: /*MC
90: PCKaczmarz - Kaczmarz iteration
92: Options Database Key:
93: . -pc_sor_lambda <1.0> - Sets damping parameter lambda
95: Level: beginner
97: Note:
98: In parallel this is block-Jacobi with Kaczmarz inner solve.
100: References:
101: . * - S. Kaczmarz, "Angenaherte Auflosing von Systemen Linearer Gleichungen",
102: Bull. Internat. Acad. Polon. Sci. C1. A, 1937.
104: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCBJACOBI`
105: M*/
107: PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
108: {
109: PC_Kaczmarz *jac;
111: PetscNew(&jac);
113: pc->ops->apply = PCApply_Kaczmarz;
114: pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz;
115: pc->ops->setup = NULL;
116: pc->ops->view = PCView_Kaczmarz;
117: pc->ops->destroy = PCDestroy_Kaczmarz;
118: pc->data = (void *)jac;
119: jac->lambda = 1.0;
120: jac->symmetric = PETSC_FALSE;
121: return 0;
122: }