Changes for an Lapack module

[dmbates]

Add base/Lapack.jl containing a module with the definitions of the direct Lapack interface functions. Code from base/{factorizations,linalg_lapack,linalg_specialized}.jl was moved to this file. The order of files in sysimg.jl was changed to accommodate these changes. New functions were added to the list of exports in export.jl.

Currently this code fails the test/lapack.jl script because the first call to eig() fails. For example

julia> a    = rand(3,3)
3x3 Float64 Array:
 0.488698  0.106265  0.0258011
 0.510604  0.646657  0.648978 
 0.878416  0.1577    0.742774 

julia> asym = a' + a
3x3 Float64 Array:
 0.977395  0.616869  0.904217
 0.616869  1.29331   0.806678
 0.904217  0.806678  1.48555 

julia> ishermitian(asym)
true

julia> eig(asym)
type error: apply: expected Function, got IntrinsicFunction

julia> eig(asym)
([0.292192, 0.625846, 2.83822],
3x3 Float64 Array:
  0.800213     0.322779  0.505443
 -0.00990314  -0.835577  0.549284
 -0.599635     0.444549  0.665443)

[JeffBezanson]

The error is my bug, and it will be fixed on master soon. Then there is a problem that in geev! the statement Rtyp = typeof(real(work[1])) occurs before work has been defined.

[dmbates]

Thanks for fixing the problem and spotting the other problem. I will fix the second and (fingers crossed) rebase after your commit. Should I then create a new pull request and close this one?

[pao]

You can force-push (i.e., git push -f mine HEAD) to your LapackMod branch, which will update this pull request with that code.

[dmbates]

I think this should be good to go now. make test-lapack succeeds on my system. Sorry for not squashing the commits but I don't want to get too tricky with git just before I am supposed to leave for the day. That seems to be a recipe for disaster.

[JeffBezanson]

This is a great reorganization and I look forward to merging very soon!

[ViralBShah]

This is awesome. I am not a fan of the names lud and such, but that should not prevent this pull request from being merged. As we add all this to the documentation and use it as we go ahead, we will figure out better names, or get used to these names.

I suggest a little more work before merging:

1. Move all the BLAS interfacing also into Lapack.jl (The whole interface is currently only 1200 lines!)
2. Merge all the linalg_* files into one linalg.jl (will add up to about a 1000 lines, which is manageable)
3. Stick all of linalg_blas.jl into linalg.jl as well, so that we only have Lapack.jl and linalg.jl. There are some abstract definitions in linalg.jl now, and there may be some merit in keeping them separate.

The code has actually become a lot more compact since we started the work on lapack, and we have a whole lot more functionality in there now.

[ViralBShah]

BTW, @StefanKarpinski has said that we should follow the convention of using only lower case characters for filenames as convention. I have not yet reverted DSP.jl and DSP_fftw.jl, because once you introduce upper case names, it is a real pain on platforms where the filesystem ignores the case (OS X, for example, and perhaps Windows also).

It would be best not to introduce Lapack.jl, and instead introduce lapack.jl in the pull request. I am not sure what is the best way to do this. For my RNG stuff, I just found it easier to junk my pull request and create a new one.

[JeffBezanson]

I would prefer to separate the blas and lapack interfaces. They are separate libraries after all.

[ViralBShah]

They really are joined at the hip, and we do ship them as one library (libopenblas). This only combines the interfaces in the Lapack module, while still allowing the flexibility to ship different libraries for BLAS and LAPACK.

I recollect discussion in some other issue where it seemed like merging these was the right thing to do. I do not have a strong opinion on this one way or another.

[JeffBezanson]

Yes, in our default setup they're linked together, but lapack calls blas and not the other way around, so there is an abstraction layer there.

[dmbates]

