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