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    isascii;

     87:   PetscFunctionBegin;
     88:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
     89:   if (isascii) 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: }