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:   PetscFunctionBegin;
 11:   PetscCall(PetscFree(pc->data));
 12:   PetscFunctionReturn(PETSC_SUCCESS);
 13: }

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

 26:   PetscFunctionBegin;
 27:   PetscCall(MatGetOwnershipRange(pc->pmat, &xs, &xe));
 28:   PetscCall(MatGetOwnershipRangeColumn(pc->pmat, &ys, &ye));
 29:   PetscCall(VecSet(y, 0.));
 30:   PetscCall(VecGetArrayRead(x, &xarray));
 31:   PetscCall(VecGetArray(y, &yarray));
 32:   for (i = xs; i < xe; i++) {
 33:     /* get the maximum row width and row norms */
 34:     PetscCall(MatGetRow(pc->pmat, i, &ncols, &cols, &vals));
 35:     r    = xarray[i - xs];
 36:     anrm = 0.;
 37:     for (j = 0; j < ncols; j++) {
 38:       if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
 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) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
 44:       }
 45:     }
 46:     PetscCall(MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals));
 47:   }
 48:   if (jac->symmetric) {
 49:     for (i = xe - 1; i >= xs; i--) {
 50:       PetscCall(MatGetRow(pc->pmat, i, &ncols, &cols, &vals));
 51:       r    = xarray[i - xs];
 52:       anrm = 0.;
 53:       for (j = 0; j < ncols; j++) {
 54:         if (cols[j] >= ys && cols[j] < ye) r -= yarray[cols[j] - ys] * vals[j];
 55:         anrm += PetscRealPart(PetscSqr(vals[j]));
 56:       }
 57:       if (anrm > 0.) {
 58:         for (j = 0; j < ncols; j++) {
 59:           if (cols[j] >= ys && cols[j] < ye) yarray[cols[j] - ys] += vals[j] * lambda * r / anrm;
 60:         }
 61:       }
 62:       PetscCall(MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals));
 63:     }
 64:   }
 65:   PetscCall(VecRestoreArray(y, &yarray));
 66:   PetscCall(VecRestoreArrayRead(x, &xarray));
 67:   PetscFunctionReturn(PETSC_SUCCESS);
 68: }

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

 74:   PetscFunctionBegin;
 75:   PetscOptionsHeadBegin(PetscOptionsObject, "Kaczmarz options");
 76:   PetscCall(PetscOptionsReal("-pc_kaczmarz_lambda", "relaxation factor (0 < lambda)", "", jac->lambda, &jac->lambda, NULL));
 77:   PetscCall(PetscOptionsBool("-pc_kaczmarz_symmetric", "apply row projections symmetrically", "", jac->symmetric, &jac->symmetric, NULL));
 78:   PetscOptionsHeadEnd();
 79:   PetscFunctionReturn(PETSC_SUCCESS);
 80: }

 82: static PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer)
 83: {
 84:   PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data;
 85:   PetscBool    iascii;

 87:   PetscFunctionBegin;
 88:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
 89:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  lambda = %g\n", (double)jac->lambda));
 90:   PetscFunctionReturn(PETSC_SUCCESS);
 91: }

 93: /*MC
 94:      PCKACZMARZ - Kaczmarz iteration {cite}`kaczmarz1937angenaherte`

 96:    Options Database Keys:
 97: +  -pc_kaczmarz_lambda <lambda> - Sets damping parameter defaults to 1.0
 98: -  -pc_kaczmarz_symmetric       - Apply the row projections symmetrically

100:    Level: beginner

102:    Note:
103:    In parallel this is block-Jacobi with Kaczmarz inner solve on each processor.

105: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCBJACOBI`
106: M*/

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

112:   PetscFunctionBegin;
113:   PetscCall(PetscNew(&jac));

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