1 #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/ 2 3 typedef struct { 4 PetscReal lambda; /* damping parameter */ 5 } PC_Kaczmarz; 6 7 #undef __FUNCT__ 8 #define __FUNCT__ "PCDestroy_Kaczmarz" 9 static PetscErrorCode PCDestroy_Kaczmarz(PC pc) 10 { 11 PetscErrorCode ierr; 12 13 PetscFunctionBegin; 14 ierr = PetscFree(pc->data);CHKERRQ(ierr); 15 PetscFunctionReturn(0); 16 } 17 18 #undef __FUNCT__ 19 #define __FUNCT__ "PCApply_Kaczmarz" 20 static PetscErrorCode PCApply_Kaczmarz(PC pc,Vec x,Vec y) 21 { 22 PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data; 23 PetscInt xs,xe,ys,ye,ncols,i,j; 24 const PetscInt *cols; 25 const PetscScalar *vals; 26 PetscErrorCode ierr; 27 PetscScalar r; 28 PetscReal anrm; 29 PetscScalar *xarray,*yarray; 30 PetscReal lambda=jac->lambda; 31 32 PetscFunctionBegin; 33 ierr = MatGetOwnershipRange(pc->pmat,&xs,&xe);CHKERRQ(ierr); 34 ierr = MatGetOwnershipRangeColumn(pc->pmat,&ys,&ye);CHKERRQ(ierr); 35 ierr = VecSet(y,0);CHKERRQ(ierr); 36 ierr = VecGetArray(x,&xarray);CHKERRQ(ierr); 37 ierr = VecGetArray(y,&yarray);CHKERRQ(ierr); 38 for (i=xs;i<xe;i++) { 39 /* get the maximum row width and row norms */ 40 ierr = MatGetRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr); 41 /* form the full residual */ 42 r = xarray[i-xs]; 43 anrm = 0.; 44 for (j=0;j<ncols;j++) { 45 if (cols[j] >= ys && cols[j] < ye) { 46 r -= yarray[cols[j]-ys]*vals[j]; 47 } 48 anrm += PetscRealPart(PetscPowScalar(vals[j],2)); 49 } 50 if (anrm > 0) { 51 for (j=0;j<ncols;j++) { 52 if (cols[j] >= ys && cols[j] < ye) { 53 yarray[cols[j]-ys] += vals[j]*lambda*r/anrm; 54 } 55 } 56 } 57 ierr = MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr); 58 } 59 ierr = VecRestoreArray(y,&yarray);CHKERRQ(ierr); 60 ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr); 61 PetscFunctionReturn(0); 62 } 63 64 #undef __FUNCT__ 65 #define __FUNCT__ "PCSetFromOptions_Kaczmarz" 66 PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc) 67 { 68 PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data; 69 PetscErrorCode ierr; 70 71 PetscFunctionBegin; 72 ierr = PetscOptionsHead("Kaczmarz options");CHKERRQ(ierr); 73 ierr = PetscOptionsReal("-pc_kaczmarz_lambda","relaxation factor (0 < lambda)","",jac->lambda,&jac->lambda,0);CHKERRQ(ierr); 74 ierr = PetscOptionsTail();CHKERRQ(ierr); 75 PetscFunctionReturn(0); 76 } 77 78 #undef __FUNCT__ 79 #define __FUNCT__ "PCView_Kaczmarz" 80 PetscErrorCode PCView_Kaczmarz(PC pc,PetscViewer viewer) 81 { 82 PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data; 83 PetscErrorCode ierr; 84 PetscBool iascii; 85 86 PetscFunctionBegin; 87 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr); 88 if (iascii) { 89 ierr = PetscViewerASCIIPrintf(viewer," Kaczmarz: lambda = %G\n",jac->lambda);CHKERRQ(ierr); 90 } 91 PetscFunctionReturn(0); 92 } 93 94 /*MC 95 PCKaczmarz - Kaczmarz iteration 96 97 Options Database Keys: 98 . -pc_sor_lambda <1.0> - Sets damping parameter lambda 99 100 Level: beginner 101 102 Concepts: Kaczmarz, preconditioners, row projection 103 104 Notes: In parallel this computes off-processor updates at the end of every local sweep 105 106 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC 107 108 M*/ 109 110 #undef __FUNCT__ 111 #define __FUNCT__ "PCCreate_Kaczmarz" 112 PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc) 113 { 114 PetscErrorCode ierr; 115 PC_Kaczmarz *jac; 116 117 PetscFunctionBegin; 118 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr); 119 120 pc->ops->apply = PCApply_Kaczmarz; 121 pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz; 122 pc->ops->setup = 0; 123 pc->ops->view = PCView_Kaczmarz; 124 pc->ops->destroy = PCDestroy_Kaczmarz; 125 pc->data = (void*)jac; 126 jac->lambda = 1.0; 127 PetscFunctionReturn(0); 128 } 129