SVD still enters infinite loops in some cases

[translunar]

I can't seem to reopen #5 myself, so I opened a new one.

Bad news. Still seems to be happening with the following matrix:

{ {7044.7734691220212, 0, 0}, {0, -1.284570679187241e-322, 57.264113734770199}, {0, 0, 0} }

[translunar]

The good news is I think it's an easy fix this time. I tested for S being less than epsilon instead of S equal to 0:

    if (n - k0 == 2 && fabs(S(k0, k0))<eps && fabs(S(k0 + 1, k0 + 1))<eps) { // Compute mu

[wenyan4work]

Sorry to interrupt the discussion about SVD because I am also interested in the topic and did some tests a while ago. I have been a loyal user of PVFMM for a long time. For the matrix types used in the FMM algorithms, pvfmmSVD is consistent and robust.

To write a general purpose numerically robust SVD and guarantee stability over all range of double precision is really hard, and for that purpose I would recommend some pure linear algebra libraries like eigen, lapack or MKL. They usually offer a bunch of different SVD algorithms to fit different matrices. In particular, they all implemented the QR-preconditioned JacobiSVD which is slow but can be really accurate for some types of matrices. BDCsvd in eigen and dgesdd in lapack are also generally stable.

Thanks for the bug notice. I am merging the fixes in my fork of pvfmm.

[translunar]

@wenyan4work Sure. :) I found it because it was included in a StackOverflow answer, and I needed something I could use quickly with raw 2D arrays. I'm not using the whole library — just the algorithm itself. I do want to say how much I appreciate the quick response time from @dmalhotra yesterday.

As a totally random side-note, one of the authors of the chapter this algorithm comes from — Inderjit Dhillon — co-authored a paper with me once at the University of Texas. It was not on this topic. :)

[dmalhotra]

@mohawkjohn I have included your fix in 85534cfbdff62cc26c44c2bba270df6373f63fda. I will leave the issue open for now, in case you encounter more such cases. I think @wenyan4work may be right, this may not be robust enough for all applications.