Permalink
Browse files

Merge branch 'prbrune/pckaczmarz'

* prbrune/pckaczmarz:
  PCKACZMARZ: Sign mistake for symmetric application
  PCKACZMARZ: Added symmetric application option.
  PCKACZMARZ: Added a basic Kaczmarz row-projection smoother as a PC

Conflicts:
	src/ksp/pc/impls/makefile
  • Loading branch information...
2 parents ed7d4a5 + d6d3918 commit 4c2751e9010174328ad406939f419b25d54f92a8 @prbrune prbrune committed Feb 19, 2014
View
@@ -78,6 +78,7 @@ typedef const char* PCType;
#define PCBICGSTABCUSP "bicgstabcusp"
#define PCAINVCUSP "ainvcusp"
#define PCBDDC "bddc"
+#define PCKACZMARZ "kaczmarz"
/* Logging support */
PETSC_EXTERN PetscClassId PC_CLASSID;
@@ -0,0 +1,155 @@
+#include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
+
+typedef struct {
+ PetscReal lambda; /* damping parameter */
+ PetscBool symmetric; /* apply the projections symmetrically */
+} PC_Kaczmarz;
+
+#undef __FUNCT__
+#define __FUNCT__ "PCDestroy_Kaczmarz"
+static PetscErrorCode PCDestroy_Kaczmarz(PC pc)
+{
+ PetscErrorCode ierr;
+
+ PetscFunctionBegin;
+ ierr = PetscFree(pc->data);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PCApply_Kaczmarz"
+static PetscErrorCode PCApply_Kaczmarz(PC pc,Vec x,Vec y)
+{
+ PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
+ PetscInt xs,xe,ys,ye,ncols,i,j;
+ const PetscInt *cols;
+ const PetscScalar *vals;
+ PetscErrorCode ierr;
+ PetscScalar r;
+ PetscReal anrm;
+ PetscScalar *xarray,*yarray;
+ PetscReal lambda=jac->lambda;
+
+ PetscFunctionBegin;
+ ierr = MatGetOwnershipRange(pc->pmat,&xs,&xe);CHKERRQ(ierr);
+ ierr = MatGetOwnershipRangeColumn(pc->pmat,&ys,&ye);CHKERRQ(ierr);
+ ierr = VecSet(y,0.);CHKERRQ(ierr);
+ ierr = VecGetArray(x,&xarray);CHKERRQ(ierr);
+ ierr = VecGetArray(y,&yarray);CHKERRQ(ierr);
+ for (i=xs;i<xe;i++) {
+ /* get the maximum row width and row norms */
+ ierr = MatGetRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
+ r = xarray[i-xs];
+ anrm = 0.;
+ for (j=0;j<ncols;j++) {
+ if (cols[j] >= ys && cols[j] < ye) {
+ r -= yarray[cols[j]-ys]*vals[j];
+ }
+ anrm += PetscRealPart(PetscSqr(vals[j]));
+ }
+ if (anrm > 0.) {
+ for (j=0;j<ncols;j++) {
+ if (cols[j] >= ys && cols[j] < ye) {
+ yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
+ }
+ }
+ }
+ ierr = MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
+ }
+ if (jac->symmetric) {
+ for (i=xe-1;i>=xs;i--) {
+ ierr = MatGetRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
+ r = xarray[i-xs];
+ anrm = 0.;
+ for (j=0;j<ncols;j++) {
+ if (cols[j] >= ys && cols[j] < ye) {
+ r -= yarray[cols[j]-ys]*vals[j];
+ }
+ anrm += PetscRealPart(PetscSqr(vals[j]));
+ }
+ if (anrm > 0.) {
+ for (j=0;j<ncols;j++) {
+ if (cols[j] >= ys && cols[j] < ye) {
+ yarray[cols[j]-ys] += vals[j]*lambda*r/anrm;
+ }
+ }
+ }
+ ierr = MatRestoreRow(pc->pmat,i,&ncols,&cols,&vals);CHKERRQ(ierr);
+ }
+ }
+ ierr = VecRestoreArray(y,&yarray);CHKERRQ(ierr);
+ ierr = VecRestoreArray(x,&xarray);CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PCSetFromOptions_Kaczmarz"
+PetscErrorCode PCSetFromOptions_Kaczmarz(PC pc)
+{
+ PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
+ PetscErrorCode ierr;
+
+ PetscFunctionBegin;
+ ierr = PetscOptionsHead("Kaczmarz options");CHKERRQ(ierr);
+ ierr = PetscOptionsReal("-pc_kaczmarz_lambda","relaxation factor (0 < lambda)","",jac->lambda,&jac->lambda,0);CHKERRQ(ierr);
+ ierr = PetscOptionsBool("-pc_kaczmarz_symmetric","apply row projections symmetrically","",jac->symmetric,&jac->symmetric,0);CHKERRQ(ierr);
+ ierr = PetscOptionsTail();CHKERRQ(ierr);
+ PetscFunctionReturn(0);
+}
+
+#undef __FUNCT__
+#define __FUNCT__ "PCView_Kaczmarz"
+PetscErrorCode PCView_Kaczmarz(PC pc,PetscViewer viewer)
+{
+ PC_Kaczmarz *jac = (PC_Kaczmarz*)pc->data;
+ PetscErrorCode ierr;
+ PetscBool iascii;
+
+ PetscFunctionBegin;
+ ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
+ if (iascii) {
+ ierr = PetscViewerASCIIPrintf(viewer," Kaczmarz: lambda = %G\n",jac->lambda);CHKERRQ(ierr);
+ }
+ PetscFunctionReturn(0);
+}
+
+/*MC
+ PCKaczmarz - Kaczmarz iteration
+
+ Options Database Keys:
+. -pc_sor_lambda <1.0> - Sets damping parameter lambda
+
+ Level: beginner
+
+ Concepts: Kaczmarz, preconditioners, row projection
+
+ Notes: In parallel this is block-Jacobi with Kaczmarz inner solve.
+
+ References:
+ S. Kaczmarz, “Angenaherte Auflosing von Systemen Linearer Gleichungen”,
+ Bull. Internat. Acad. Polon. Sci. C1. A, pp.335-357, 1937.
+
+.seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC
+
+M*/
+
+#undef __FUNCT__
+#define __FUNCT__ "PCCreate_Kaczmarz"
+PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC pc)
+{
+ PetscErrorCode ierr;
+ PC_Kaczmarz *jac;
+
+ PetscFunctionBegin;
+ ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
+
+ pc->ops->apply = PCApply_Kaczmarz;
+ pc->ops->setfromoptions = PCSetFromOptions_Kaczmarz;
+ pc->ops->setup = 0;
+ pc->ops->view = PCView_Kaczmarz;
+ pc->ops->destroy = PCDestroy_Kaczmarz;
+ pc->data = (void*)jac;
+ jac->lambda = 1.0;
+ jac->symmetric = PETSC_FALSE;
+ PetscFunctionReturn(0);
+}
@@ -0,0 +1,15 @@
+
+ALL: lib
+
+CFLAGS = ${GAMG_INCLUDE}
+FFLAGS =
+SOURCEC = kaczmarz.c
+SOURCEF =
+SOURCEH =
+LIBBASE = libpetscksp
+MANSEC = PC
+LOCDIR = src/ksp/pc/impls/kaczmarz/
+
+include ${PETSC_DIR}/conf/variables
+include ${PETSC_DIR}/conf/rules
+include ${PETSC_DIR}/conf/test
@@ -4,7 +4,7 @@ ALL: lib
LIBBASE = libpetscksp
DIRS = jacobi none sor shell bjacobi mg eisens asm ksp composite redundant spai is pbjacobi ml\
mat hypre tfs fieldsplit factor galerkin cp wb python ainvcusp sacusp bicgstabcusp\
- lsc redistribute gasm svd gamg parms bddc
+ lsc redistribute gasm svd gamg parms bddc kaczmarz
LOCDIR = src/ksp/pc/impls/
include ${PETSC_DIR}/conf/variables
@@ -27,6 +27,7 @@ PETSC_EXTERN PetscErrorCode PCCreate_LSC(PC);
PETSC_EXTERN PetscErrorCode PCCreate_Redistribute(PC);
PETSC_EXTERN PetscErrorCode PCCreate_SVD(PC);
PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC);
+PETSC_EXTERN PetscErrorCode PCCreate_Kaczmarz(PC);
#if defined(PETSC_HAVE_ML)
PETSC_EXTERN PetscErrorCode PCCreate_ML(PC);
@@ -106,6 +107,7 @@ PetscErrorCode PCRegisterAll(void)
ierr = PCRegister(PCREDISTRIBUTE ,PCCreate_Redistribute);CHKERRQ(ierr);
ierr = PCRegister(PCSVD ,PCCreate_SVD);CHKERRQ(ierr);
ierr = PCRegister(PCGAMG ,PCCreate_GAMG);CHKERRQ(ierr);
+ ierr = PCRegister(PCKACZMARZ ,PCCreate_Kaczmarz);CHKERRQ(ierr);
#if defined(PETSC_HAVE_ML)
ierr = PCRegister(PCML ,PCCreate_ML);CHKERRQ(ierr);
#endif

0 comments on commit 4c2751e

Please sign in to comment.