I agree that Lapack subroutines call BLAS subroutines but I don't think that is relevant here. In Julia you just see a ccall to an Lapack subroutine and the result of that call, Whether or not that call involves calls to BLAS subroutines is irrelevant at the Julia level. The only thing that is important is that the Lapack DLL also links to the BLAS DLL and that is taken care of during the dlopen calls in start_img.jl.

I have committed changes to the LapackMod branch on my repository to replace $string(foo) with $(string(foo)), to rename Lapack.jl to lapack.jl and also to fix a few calls in testt/factorizations.jl.

I agree with moving all the linalg_* files into linalg.jl. If there is agreement on moving the BLAS interface calls into lapack.jl I will do that too and prepare a new pull request.

My LapackMod branch is becoming a mess. If there are no objections I will close this pull request, start from a clean version of master and create a new pull request. My relationship with git is improving but still has a long way to go before we can be considered friends.

[timholy]

This is just great. I really like how you centralized some of the checks, and creating the SymTridiagonal type is a great idea. I also liked the one-line comments in some places on what the different Lapack routines do; once this gets merged I may contribute more of those over time.

@ViralBShah , on the lud issue: I'd urge you to suggest alternate names now, if you can come up with something you like. There's nothing good about changing the API too often.

Do we want to create a SVDDense type? Again, I'd say something to consider after merging...this is too good to let languish!

[ViralBShah]

Let us keep BLAS separate for now, but in a Blas module. Feel free to create a new pull request.

-viral

On 16-Sep-2012, at 7:51 PM, dmbates [email protected] wrote:

I agree that Lapack subroutines call BLAS subroutines but I don't think that is relevant here. In Julia you just see a ccall to an Lapack subroutine and the result of that call, Whether or not that call involves calls to BLAS subroutines is irrelevant at the Julia level. The only thing that is important is that the Lapack DLL also links to the BLAS DLL and that is taken care of during the dlopen calls in start_img.jl.

I have committed changes to the LapackMod branch on my repository to replace $string(foo) with $(string(foo)), to rename Lapack.jl to lapack.jl and also to fix a few calls in testt/factorizations.jl.

I agree with moving all the linalg_* files into linalg.jl. If there is agreement on moving the BLAS interface calls into lapack.jl I will do that too and prepare a new pull request.

My LapackMod branch is becoming a mess. If there are no objections I will close this pull request, start from a clean version of master and create a new pull request. My relationship with git is improving but still has a long way to go before we can be considered friends.

—
Reply to this email directly or view it on GitHub.

[ViralBShah]

I will think of alternate names.

-viral

On 16-Sep-2012, at 8:22 PM, Tim Holy [email protected] wrote:

This is just great. I really like how you centralized some of the checks, and creating the SymTridiagonal type is a great idea. I also liked the one-line comments in some places on what the different Lapack routines do; once this gets merged I may contribute more of those over time.

@ViralBShah , on the lud issue: I'd urge you to suggest alternate names now, if you can come up with something you like. There's nothing good about changing the API too often.

Do we want to create a SVDDense type? Again, I'd say something to consider after merging...this is too good to let languish!

—
Reply to this email directly or view it on GitHub.

[dmbates]

Looking at base/linalg_blas.jl I'm not sure how it can be loaded after the changes in the semantics of the $ operator. That file still has instances of $string(fname)

Right now I can't test changes on a fresh pull of the master branch because of the segfault issue I mentioned on julia-dev.

[andreasnoack]

Some linalg commments!

I have translated some code from Matlab and have some suggestions for the linear algebra module in Julia. My understanding is that a major revision of the module is underway by Douglas Bates and hence I guess it is easier to write my comments here. Please correct me if I am wrong as I am new to Julia and contributing in general.

I found a bug in _jl_lapack_gels and I have made a pull request with fix which pao was kind to redirect to this place.

