change implementation of svd() to use Lapack dgesdd instead of dgesvd

[Douglas Bates]

Currently the svd function is implemented using Lapack's dgesvd et al.  The dgesdd subroutine is preferred now as it is considerably faster (http://www.netlib.org/lapack/lug/node71.html)

If I ever get my git repository straightened out I will submit a pull request.

[Douglas Bates]

Actually I have tried both *gesvd and *gesdd and have been unable to show that *gesdd is faster.  They seem to be close to the same speed, which may indicate that later releases of Lapack have used the algorithm for *gesdd in *gesvd

[Jack Poulson]

Hi Douglas,

The routine dgesdd uses a divide and conquer algorithm for the bidiagonal SVD, whereas dgesvd uses a QR algorithm. The former is faster in cases where there is significant deflation and should be generally preferred for large matrices.

Best regards,
Jack Poulson

[Stefan Karpinski]

On Mon, May 14, 2012 at 4:47 PM, Jack Poulson <jack.p...@gmail.com> wrote:

Hi Douglas,

The routine dgesdd uses a divide and conquer algorithm for the bidiagonal SVD, whereas dgesvd uses a QR algorithm. The former is faster in cases where there is significant deflation and should be generally preferred for large matrices.

Best regards,
Jack Poulson

Thank you for the clarification! Perhaps svd should be a polyalgorithm that chooses which routine to use based on matrix size? Would could of course also expose the individual algorithms too.

[Jack Poulson]

The D&C algorithm actually calls the QR algorithm for small subproblems. Its performance depends on the spectrum, so its not easy to know a priori if it will be significantly faster. As soon as I am at a computer I can point everyone to the appropriate literature.

Jack

[Stefan Karpinski]

