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: {

 13:   PetscFree(pc->data);
 14:   return(0);
 15: }

 17: static PetscErrorCode PCApply_Kaczmarz(PC pc,Vec x,Vec y)
 18: {
 19:   PC_Kaczmarz       *jac = (PC_Kaczmarz*)pc->data;
 20:   PetscInt          xs,xe,ys,ye,ncols,i,j;
 21:   const PetscInt    *cols;
 22:   const PetscScalar *vals,*xarray;
 23:   PetscErrorCode    ierr;
 24:   PetscScalar       r;
 25:   PetscReal         anrm;
 26:   PetscScalar       *yarray;
 27:   PetscReal         lambda=jac->lambda;

 30:   MatGetOwnershipRange(pc->pmat,&xs,&xe);
 31:   MatGetOwnershipRangeColumn(pc->pmat,&ys,&ye);
 32:   VecSet(y,0.);
 33:   VecGetArrayRead(x,&xarray);
 34:   VecGetArray(y,&yarray);
 35:   for (i=xs;i<xe;i++) {
 36:     /* get the maximum row width and row norms */
 37:     MatGetRow(pc->pmat,i,&ncols,&cols,&vals);
 38:     r = xarray[i-xs];
 39:     anrm = 0.;
 40:     for (j=0;j<ncols;j++) {
 41:       if (cols[j] >= ys && cols[j] < ye) {
 42:         r -= yarray[cols[j]-ys]*vals[j];
 43:       }
 44:       anrm += PetscRealPart(PetscSqr(vals[j]));
 45:     }
 46:     if (anrm > 0.) {
 47:       for (j=0;j<ncols;j++) {
 48:         if (cols[j] >= ys && cols[j] < ye) {
 49:           yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
 50:         }
 51:       }
 52:     }
 53:     MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);
 54:   }
 55:   if (jac->symmetric) {
 56:     for (i=xe-1;i>=xs;i--) {
 57:       MatGetRow(pc->pmat,i,&ncols,&cols,&vals);
 58:       r = xarray[i-xs];
 59:       anrm = 0.;
 60:       for (j=0;j<ncols;j++) {
 61:         if (cols[j] >= ys && cols[j] < ye) {
 62:           r -= yarray[cols[j]-ys]*vals[j];
 63:         }
 64:         anrm += PetscRealPart(PetscSqr(vals[j]));
 65:       }
 66:       if (anrm > 0.) {
 67:         for (j=0;j<ncols;j++) {
 68:           if (cols[j] >= ys && cols[j] < ye) {
 69:             yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
 70:           }
 71:         }
 72:       }
 73:       MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);
 74:     }
 75:   }
 76:   VecRestoreArray(y,&yarray);
 77:   VecRestoreArrayRead(x,&xarray);
 78:   return(0);
 79: }

 81: PetscErrorCode PCSetFromOptions_Kaczmarz(PetscOptionItems *PetscOptionsObject,PC pc)
 82: {
 83:   PC_Kaczmarz    *jac = (PC_Kaczmarz*)pc->data;

 87:   PetscOptionsHead(PetscOptionsObject,"Kaczmarz options");
 88:   PetscOptionsReal("-pc_kaczmarz_lambda","relaxation factor (0 < lambda)","",jac->lambda,&jac->lambda,NULL);
 89:   PetscOptionsBool("-pc_kaczmarz_symmetric","apply row projections symmetrically","",jac->symmetric,&jac->symmetric,NULL);
 90:   PetscOptionsTail();
 91:   return(0);
 92: }

 94: PetscErrorCode PCView_Kaczmarz(PC pc,PetscViewer viewer)
 95: {
 96:   PC_Kaczmarz    *jac = (PC_Kaczmarz*)pc->data;
 98:   PetscBool      iascii;

101:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
102:   if (iascii) {
103:     PetscViewerASCIIPrintf(viewer,"  lambda = %g\n",(double)jac->lambda);
104:   }
105:   return(0);
106: }

108: /*MC
109:      PCKaczmarz - Kaczmarz iteration

111:    Options Database Keys:
112: .  -pc_sor_lambda <1.0> - Sets damping parameter lambda

114:    Level: beginner

116:    Notes:
117:     In parallel this is block-Jacobi with Kaczmarz inner solve.

119:    References:
120: .  1. - S. Kaczmarz, "Angenaherte Auflosing von Systemen Linearer Gleichungen",
121:    Bull. Internat. Acad. Polon. Sci. C1. A, 1937.

123: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC

125: M*/

127: PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
128: {
130:   PC_Kaczmarz    *jac;

133:   PetscNewLog(pc,&jac);

135:   pc->ops->apply           = PCApply_Kaczmarz;
136:   pc->ops->setfromoptions  = PCSetFromOptions_Kaczmarz;
137:   pc->ops->setup           = NULL;
138:   pc->ops->view            = PCView_Kaczmarz;
139:   pc->ops->destroy         = PCDestroy_Kaczmarz;
140:   pc->data                 = (void*)jac;
141:   jac->lambda              = 1.0;
142:   jac->symmetric           = PETSC_FALSE;
143:   return(0);
144: }