When mulitplying zero-dimension matrix as fx randn(4, 0) * randn(0, 4) a 4x4 zero matrix is returned which i very convenient. However, _jl_gemm at line 430 in linalg_blas.jl passes stride(A, 2) and stride(B, 2) which may be zero. The BLAS subroutine should have a value of at least one so I guess the easiest fix would be just to write max(stride(B, 2), 1). It works on my machine. The error is only shown when Julia is build with MKL, but the requirement is general BLAS.

I often use thin SVDs and hence I have made the method

function svd{T<:Union(Float64,Float32,Complex128,Complex64)}(A::StridedMatrix{T}, vecs::Int64)
if vecs != 0
error("Not supported")
else
Base._jl_lapack_gesvd('S', 'A', copy(A))
end
end

to avoid the computation of the large U matrix. I believe it would be a good idea with such a method in Julia. A side note. I read some comments about the two SVD implementations. I vote for just one SVD function at the user level. I have no opinions on which of the Lapack routines it should call.

There is a svdvals but no eigvals. I think it would be nice with a way to avoid the computation of the vectors.

[dmbates]

@ViralBShah Have you considered alternative names for lud, chold, etc?

[ViralBShah]

I actually prefer making lu return the factorization object, even though it breaks matlab compatibility.

Once we have lud for dense, would we use lus for sparse?

-viral

On 18-Sep-2012, at 9:52 PM, dmbates [email protected] wrote:

@ViralBShah Have you considered alternative names for lud, chold, etc?

—
Reply to this email directly or view it on GitHub.

[dmbates]

Actually the proposed convention used the 'd' for 'decomposition', not 'dense'. So

lu - Matlab compatible, lu(A::StridedMatrix) = factors(lud(A))
lud - LU decomposition of (a copy of) A
lud! - LU decomposition overwriting A

I don't have strong opinions on whether or not there should be both an lu and an lud function. When an appropriate factors method is defined it seems redundant to have both.

By the way, I haven't been able to write a factors method for the BunchKaufman type because I can't work out what the storage is. The descriptions in the Lapack code make it sound like it is a generalization of the LDLt where D is symmetric block-diagonal with blocks of size 1 by 1 or 2 by 2. When the symmetric matrix is positive definite the blocks are all 1 by 1 but I can't make sense of the diagonal.

[timholy]

@ViralBShah, aside from the general problems from incompatibility, we've realized that the same logic applied to chol would introduce bugs, not errors. That seems like a particularly bad way to break compatibility, hence the new names. See this discussion:
826c42f#commitcomment-1823621

Would ludcmp be better? I came to this by realizing that if we do this for svd, then svdd seems dumb, but svdcmp would be clear as the analog of the Matlab-compatible svd. But maybe the SVD doesn't have a representation for specialized types that's any better than dense, anyway, in which case there is no need to do this for SVD.

[timholy]

@dmbates, presumably the first person who needs factors of the BunchKaufman decomposition will figure this out :-). Thanks for adding it, I'd never heard about it before.

[ViralBShah]

Ok, I recollect the discussion now. How about lufact? Should we standardize on using factorizations as the terminology in names, and avoid the name decomposition? Basically, a consistent naming scheme one way or another.

-viral

On 18-Sep-2012, at 11:49 PM, dmbates [email protected] wrote:

Actually the proposed convention used the 'd' for 'decomposition', not 'dense'. So

lu - Matlab compatible, lu(A::StridedMatrix) = factors(lud(A))
lud - LU decomposition of (a copy of) A
lud! - LU decomposition overwriting A

I don't have strong opinions on whether or not there should be both an lu and an lud function. When an appropriate factors method is defined it seems redundant to have both.

—
Reply to this email directly or view it on GitHub.

[dmbates]

I feel kind of guilty about the term factorizations rather than decompositions because I was the one who introduced the term. In many ways I think that ludecomp, choldecomp, ... sound better than lufact, cholfact, ... but I am happy to go with whatever decision is made.

I am going to close this pull request now and open up another on a clean branch.

[StefanKarpinski]

