1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/ 2 3 typedef struct { 4 PetscReal lambda; /* damping parameter */ 5 PetscBool symmetric; /* apply the projections symmetrically */ 6 } PC_Kaczmarz; 7 8 static PetscErrorCode PCDestroy_Kaczmarz(PC pc) 9 { 10 PetscFunctionBegin; 11 PetscCall(PetscFree(pc->data)); 12 PetscFunctionReturn(PETSC_SUCCESS); 13 } 14 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; 25 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 } 69 70 static PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc, PetscOptionItems PetscOptionsObject) 71 { 72 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data; 73 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 } 81 82 static PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer) 83 { 84 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data; 85 PetscBool isascii; 86 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 } 92 93 /*MC 94 PCKACZMARZ - Kaczmarz iteration {cite}`kaczmarz1937angenaherte` 95 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 99 100 Level: beginner 101 102 Note: 103 In parallel this is block-Jacobi with Kaczmarz inner solve on each processor. 104 105 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCJACOBI`, `PCBJACOBI` 106 M*/ 107 108 PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc) 109 { 110 PC_Kaczmarz *jac; 111 112 PetscFunctionBegin; 113 PetscCall(PetscNew(&jac)); 114 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 } 125