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