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) r -= yarray[cols[j] - ys] * vals[j];
 37:       anrm += PetscRealPart(PetscSqr(vals[j]));
 38:     }
 39:     if (anrm > 0.) {
 40:       for (j = 0; j < ncols; j++) {
 41:         if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
 42:       }
 43:     }
 44:     MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals);
 45:   }
 46:   if (jac->symmetric) {
 47:     for (i = xe - 1; i >= xs; i--) {
 48:       MatGetRow(pc->pmat, i, &ncols, &cols, &vals);
 49:       r    = xarray[i - xs];
 50:       anrm = 0.;
 51:       for (j = 0; j < ncols; j++) {
 52:         if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
 53:         anrm += PetscRealPart(PetscSqr(vals[j]));
 54:       }
 55:       if (anrm > 0.) {
 56:         for (j = 0; j < ncols; j++) {
 57:           if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
 58:         }
 59:       }
 60:       MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals);
 61:     }
 62:   }
 63:   VecRestoreArray(y, &yarray);
 64:   VecRestoreArrayRead(x, &xarray);
 65:   return 0;
 66: }

 68: PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc, PetscOptionItems *PetscOptionsObject)
 69: {
 70:   PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;

 72:   PetscOptionsHeadBegin(PetscOptionsObject, "Kaczmarz options");
 73:   PetscOptionsReal("-pc_kaczmarz_lambda", "relaxation factor (0 < lambda)", "", jac->lambda, &jac->lambda, NULL);
 74:   PetscOptionsBool("-pc_kaczmarz_symmetric", "apply row projections symmetrically", "", jac->symmetric, &jac->symmetric, NULL);
 75:   PetscOptionsHeadEnd();
 76:   return 0;
 77: }

 79: PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer)
 80: {
 81:   PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
 82:   PetscBool    iascii;

 84:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
 85:   if (iascii) PetscViewerASCIIPrintf(viewer, "  lambda = %g\n", (double)jac->lambda);
 86:   return 0;
 87: }

 89: /*MC
 90:      PCKaczmarz - Kaczmarz iteration

 92:    Options Database Key:
 93: .  -pc_sor_lambda <1.0> - Sets damping parameter lambda

 95:    Level: beginner

 97:    Note:
 98:     In parallel this is block-Jacobi with Kaczmarz inner solve.

100:    References:
101: .  * - S. Kaczmarz, "Angenaherte Auflosing von Systemen Linearer Gleichungen",
102:    Bull. Internat. Acad. Polon. Sci. C1. A, 1937.

104: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCBJACOBI`
105: M*/

107: PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
108: {
109:   PC_Kaczmarz *jac;

111:   PetscNew(&jac);

113:   pc->ops->apply          = PCApply_Kaczmarz;
114:   pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz;
115:   pc->ops->setup          = NULL;
116:   pc->ops->view           = PCView_Kaczmarz;
117:   pc->ops->destroy        = PCDestroy_Kaczmarz;
118:   pc->data                = (void *)jac;
119:   jac->lambda             = 1.0;
120:   jac->symmetric          = PETSC_FALSE;
121:   return 0;
122: }