[petsc-users] How does KSPSetNullSpace() work?

[Barry Smith]

On Jan 10, 2012, at 7:43 AM, Jed Brown wrote:

> On Tue, Jan 10, 2012 at 03:50, Alexander Grayver <agrayver at gfz-potsdam.de> wrote:
> http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/KSP/KSPSetNullSpace.html
> 
> You set nullspace without any particular information:
> 
> ierr = MatNullSpaceCreate(PETSC_COMM_WORLD, PETSC_TRUE, 0, PETSC_NULL, &nullsp);
> ierr = KSPSetNullSpace(ksp, nullsp);
> 
> What does it project and how works in this case?
> 
> This removes the constant null space. The has_const argument is equivalent to creating a constant vector (but more efficient).

   Alexandar,

      You need to learn to use etags or one of the other mechanisms for searching PETSc source code, see sections 

13.8 Emacs Users . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174
 13.9 Vi and Vim Users . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174 
13.10EclipseUsers ..........................................175 
13.11Qt Creator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175 13.12
DevelopersStudioUsers ....................................177 13.13
XCodeUsers(TheAppleGUIDevelopmentSystem . . . . . . . . . . . . . . . . . . . . . 177

in the users manual: http://www.mcs.anl.gov/petsc/petsc-dev/docs/manual.pdf

This way you can see for yourself exactly how PETSc is doing things. It will be quicker and more accurate than always asking us (especially since we have to look at the code to answer your question anyways). Anyways

Here is how the multiplies and preconditioner applies are done in the KSP solvers

#define KSP_MatMult(ksp,A,x,y)          (!ksp->transpose_solve) ? MatMult(A,x,y)                                                            : MatMultTranspose(A,x,y) 
#define KSP_MatMultTranspose(ksp,A,x,y) (!ksp->transpose_solve) ? MatMultTranspose(A,x,y)                                                   : MatMult(A,x,y) 
#define KSP_PCApply(ksp,x,y)            (!ksp->transpose_solve) ? (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y))                      : PCApplyTranspose(ksp->pc,x,y) 
#define KSP_PCApplyTranspose(ksp,x,y)   (!ksp->transpose_solve) ? PCApplyTranspose(ksp->pc,x,y)                                             : (PCApply(ksp->pc,x,y) || KSP_RemoveNullSpace(ksp,y)) 
#define KSP_PCApplyBAorAB(ksp,x,y,w)    (!ksp->transpose_solve) ? (PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w) || KSP_RemoveNullSpace(ksp,y)) : PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w)
#define KSP_PCApplyBAorABTranspose(ksp,x,y,w)    (!ksp->transpose_solve) ? (PCApplyBAorABTranspose(ksp->pc,ksp->pc_side,x,y,w) || KSP_RemoveNullSpace(ksp,y)) : PCApplyBAorAB(ksp->pc,ksp->pc_side,x,y,w)

while 

#define KSP_RemoveNullSpace(ksp,y) ((ksp->nullsp && ksp->pc_side == PC_LEFT) ? MatNullSpaceRemove(ksp->nullsp,y,PETSC_NULL) : 0)

and 

/*@C
   MatNullSpaceRemove - Removes all the components of a null space from a vector.

   Collective on MatNullSpace

   Input Parameters:
+  sp - the null space context
.  vec - the vector from which the null space is to be removed 
-  out - if this is requested (not PETSC_NULL) then this is a vector with the null space removed otherwise
         the removal is done in-place (in vec)

   Note: The user is not responsible for the vector returned and should not destroy it.

   Level: advanced

.keywords: PC, null space, remove

.seealso: MatNullSpaceCreate(), MatNullSpaceDestroy(), MatNullSpaceSetFunction()
@*/
PetscErrorCode  MatNullSpaceRemove(MatNullSpace sp,Vec vec,Vec *out)
{
  PetscScalar    sum;
  PetscInt       i,N;
  PetscErrorCode ierr;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(sp,MAT_NULLSPACE_CLASSID,1); 
  PetscValidHeaderSpecific(vec,VEC_CLASSID,2); 

  if (out) {
    PetscValidPointer(out,3); 
    if (!sp->vec) { 
      ierr = VecDuplicate(vec,&sp->vec);CHKERRQ(ierr); 
      ierr = PetscLogObjectParent(sp,sp->vec);CHKERRQ(ierr);
    }
    ierr = VecCopy(vec,sp->vec);CHKERRQ(ierr);
    vec = *out = sp->vec;
  }
  
  if (sp->has_cnst) {
    ierr = VecGetSize(vec,&N);CHKERRQ(ierr);
    if (N > 0) {
      ierr = VecSum(vec,&sum);CHKERRQ(ierr);
      sum  = sum/((PetscScalar)(-1.0*N));
      ierr = VecShift(vec,sum);CHKERRQ(ierr);
    }
  }
  
  if (sp->n) {
    ierr = VecMDot(vec,sp->n,sp->vecs,sp->alpha);CHKERRQ(ierr);
    for (i=0; i<sp->n; i++) sp->alpha[i] = -sp->alpha[i];
    ierr = VecMAXPY(vec,sp->n,sp->alpha,sp->vecs);CHKERRQ(ierr);
  }

  if (sp->remove){
    ierr = (*sp->remove)(sp,vec,sp->rmctx);CHKERRQ(ierr);
  }
  PetscFunctionReturn(0);
}

All I found within seconds using etags. Now you can see exactly where and how the null space is being removed.

   Barry