xref: /petsc/src/ksp/pc/impls/kaczmarz/kaczmarz.c (revision 3f93e5bd7810365176b4de538c1ab7498b758201)

#include <petsc-private/pcimpl.h>               /*I "petscpc.h" I*/

typedef struct {
  PetscReal  lambda; /* damping parameter */
} 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);
    /* form the full residual */
    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(PetscPowScalar(vals[j],2));
    }
    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 = 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 computes off-processor updates at the end of every local sweep

.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;
  PetscFunctionReturn(0);
}