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

[Alexander Grayver]

Sorry Barry, I will try to follow this way, but it sometimes scares me 
to get into PETSc source code, so I dare to use your openness. Should 
give it up, though. :)

On 10.01.2012 15:26, Barry Smith wrote:
> 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
>
>
>