Actual source code: kaczmarz.c
petsc-3.6.1 2015-08-06
1: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
3: typedef struct {
4: PetscReal lambda; /* damping parameter */
5: PetscBool symmetric; /* apply the projections symmetrically */
6: } PC_Kaczmarz;
10: static PetscErrorCode PCDestroy_Kaczmarz(PC pc)
11: {
15: PetscFree(pc->data);
16: return(0);
17: }
21: static PetscErrorCode PCApply_Kaczmarz(PC pc,Vec x,Vec y)
22: {
23: PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
24: PetscInt xs,xe,ys,ye,ncols,i,j;
25: const PetscInt *cols;
26: const PetscScalar *vals;
27: PetscErrorCode ierr;
28: PetscScalar r;
29: PetscReal anrm;
30: PetscScalar *xarray,*yarray;
31: PetscReal lambda=jac->lambda;
34: MatGetOwnershipRange(pc->pmat,&xs,&xe);
35: MatGetOwnershipRangeColumn(pc->pmat,&ys,&ye);
36: VecSet(y,0.);
37: VecGetArray(x,&xarray);
38: VecGetArray(y,&yarray);
39: for (i=xs;i<xe;i++) {
40: /* get the maximum row width and row norms */
41: MatGetRow(pc->pmat,i,&ncols,&cols,&vals);
42: r = xarray[i-xs];
43: anrm = 0.;
44: for (j=0;j<ncols;j++) {
45: if (cols[j] >= ys && cols[j] < ye) {
46: r -= yarray[cols[j]-ys]*vals[j];
47: }
48: anrm += PetscRealPart(PetscSqr(vals[j]));
49: }
50: if (anrm > 0.) {
51: for (j=0;j<ncols;j++) {
52: if (cols[j] >= ys && cols[j] < ye) {
53: yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
54: }
55: }
56: }
57: MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);
58: }
59: if (jac->symmetric) {
60: for (i=xe-1;i>=xs;i--) {
61: MatGetRow(pc->pmat,i,&ncols,&cols,&vals);
62: r = xarray[i-xs];
63: anrm = 0.;
64: for (j=0;j<ncols;j++) {
65: if (cols[j] >= ys && cols[j] < ye) {
66: r -= yarray[cols[j]-ys]*vals[j];
67: }
68: anrm += PetscRealPart(PetscSqr(vals[j]));
69: }
70: if (anrm > 0.) {
71: for (j=0;j<ncols;j++) {
72: if (cols[j] >= ys && cols[j] < ye) {
73: yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
74: }
75: }
76: }
77: MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);
78: }
79: }
80: VecRestoreArray(y,&yarray);
81: VecRestoreArray(x,&xarray);
82: return(0);
83: }
87: PetscErrorCode PCSetFromOptions_Kaczmarz(PetscOptions *PetscOptionsObject,PC pc)
88: {
89: PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
93: PetscOptionsHead(PetscOptionsObject,"Kaczmarz options");
94: PetscOptionsReal("-pc_kaczmarz_lambda","relaxation factor (0 < lambda)","",jac->lambda,&jac->lambda,NULL);
95: PetscOptionsBool("-pc_kaczmarz_symmetric","apply row projections symmetrically","",jac->symmetric,&jac->symmetric,NULL);
96: PetscOptionsTail();
97: return(0);
98: }
102: PetscErrorCode PCView_Kaczmarz(PC pc,PetscViewer viewer)
103: {
104: PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
106: PetscBool iascii;
109: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
110: if (iascii) {
111: PetscViewerASCIIPrintf(viewer," Kaczmarz: lambda = %G\n",jac->lambda);
112: }
113: return(0);
114: }
116: /*MC
117: PCKaczmarz - Kaczmarz iteration
119: Options Database Keys:
120: . -pc_sor_lambda <1.0> - Sets damping parameter lambda
122: Level: beginner
124: Concepts: Kaczmarz, preconditioners, row projection
126: Notes: In parallel this is block-Jacobi with Kaczmarz inner solve.
128: References:
129: S. Kaczmarz, �Angenaherte Auflosing von Systemen Linearer Gleichungen�,
130: Bull. Internat. Acad. Polon. Sci. C1. A, pp.335-357, 1937.
132: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
134: M*/
138: PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
139: {
141: PC_Kaczmarz *jac;
144: PetscNewLog(pc,&jac);
146: pc->ops->apply = PCApply_Kaczmarz;
147: pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz;
148: pc->ops->setup = 0;
149: pc->ops->view = PCView_Kaczmarz;
150: pc->ops->destroy = PCDestroy_Kaczmarz;
151: pc->data = (void*)jac;
152: jac->lambda = 1.0;
153: jac->symmetric = PETSC_FALSE;
154: return(0);
155: }