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) {
37: r -= yarray[cols[j]-ys]*vals[j];
38: }
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) {
44: yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
45: }
46: }
47: }
48: MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);
49: }
50: if (jac->symmetric) {
51: for (i=xe-1;i>=xs;i--) {
52: MatGetRow(pc->pmat,i,&ncols,&cols,&vals);
53: r = xarray[i-xs];
54: anrm = 0.;
55: for (j=0;j<ncols;j++) {
56: if (cols[j] >= ys && cols[j] < ye) {
57: r -= yarray[cols[j]-ys]*vals[j];
58: }
59: anrm += PetscRealPart(PetscSqr(vals[j]));
60: }
61: if (anrm > 0.) {
62: for (j=0;j<ncols;j++) {
63: if (cols[j] >= ys && cols[j] < ye) {
64: yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
65: }
66: }
67: }
68: MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);
69: }
70: }
71: VecRestoreArray(y,&yarray);
72: VecRestoreArrayRead(x,&xarray);
73: return 0;
74: }
76: PetscErrorCode PCSetFromOptions_Kaczmarz(PetscOptionItems *PetscOptionsObject,PC pc)
77: {
78: PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
80: PetscOptionsHead(PetscOptionsObject,"Kaczmarz options");
81: PetscOptionsReal("-pc_kaczmarz_lambda","relaxation factor (0 < lambda)","",jac->lambda,&jac->lambda,NULL);
82: PetscOptionsBool("-pc_kaczmarz_symmetric","apply row projections symmetrically","",jac->symmetric,&jac->symmetric,NULL);
83: PetscOptionsTail();
84: return 0;
85: }
87: PetscErrorCode PCView_Kaczmarz(PC pc,PetscViewer viewer)
88: {
89: PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
90: PetscBool iascii;
92: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
93: if (iascii) {
94: PetscViewerASCIIPrintf(viewer," lambda = %g\n",(double)jac->lambda);
95: }
96: return 0;
97: }
99: /*MC
100: PCKaczmarz - Kaczmarz iteration
102: Options Database Keys:
103: . -pc_sor_lambda <1.0> - Sets damping parameter lambda
105: Level: beginner
107: Notes:
108: In parallel this is block-Jacobi with Kaczmarz inner solve.
110: References:
111: . * - S. Kaczmarz, "Angenaherte Auflosing von Systemen Linearer Gleichungen",
112: Bull. Internat. Acad. Polon. Sci. C1. A, 1937.
114: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
116: M*/
118: PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
119: {
120: PC_Kaczmarz *jac;
122: PetscNewLog(pc,&jac);
124: pc->ops->apply = PCApply_Kaczmarz;
125: pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz;
126: pc->ops->setup = NULL;
127: pc->ops->view = PCView_Kaczmarz;
128: pc->ops->destroy = PCDestroy_Kaczmarz;
129: pc->data = (void*)jac;
130: jac->lambda = 1.0;
131: jac->symmetric = PETSC_FALSE;
132: return 0;
133: }