[petsc-users] MatSetNullSpace

[Bikash Kanungo]

Thank you Matthew for the quick response.

I should provide some more details.

1) This implies that A is symmetric


Yes A is symmetric

2) I think you mean the b is orthogonal to Q


No. I've an extra condition on the solution x - that it has to be
orthogonal to Q. Without this condition the problem will have multiple
solutions. It's same as finding the minimum norm solution


> Its not MatSetNullSpace() that fixes this issue, it is the linear solver
> itself. Many Krylov solvers can converge to the

minimum norm solution of this rank deficient problem.


I was unable to get convergence without using MatSetNullSpace(). So it
seems that the Krlov solvers are inadequate just by themselves.

Second, we remove components in the nullspace from each iterate, which is
> not what you are doing above.
>

This might be the reason why without MatSetNullSpace() I'm unable to attain
convergence.


> It seems easier to just give your Q as the nullspace.

Yeah using Q as nullspace is the easiest option. But my basis is
non-orthogonal and has an overlap matrix S. So in order to provide
orthonormal nullvectors as Q to Petsc, I need to evaluate S^{-1/2}. As of
now, for small problems I can evalaute S^{-1/2}, however, it's a show
stopper for large ones. I can possibly use MatNullSpaceSetFunction() to
circumvent the requirement of orthonormal nullvectors in MatSetNullSpace().
I would appreciate your input on this.

Lastly, I would like to know what operations do MatSetNullSpace enforce
within a KSP solve.

Thanks,
Bikash

On Thu, Jun 1, 2017 at 5:25 PM, Matthew Knepley <knepley at gmail.com> wrote:

> On Thu, Jun 1, 2017 at 2:07 PM, Bikash Kanungo <bikash at umich.edu> wrote:
>
>> Hi,
>>
>> I'm trying to solve a linear system of equations  Ax=b, where A has a
>> null space (say Q) and x is known to be orthogonal to Q. In order to avoid
>> ill-conditioning, I was trying to do the following:
>>
>
> 1) This implies that A is symmetric
>
> 2) I think you mean the b is orthogonal to Q
>
>
>>
>>    1.  Create A as a shell matrix
>>    2.  Overload the MATOP_MULT operation for with my own function which
>>    returns
>>    y = A*(I - QQ^T)x instead of y = Ax
>>    3. Upon convergence, solution = (I-QQ^T)x instead of x.
>>
>> However, I realized that the linear solver can make x have any arbitrary
>> component along Q and still y = A*(I-QQ^T)x will remain unaffected, and
>> hence can cause convergence issues. Indeed, I saw such convergence
>> problems. What fixed the problem was using MatSetNullSpace for A with Q as
>> the nullspace, in addition to the above three steps.
>>
>> So my question is what exactly is MatSetNullSpace doing? And since the
>> full A information is not present and A is only accessed through
>> MAT_OP_MULT, I'm confused as how MatSetNullSpace might be fixing the
>> convergence issue.
>>
>
> Its not MatSetNullSpace() that fixes this issue, it is the linear solver
> itself. Many Krylov solvers can converge to the
> minimum norm solution of this rank deficient problem. Second, we remove
> components in the nullspace from each
> iterate, which is not what you are doing above. It seems easier to just
> give your Q as the nullspace.
>
>   Thanks,
>
>      Matt
>
>
>> Thanks,
>> Bikash
>>
>> --
>> Bikash S. Kanungo
>> PhD Student
>> Computational Materials Physics Group
>> Mechanical Engineering
>> University of Michigan
>>
>>
>
>
> --
> What most experimenters take for granted before they begin their
> experiments is infinitely more interesting than any results to which their
> experiments lead.
> -- Norbert Wiener
>
> http://www.caam.rice.edu/~mk51/
>



-- 
Bikash S. Kanungo
PhD Student
Computational Materials Physics Group
Mechanical Engineering
University of Michigan
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20170601/5a975c94/attachment-0001.html>