Would it be nuts to use the uppercase version to indicate that you want a decomposition object? So LU(X) versus lu(X). The idea (which may not currently work) would be that LU could be an abstract supertype of the LUDense and LUSparse, etc. decomposition types, but also a constructor for such decompositions. The problem is that we would need support for defining methods on abstract types. cc: @JeffBezanson.

Alternative idea would be to have a decomposition package that includes the decomposition versions of all the Matlab style calls. Don't really care for that, however. I'd be more inclined to just use the Matlab names and have people call them differently. The one that's problematic — chol, right? — can just be called chold to avoid the unobvious breakage.

[ViralBShah]

We have reserved the uppercase for things like type and module names, and
do not really do it for basic library functions. I would hate to press
Shift that many times, if we start allowing uppercase for function names
too - but I like the current consistency.

I'd be ok with either name, and we can still switch from factorizations to
decompositions. I only like factorizations better because fact is shorter
than suffixing decomp. @dmbates - your call.

If anything, the decomposition versions should be the default, since they
are much more efficient. In order to be matlab compatible, you actually end
up doing much more work, and have an interface that encourages a style of
programming that leads to poor performance.

I would think that the simplest and cleanest thing to do right now is to
use the decomp/fact suffixes, and retain matlab behaviour for
matlab-compatible names. It would be nice not to muck with matlab behaviour
here, since this is what matlab is really popular for.

-viral

On Wed, Sep 19, 2012 at 10:46 PM, Stefan Karpinski <[email protected]

wrote:

Would it be nuts to use the uppercase version to indicate that you want a
decomposition object? So LU(X) versus lu(X). The idea (which may not
currently work) would be that LU could be an abstract supertype of the
LUDense and LUSparse, etc. decomposition types, but also a constructor
for such decompositions. The problem is that we would need support for
defining methods on abstract types. cc: @JeffBezansonhttps://github.com/JeffBezanson
.

Alternative idea would be to have a decomposition package that includes
the decomposition versions of all the Matlab style calls. Don't really care
for that, however. I'd be more inclined to just use the Matlab names and
have people call them differently. The one that's problematic — chol,
right? — can just be called chold to avoid the unobvious breakage.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/1281#issuecomment-8698677.

-viral

[ViralBShah]

@timholy Even in SVD, IIRC, you can do the decomposition, but choose not to
construct U and V, or construct them later, etc. It would be nice to expose
some of the more general features of the lapack design in the decomposition
interface.

-viral

On Wed, Sep 19, 2012 at 12:08 AM, Tim Holy [email protected] wrote:

@ViralBShah https://github.com/ViralBShah, aside from the general
problems from incompatibility, we've realized that the same logic applied
to chol would introduce bugs, not errors. That seems like a particularly
bad way to break compatibility, hence the new names. See this discussion:
826c42f#commitcomment-1823621http://JuliaLang/julia/commit/826c42f4eca1f8a939df5484036fcc4412c55b59#commitcomment-1823621

Would ludcmp be better? I came to this by realizing that if we do this
for svd, then svdd seems dumb, but svdcmp would be clear as the analog of
the Matlab-compatible svd. But maybe the SVD doesn't have a
representation for specialized types that's any better than dense, anyway,
in which case there is no need to do this for SVD.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/1281#issuecomment-8665307.

-viral

[dmbates]

@ViralBShah I don't think it is possible to create an SVD in a compact form and later use that form to generate U and/or V' Neither *gesvd nor *gesdd allow that. You can do the reduction to bidiagonal form like that in *gebrd, which I now see I have wrong in #1290. Grrr. It is like that tiny window when you publish a book between getting the printed copy in your hands and receiving the first errata report.

[ViralBShah]

@dmbates, thanks for pointing out. It is late, and I got lazy to look up
the the lapack definitions.

On Wed, Sep 19, 2012 at 11:32 PM, dmbates [email protected] wrote:

