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 PetscFunctionBegin; 10 PetscCall(PetscFree(pc->data)); 11 PetscFunctionReturn(0); 12 } 13 14 static PetscErrorCode PCApply_Kaczmarz(PC pc, Vec x, Vec y) { 15 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data; 16 PetscInt xs, xe, ys, ye, ncols, i, j; 17 const PetscInt *cols; 18 const PetscScalar *vals, *xarray; 19 PetscScalar r; 20 PetscReal anrm; 21 PetscScalar *yarray; 22 PetscReal lambda = jac->lambda; 23 24 PetscFunctionBegin; 25 PetscCall(MatGetOwnershipRange(pc->pmat, &xs, &xe)); 26 PetscCall(MatGetOwnershipRangeColumn(pc->pmat, &ys, &ye)); 27 PetscCall(VecSet(y, 0.)); 28 PetscCall(VecGetArrayRead(x, &xarray)); 29 PetscCall(VecGetArray(y, &yarray)); 30 for (i = xs; i < xe; i++) { 31 /* get the maximum row width and row norms */ 32 PetscCall(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 PetscCall(MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals)); 45 } 46 if (jac->symmetric) { 47 for (i = xe - 1; i >= xs; i--) { 48 PetscCall(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 PetscCall(MatRestoreRow(pc->pmat, i, &ncols, &cols, &vals)); 61 } 62 } 63 PetscCall(VecRestoreArray(y, &yarray)); 64 PetscCall(VecRestoreArrayRead(x, &xarray)); 65 PetscFunctionReturn(0); 66 } 67 68 PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc, PetscOptionItems *PetscOptionsObject) { 69 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data; 70 71 PetscFunctionBegin; 72 PetscOptionsHeadBegin(PetscOptionsObject, "Kaczmarz options"); 73 PetscCall(PetscOptionsReal("-pc_kaczmarz_lambda", "relaxation factor (0 < lambda)", "", jac->lambda, &jac->lambda, NULL)); 74 PetscCall(PetscOptionsBool("-pc_kaczmarz_symmetric", "apply row projections symmetrically", "", jac->symmetric, &jac->symmetric, NULL)); 75 PetscOptionsHeadEnd(); 76 PetscFunctionReturn(0); 77 } 78 79 PetscErrorCode PCView_Kaczmarz(PC pc, PetscViewer viewer) { 80 PC_Kaczmarz *jac = (PC_Kaczmarz *)pc->data; 81 PetscBool iascii; 82 83 PetscFunctionBegin; 84 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii)); 85 if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, " lambda = %g\n", (double)jac->lambda)); } 86 PetscFunctionReturn(0); 87 } 88 89 /*MC 90 PCKaczmarz - Kaczmarz iteration 91 92 Options Database Keys: 93 . -pc_sor_lambda <1.0> - Sets damping parameter lambda 94 95 Level: beginner 96 97 Notes: 98 In parallel this is block-Jacobi with Kaczmarz inner solve. 99 100 References: 101 . * - S. Kaczmarz, "Angenaherte Auflosing von Systemen Linearer Gleichungen", 102 Bull. Internat. Acad. Polon. Sci. C1. A, 1937. 103 104 .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC` 105 106 M*/ 107 108 PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc) { 109 PC_Kaczmarz *jac; 110 111 PetscFunctionBegin; 112 PetscCall(PetscNewLog(pc, &jac)); 113 114 pc->ops->apply = PCApply_Kaczmarz; 115 pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz; 116 pc->ops->setup = NULL; 117 pc->ops->view = PCView_Kaczmarz; 118 pc->ops->destroy = PCDestroy_Kaczmarz; 119 pc->data = (void *)jac; 120 jac->lambda = 1.0; 121 jac->symmetric = PETSC_FALSE; 122 PetscFunctionReturn(0); 123 } 124