Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add xggsvd3 wrappers. #14444

Merged
merged 1 commit into from
Dec 28, 2015
Merged

Add xggsvd3 wrappers. #14444

merged 1 commit into from
Dec 28, 2015

Conversation

andreasnoack
Copy link
Member

The subroutines replace xggsvd in LAPACK 3.6.0 so calls are made conditionally on the existence of the symbols.

Update: I've also changed the blasfunc function to a macro because the function was awkward to use without a loop because of the @eval block. E.g. the new LAPACK.laver() definition.

cc: @yuyichao

@@ -77,7 +77,12 @@ end
GeneralizedSVD{T}(U::AbstractMatrix{T}, V::AbstractMatrix{T}, Q::AbstractMatrix{T}, a::Vector, b::Vector, k::Int, l::Int, R::AbstractMatrix{T}) = GeneralizedSVD{T,typeof(U)}(U, V, Q, a, b, k, l, R)

function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T})
U, V, Q, a, b, k, l, R = LAPACK.ggsvd!('U', 'V', 'Q', A, B)
# xggsvd3 replaced xggsvd in LAPACK 3.6.0
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can we initialize this in __init__ or use cglobal to avoid a dlsym on every function call?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That would be cleaner even though this function is so expensive that I don't think it matters much. I'll update the PR.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have now defined a constant with the LAPACK version number instead.


import ..LinAlg: BlasFloat, Char, BlasInt, LAPACKException,
DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare

function __init__()
const global VERSION = VersionNumber(laver()...)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is actually bad .... For one thing, pre-compilation won't treat the version number as a constant.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can you explain a bit more? When is the constant set? It should correspond to the version actually linked when Julia is launched. What is the right solution?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ref #12010 and #13817

and `V'*B*Q = D2*R`. `D1` has `alpha` on its diagonal and `D2` has `beta` on its
diagonal. If `jobu = U`, the orthogonal/unitary matrix `U` is computed. If
`jobv = V` the orthogonal/unitary matrix `V` is computed. If `jobq = Q`,
the orthogonal/unitary matrix `Q` is computed. If `job{u,v,q} = N`, that
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

curly braces are a little confusing here since they mean type parameters in julia code

@tkelman
Copy link
Contributor

tkelman commented Dec 19, 2015

What error do you get in the case that your lapack is too new or too old and you try to call the LAPACK. wrapper function directly? If these are interchangeable and we won't have both available at the same time, should we put them in the same Julia wrapper method?

(Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
major, minor, patch)
return major[1], minor[1], patch[1]
end
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Although this function is only used during initialization time, maybe it's better to just return a VersionNumber from this function? Or is a tuple easier in some cases? Also I think allocating a Ref{BlasInt}() is faster than BlasInt[0].

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The latest version uses Refs, but I've kept the tuple return because, in general, the return from the LAPACK is similar to what is returned from the subroutines. Sometimes we make some postprocessing in the wrappers, but most often this is done down the pipeline and I don't think we return "special" Julia types from any of the LAPACK wrappers.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Fair enough.

I have no idea about the function you are wrapping but it otherwise LGTM. =)

@andreasnoack
Copy link
Member Author

@tkelman The old unblocked and the new blocked API are very similar. However, they have added a build flag for deprecated functions which would give both versions and I'm not sure how MKL and Accelerate will handle a deprecation like this. It is also slightly clearner to have two functions instead of one with more branching but, of course, with the cost of some code duplication.

@tkelman
Copy link
Contributor

tkelman commented Dec 23, 2015

If they're completely redundant then I disagree that two functions is cleaner than one, but it sounds like under some configurations they might not be.

Don't use job{u,v,q} in the docstring, that's too easily confusable with Julia type parameter syntax which isn't what you mean there.

…so calls are made conditionally

 on the existence of the symbols by defining a constant with the LAPACK version number.

Make blasfunc a macro such that it can be called without @eval $name.
@andreasnoack
Copy link
Member Author

Don't use job{u,v,q} in the docstring, that's too easily confusable with Julia type parameter syntax which isn't what you mean there.

I forgot your comment to the documentation. The entry was copied verbatim from ggsvd!, but I have now changed both versions.

andreasnoack added a commit that referenced this pull request Dec 28, 2015
@andreasnoack andreasnoack merged commit 57eb9c1 into master Dec 28, 2015
@andreasnoack andreasnoack deleted the anj/ggsvd3 branch December 28, 2015 13:39
major = BlasInt[0]
minor = BlasInt[0]
patch = BlasInt[0]
ccall((@blasfunc(ilaver_), liblapack), Void,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be :ilaver_ ? It's inside a macro call so apparently works anyway, but looks a little strange.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think a macro can handle :ilaver_ and ilaver_ equally well (probably even easier without the quote) so may as well save some typing?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just looks strange to be using what looks like, but isn't, a defined Julia object.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

julia> macro whatisthis(thing)
           return typeof(thing)
       end

julia> @whatisthis(:ilaver_)
Expr

julia> @whatisthis(ilaver_)
Symbol

julia> for s in (:ilaver_,)
           @eval println(@whatisthis($s))
       end
Symbol

I was also puzzled by this. The macro could be modified to work for Expr as well.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If we can get away without the $ on the @eval version that would seem cleaner to me

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

3 participants