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(PetscOptionItems *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: .  1. - S. Kaczmarz, “Angenaherte Auflosing von Systemen Linearer Gleichungen”,
130:    Bull. Internat. Acad. Polon. Sci. C1. A, 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: }