Great. I really appreciate this feedback (otherwise, we're just shooting in the dark).

[Jack Poulson]

The appropriate reference is LAPACK Working Note 88:
http://www.netlib.org/lapack/lawnspdf/lawn88.pdf

Though it is not explicitly discussed in that paper (AFAIK), the main reason one would choose the implementation based upon the QR algorithm (dgesvd) rather than the one which uses a Divide and Conquer algorithm (dgesdd) is that the latter one requires O(min(m,n)^2) workspace, where k=min(m,n), while the former only requires O(max(m,n)). See the descriptions of LWORK in the following files:
http://www.netlib.org/lapack/double/dgesvd.f
http://www.netlib.org/lapack/double/dgesdd.f

In summary, I recommend always calling the dgesdd if the memory is available, otherwise dgesvd.

Jack

[Viral Shah]

When I originally implemented this, I did what I thought was the safest thing - a lot of people are still using very old versions of LAPACK. We should write the wrappers for dgesdd  in any case. Perhaps, we can pass options to the {svd} function so that it can select the divide and conquer method. The default can be `dgesvd`, but the advanced user can select `dgesdd` if they want it.

-viral

[Jack Poulson]

I was under the impression that dgesdd and friends had been in LAPACK for 12 years (LAPACK 3.0). It seems strange to me that someone who would be using Julia would be dependent upon such an old math library.

On the same subject, I looked in base/linalg_lapack.jl and LAPACK's QR algorithm approach to the Hermitian eigenvalue problem is being used (ssyev/dsyev/cheev/zheev). While there is a much faster Divide and Conquer approach (dsyevd et al) with similar issues as for SVD, the approach based upon Multiple Relatively Robust Representations (MRRR), and available in ssyevr/dsyevr/cheevr/zheevr are often significantly faster, allow for computing subsets of the spectrum, and require less workspace.

I haven't looked into the guts of Julia's build system, but how easy would it be to perform the cmake equivalent of "check_function_exists" on dgesdd/dsyevd/dsyevr to see if they're available?

If anyone is interested, this paper provides a nice overview of the tradeoffs between performance, accuracy, and memory requirements of the symmetric tridiagonal eigensolver approaches:
http://epubs.siam.org/sisc/resource/1/sjoce3/v30/i3/p1508_s1

Jack

[Patrick O'Leary]

On Sunday, May 20, 2012 3:41:55 PM UTC-5, Jack Poulson wrote:

I haven't looked into the guts of Julia's build system, but how easy would it be to perform the cmake equivalent of "check_function_exists" on dgesdd/dsyevd/dsyevr to see if they're available?

Right now it's just raw GNU Make, though Keno Fischer has done some work on CMake for the Windows port, and that may eventually be where the whole build ends up. This would certainly argue for that approach.

[Jeff Bezanson]

I can't imagine needing to work with a very old lapack. lapack is one
of the most open and readily-available, build-anywhere pieces of
software out there.

[Viral Shah]

Ok, my bad. I somehow thought that dgesdd was a recent addition. In that case, we should just go ahead and add it. For now, we could call it using an Option, as we figure out the right polyalgorithm. We should do the same for QR as well. Let's add the wrappers and get these to work to start with.

-viral

[Stefan Karpinski]

I think that Jack mentioned that LAPACK already uses basically a polyalgorithm. Not positive on that.

[Jack Poulson]

Yes, the LAPACK routine dgesdd uses a Divide and Conquer approach at its core, which recurses on smaller and smaller problems; at the smallest level, it uses a QR algorithm which is identical to that of dgesvd.

The only real reason to use dgesvd instead is for accuracy and workspace reasons. The former should only be important if you care about high relative accuracy for tiny singular values.

Jack

[Stefan Karpinski]

Ok, seems clear that we should use dgesdd and just expose dgesvd as a separate function for particular uses.

[John Cowan]

On Sun, May 20, 2012 at 4:58 PM, Patrick O'Leary
<patrick...@gmail.com> wrote:

> Right now it's just raw GNU Make, though Keno Fischer has done some work on
> CMake for the Windows port, and that may eventually be where the whole build
> ends up. This would certainly argue for that approach.

Do not, if you value your life and your soul, go there. Autotools and
CMake alike are inventions of the Devil, but CMake is worse.

As I may have said before, the only environments that matter any more
are Linux, BSD, Mac OS X, Solaris, Haiku, Cygwin, MinGW, MinGW+Msys.
Create an appopriate .h file containing the correct settings of HAVE_*
for whatever names you actually require, and leave it at that.

--
GMail doesn't have rotating .sigs, but you can see mine at
http://www.ccil.org/~cowan/signatures

[Stefan Karpinski]

On Mon, May 21, 2012 at 2:16 AM, John Cowan <johnw...@gmail.com> wrote:

On Sun, May 20, 2012 at 4:58 PM, Patrick O'Leary
<patrick...@gmail.com> wrote:

> Right now it's just raw GNU Make, though Keno Fischer has done some work on
> CMake for the Windows port, and that may eventually be where the whole build
> ends up. This would certainly argue for that approach.

Do not, if you value your life and your soul, go there.  Autotools and
CMake alike are inventions of the Devil, but CMake is worse.

This may be my favorite email ever posted on this list. 

[Patrick O'Leary]

On Monday, May 21, 2012 1:16:50 AM UTC-5, John Cowan wrote:

On Sun, May 20, 2012 at 4:58 PM, Patrick O'Leary
<patrick...@gmail.com> wrote:

> Right now it's just raw GNU Make, though Keno Fischer has done some work on
> CMake for the Windows port, and that may eventually be where the whole build
> ends up. This would certainly argue for that approach.

Do not, if you value your life and your soul, go there.  Autotools and
CMake alike are inventions of the Devil, but CMake is worse.

And plain Makefiles are also inherently broken, particularly with the build-the-world approach taken by julia. All build systems suck; find the one that sucks least for your application and run with it. (Disclaimer: I don't know what the least-bad build system for julia is.)

As I may have said before, the only environments that matter any more
are Linux, BSD, Mac OS X, Solaris, Haiku, Cygwin, MinGW, MinGW+Msys.
Create an appopriate .h file containing the correct settings of HAVE_*
for whatever names you actually require, and leave it at that.

Haiku is more important than MSVC? I'll have to disagree with you there. MSVC generates good code, for better or worse. There's a reason Mozilla puts in the effort to use MSVC for their Windows builds of Firefox, etc.

Of course what would be great would be to have julia building entirely on LLVM. Now that would be nice. But not much to do with the build system.

[John Cowan]

On Mon, May 21, 2012 at 2:35 AM, Patrick O'Leary
<patrick...@gmail.com> wrote:

> Haiku is more important than MSVC? I'll have to disagree with you there.
> MSVC generates good code, for better or worse.

Sure. I just didn't think MSVC was practical, given that you depend
on Fortran so much. If it is, all the better!

Haiku is just a target of opportunity, and I'd have no problem with
not supporting it.

> Of course what would be great would be to have julia building entirely on
> LLVM. Now that would be nice. But not much to do with the build system.

Indeed.

[Neuwirth Erich]

I am trying to build on my Mac, but the build fails.
I was able to build before.


Mahler:julia neuwirth$ make
CC src/jltypes.o
In file included from support/ios.h:6,
from support/libsupport.h:9,
from julia.h:8,
from jltypes.c:13:
support/../../usr/include/uv.h:58:18: error: ares.h: No such file or directory
In file included from support/ios.h:6,
from support/libsupport.h:9,
from julia.h:8,
from jltypes.c:13:
support/../../usr/include/uv.h:1052: error: expected declaration specifiers or ‘...’ before ‘ares_channel’
support/../../usr/include/uv.h:1052: warning: ‘struct ares_options’ declared inside parameter list
support/../../usr/include/uv.h:1052: warning: its scope is only this definition or declaration, which is probably not what you want
support/../../usr/include/uv.h:1055: error: expected declaration specifiers or ‘...’ before ‘ares_channel’
support/../../usr/include/uv.h:1528: error: expected specifier-qualifier-list before ‘ares_channel’
make[2]: *** [jltypes.o] Error 1
make[1]: *** [julia-release] Error 2
make: *** [release] Error 2

[Jameson Nash]

try `make -C deps clean-uv`, then try `make` again

[Stefan Karpinski]

I think that libuv needs a forced rebuild. Try "make -C deps clean-uv && make testall"

[Stefan Karpinski]

jinx.

[Jack Poulson]

On Mon, May 21, 2012 at 1:35 AM, Patrick O'Leary <patrick...@gmail.com> wrote:

On Monday, May 21, 2012 1:16:50 AM UTC-5, John Cowan wrote:

On Sun, May 20, 2012 at 4:58 PM, Patrick O'Leary
<patrick...@gmail.com> wrote:

> Right now it's just raw GNU Make, though Keno Fischer has done some work on
> CMake for the Windows port, and that may eventually be where the whole build
> ends up. This would certainly argue for that approach.

Do not, if you value your life and your soul, go there.  Autotools and
CMake alike are inventions of the Devil, but CMake is worse.

And plain Makefiles are also inherently broken, particularly with the build-the-world approach taken by julia. All build systems suck; find the one that sucks least for your application and run with it. (Disclaimer: I don't know what the least-bad build system for julia is.)

Regardless of everyone's personal least-hated build system, I have found the functionality provided by check_function_exists invaluable for various situations, including various tests on user-specified BLAS/LAPACK libraries:

1) Is there an underscore at the end of every BLAS symbol? (varies on the same platform)

2) Is there an underscore at the end of every LAPACK symbol? (varies on the same platform, INDEPENDENTLY of the BLAS)

3) Was BLAS/LAPACK provided via MKL? If so, thread safety is different in certain situations.

Each of these issues is easily handled with check_function_exists. I am not sure how a Makefile would handle this without specifying the particular tools used on each platform.

Jack

[Viral Shah]

Thankfully, the issue of the double underscores has not yet come up. I have been silently waiting in anticipation. Also, people have linked MKL without much trouble, which is nice. I would prefer continuing with our creaky build system until it becomes completely unbearable. We'll have cmake anyways for Windows, and that will allow us to experiment and do a smooth transition if we do want to migrate to cmake.

-viral

[Erich Neuwirth]

Now (after quite a fee warnings related to ares) I get

/usr/bin/ranlib: file: uv.a(inet_ntop.o) has no symbols

    FLISP /src/julia_flisp.boot

arg-error: io.tostring!: requires memory stream

in file ./mk_julia_flisp_boot.scm

#0 fatal error:

(arg-error "io.tostring!: requires memory stream")

make[2]: *** [julia_flisp.boot] Error 1

make[1]: *** [julia-release] Error 2

[Jameson Nash]

You also need to do `make clean` in the top level (or `make -C src clean`, if that doesn't work)

[John Myles White]

If none of that works, I can say that just deleting the libuv directory worked for me on the first try.

 -- John

[Federico Calboli]

I tried all of the above but nothing works (on Ubuntu 12.10), but I still get:

    FLISP src/julia_flisp.boot

arg-error: io.tostring!: requires memory stream

in file ./mk_julia_flisp_boot.scm

#0 fatal error:

(arg-error "io.tostring!: requires memory stream")

make[2]: *** [julia_flisp.boot] Error 1

make[1]: *** [julia-release] Error 2

make: *** [release] Error 2

BW

F

[Jameson Nash]

When FLISP is complaining it can't compile something, I find `make cleanall` is the best medicine :)

--