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 > > >
