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(0); 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) { 39 r -= yarray[cols[j]-ys]*vals[j]; 40 } 41 anrm += PetscRealPart(PetscSqr(vals[j])); 42 } 43 if (anrm > 0.) { 44 for (j=0;j<ncols;j++) { 45 if (cols[j] >= ys && cols[j] < ye) { 46 yarray[cols[j]-ys] += vals[j]*lambda*r/anrm; 47 } 48 } 49 } 50 PetscCall(MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals)); 51 } 52 if (jac->symmetric) { 53 for (i=xe-1;i>=xs;i--) { 54 PetscCall(MatGetRow(pc->pmat,i,&ncols,&cols,&vals)); 55 r = xarray[i-xs]; 56 anrm = 0.; 57 for (j=0;j<ncols;j++) { 58 if (cols[j] >= ys && cols[j] < ye) { 59 r -= yarray[cols[j]-ys]*vals[j]; 60 } 61 anrm += PetscRealPart(PetscSqr(vals[j])); 62 } 63 if (anrm > 0.) { 64 for (j=0;j<ncols;j++) { 65 if (cols[j] >= ys && cols[j] < ye) { 66 yarray[cols[j]-ys] += vals[j]*lambda*r/anrm; 67 } 68 } 69 } 70 PetscCall(MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals)); 71 } 72 } 73 PetscCall(VecRestoreArray(y,&yarray)); 74 PetscCall(VecRestoreArrayRead(x,&xarray)); 75 PetscFunctionReturn(0); 76 } 77 78 PetscErrorCode PCSetFromOptions_Kaczmarz(PetscOptionItems *PetscOptionsObject,PC pc) 79 { 80 PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data; 81 82 PetscFunctionBegin; 83 PetscCall(PetscOptionsHead(PetscOptionsObject,"Kaczmarz options")); 84 PetscCall(PetscOptionsReal("-pc_kaczmarz_lambda","relaxation factor (0 < lambda)","",jac->lambda,&jac->lambda,NULL)); 85 PetscCall(PetscOptionsBool("-pc_kaczmarz_symmetric","apply row projections symmetrically","",jac->symmetric,&jac->symmetric,NULL)); 86 PetscCall(PetscOptionsTail()); 87 PetscFunctionReturn(0); 88 } 89 90 PetscErrorCode PCView_Kaczmarz(PC pc,PetscViewer viewer) 91 { 92 PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data; 93 PetscBool iascii; 94 95 PetscFunctionBegin; 96 PetscCall(PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii)); 97 if (iascii) { 98 PetscCall(PetscViewerASCIIPrintf(viewer," lambda = %g\n",(double)jac->lambda)); 99 } 100 PetscFunctionReturn(0); 101 } 102 103 /*MC 104 PCKaczmarz - Kaczmarz iteration 105 106 Options Database Keys: 107 . -pc_sor_lambda <1.0> - Sets damping parameter lambda 108 109 Level: beginner 110 111 Notes: 112 In parallel this is block-Jacobi with Kaczmarz inner solve. 113 114 References: 115 . * - S. Kaczmarz, "Angenaherte Auflosing von Systemen Linearer Gleichungen", 116 Bull. Internat. Acad. Polon. Sci. C1. A, 1937. 117 118 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 119 120 M*/ 121 122 PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc) 123 { 124 PC_Kaczmarz *jac; 125 126 PetscFunctionBegin; 127 PetscCall(PetscNewLog(pc,&jac)); 128 129 pc->ops->apply = PCApply_Kaczmarz; 130 pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz; 131 pc->ops->setup = NULL; 132 pc->ops->view = PCView_Kaczmarz; 133 pc->ops->destroy = PCDestroy_Kaczmarz; 134 pc->data = (void*)jac; 135 jac->lambda = 1.0; 136 jac->symmetric = PETSC_FALSE; 137 PetscFunctionReturn(0); 138 } 139