> On Oct 12, 2016, at 10:18 PM, Fande Kong <fdkong...@gmail.com> wrote:
>
>
>
> On Wed, Oct 12, 2016 at 9:12 PM, Barry Smith <bsm...@mcs.anl.gov> wrote:
>
> > On Oct 12, 2016, at 10:00 PM, Amneet Bhalla <mail2amn...@gmail.com> wrote:
> >
> >
> >
> > On Monday, October 10, 2016, Barry Smith <bsm...@mcs.anl.gov> wrote:
> >
> > > On Oct 10, 2016, at 4:01 PM, Kong, Fande <fande.k...@inl.gov> wrote:
> > >
> > > Hi All,
> > >
> > > I know how to remove the null spaces from a singular system using
> > > creating a MatNullSpace and attaching it to Mat.
> > >
> > > I was really wondering what is the philosophy behind this? The exact
> > > algorithms we are using in PETSc right now? Where we are dealing with
> > > this, preconditioner, linear solver, or nonlinear solver?
> >
> > It is in the Krylov solver.
> >
> > The idea is very simple. Say you have a singular A with null space N
> > (that all values Ny are in the null space of A. So N is tall and skinny)
> > and you want to solve A x = b where b is in the range of A. This problem
> > has an infinite number of solutions Ny + x* since A (Ny + x*) = ANy +
> > Ax* = Ax* = b where x* is the "minimum norm solution; that is Ax* = b and
> > x* has the smallest norm of all solutions.
> >
> > With left preconditioning B A x = B b GMRES, for example, normally
> > computes the solution in the as alpha_1 Bb + alpha_2 BABb + alpha_3
> > BABABAb + .... but the B operator will likely introduce some component
> > into the direction of the null space so as GMRES continues the "solution"
> > computed will grow larger and larger with a large component in the null
> > space of A. Hence we simply modify GMRES a tiny bit by building the
> > solution from alpha_1 (I-N)Bb + alpha_2 (I-N)BABb + alpha_3 (I-N)BABABAb
> > + .... that is we remove from each new direction anything in the direction
> > of the null space. Hence the null space doesn't directly appear in the
> > preconditioner, just in the KSP method. If you attach a null space to the
> > matrix, the KSP just automatically uses it to do the removal above.
> >
> > Barry, if identity matrix I is of size M x M (which is also the size of A)
> > then are you augmenting N (size M x R; R < M) by zero colums to make I - N
> > possible? If so it means that only first R values of vector Bb are used for
> > scaling zero Eigenvectors of A. Does this choice affect iteration count,
> > meaning one can arbitrarily choose any R values of the vector Bb to scale
> > zero eigenvectors of A?
>
> Yes I wasn't very precise here. I should have written it as as the
> projection of the vector onto the complement of the null space which I think
> can be written as I - N(N'N)^-1N'
> This is done with the routine MatNullSpaceRemove(). The basis of the null
> space you provide to MatNullSpaceCreate() should not effect the convergence
> of the Krylov method at all since the same null space is removed regardless
> of the basis.
>
> I think we need to make sure that the basis vectors are orthogonal to each
> other and they are normalized. Right?
Yes, MatNullSpaceCreate() should report an error if all the vectors are not
orthonormal.
#if defined(PETSC_USE_DEBUG)
if (n) {
PetscScalar *dots;
for (i=0; i<n; i++) {
PetscReal norm;
ierr = VecNorm(vecs[i],NORM_2,&norm);CHKERRQ(ierr);
if (PetscAbsReal(norm - 1.0) > PETSC_SQRT_MACHINE_EPSILON)
SETERRQ2(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D
must have 2-norm of 1.0, it is %g",i,(double)norm);
}
if (has_cnst) {
for (i=0; i<n; i++) {
PetscScalar sum;
ierr = VecSum(vecs[i],&sum);CHKERRQ(ierr);
if (PetscAbsScalar(sum) > PETSC_SQRT_MACHINE_EPSILON)
SETERRQ2(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D
must be orthogonal to constant vector, inner product is
%g",i,(double)PetscAbsScalar(sum));
}
}
ierr = PetscMalloc1(n-1,&dots);CHKERRQ(ierr);
for (i=0; i<n-1; i++) {
PetscInt j;
ierr = VecMDot(vecs[i],n-i-1,vecs+i+1,dots);CHKERRQ(ierr);
for (j=0;j<n-i-1;j++) {
if (PetscAbsScalar(dots[j]) > PETSC_SQRT_MACHINE_EPSILON)
SETERRQ3(PetscObjectComm((PetscObject)vecs[i]),PETSC_ERR_ARG_WRONG,"Vector %D
must be orthogonal to vector %D, inner product is
%g",i,i+j+1,(double)PetscAbsScalar(dots[j]));
}
}
PetscFree(dots);CHKERRQ(ierr);
}
#endif
>
> Fande,
>
>
> Barry
>
>
> >
> > With right preconditioning the solution is built from alpha_1 b +
> > alpha_2 ABb + alpha_3 ABABb + .... and again we apply (I-N) to each term to
> > remove any part that is in the null space of A.
> >
> > Now consider the case A y = b where b is NOT in the range of A. So the
> > problem has no "true" solution, but one can find a least squares solution
> > by rewriting b = b_par + b_perp where b_par is in the range of A and b_perp
> > is orthogonal to the range of A and solve instead A x = b_perp. If you
> > provide a MatSetTransposeNullSpace() then KSP automatically uses it to
> > remove b_perp from the right hand side before starting the KSP iterations.
> >
> > The manual pages for MatNullSpaceAttach() and
> > MatTranposeNullSpaceAttach() discuss this an explain how it relates to the
> > fundamental theorem of linear algebra.
> >
> > Note that for symmetric matrices the two null spaces are the same.
> >
> > Barry
> >
> >
> > A different note: This "trick" is not a "cure all" for a totally
> > inappropriate preconditioner. For example if one uses for a preconditioner
> > a direct (sparse or dense) solver or an ILU(k) one can end up with a very
> > bad solver because the direct solver will likely produce a very small pivot
> > at some point thus the triangular solver applied in the precondition can
> > produce HUGE changes in the solution (that are not physical) and so the
> > preconditioner basically produces garbage. On the other hand sometimes it
> > works out ok.
> >
> >
> > >
> > >
> > > Fande Kong,
> >
> >
> >
> > --
> > --Amneet
> >
> >
> >
> >
>
>
