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