@ViralBShah https://github.com/ViralBShah I don't think it is possible
to create an SVD in a compact form and later use that form to generate U
and/or V' Neither *gesvd nor *gesdd allow that. You can do the reduction
to bidiagonal form like that in *gebrd, which I now see I have wrong in
#1290 #1290. Grrr. It is like
that tiny window when you publish a book between getting the printed copy
in your hands and receiving the first errata report.

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/1281#issuecomment-8700363.

-viral

[StefanKarpinski]

For SVD, I would argue (tentatively) that providing left and right SVD variants is probably the best way to go. Maybe

U,S = lsvd(X)
S,V = rsvd(X)

or

U,S = svdl(X)
S,V = svdr(X)

[dmbates]

Why not use multiple dispatch and create an svd method taking two additional Char arguments that get passed to Lapack.gesvd

[StefanKarpinski]

Passing Char arguments to change behavior is not a good interface, despite what Fortran has led people to believe.

[dmbates]

How about Bool arguments? They don't quite work in this case because there are different versions of U, for example. However, much of the usage of Char arguments in Lapack is for binary switches ('U'pper/'L'ower triangles, 'U'nit/'N'onunit diagonal, 'L'eft/'R'ight multipliers, 'T'ranspose/'N'oTranspose, etc.) To me it would be cleaner to make these switches Bool variables right up to the time they get passed to the Lapack routines. Then you don't need all the checking if only the two possible states are represented and that upper case is used, etc.

[StefanKarpinski]

I'm not quite clear on what we're talking about anymore. Just whether to compute the left and/or right factors? I guess there are three flags you could potentially pass in: boolean for each of the three components U, S and V. Giving false for any of them would result in not returning that one. This is problematic for type inference though, where writing something in the spirit of U,S = svd(X,right=false), pretending we have keyword args for the moment, would require analysis of values to determine the types of the output arguments (assuming that S is returned as a compact DiagMatrix type or something like that). If you use lsvd and rsvd then each function always returns the same types and the type inference system can figure that out.

[dmbates]

Okay, I can see why making type inference easier is a valid reason.

In my comment on different versions of U I meant that for an m by n matrix A with m > n you can generate either an m by n matrix U or an m by m matrix U. Most of the time you want the smaller one (m by n).

[JeffBezanson]

Just asked Alan about terminology and he says "singular value decomposition" but "LU factorization". History being what it is.

My big question is: would the performance vs. compatibility problem be totally solved if the syntax U,S,V = svd_object(A) could be used to request full matrices from a decomposition object? Or are there other important cases that would also break? If this would totally fix things, it would be worth figuring out how to get it.

[nolta]

For what it's worth, mathematica decided to call everything decompositions: LUDecomposition, QRDecomposition, etc.

[dmbates]

wikipedia.org uses "decomposition" for Cholesky, LU, QR and SVD and they don't have anything to say about Bunch-Kaufman. Not that I would hold wikipedia.org to be the definitive authority on usage.

As far as I can recall, the SVD is the only decomposition that has so many options. If you look at the comments in deps/openblas-v0.2.3/lapack-3.4.1/SRC/dgesvd.f the left or right singular vectors can each be option 'A' (all of them - the result is square), 'S' (small version - number of columns is min(m,n)), 'N' for none and 'O' for overwrite which is probably not relevant here. To get all these options would require

svdvals - 'N', 'N' - singular values only
svd - 'S', 'S' - both left and right
svdl - 'S', 'N' - left only
svdr - 'N', 'S' - right only

plus some version of the last two for combinations like 'A', 'N', I would suggest svdL for 'A','N' but that brings upper-case letters into the name.

To me it seems that we are getting too detailed. If someone really wants one of these arcane combinations they should call Lapack.gesvd! or Lapack.gesdd! directly.

[JeffBezanson]

I agree. A big reason to have easily-accessible direct bindings to lapack is to get those extra behaviors in case you want them.