|
2 | 2 | using LinearAlgebra: BlasFloat, Char, BlasInt, LAPACK, LAPACKException,
|
3 | 3 | DimensionMismatch, SingularException, PosDefException, chkstride1, checksquare,
|
4 | 4 | triu!
|
5 |
| -using LinearAlgebra.BLAS: @blasfunc, libblas, BlasReal, BlasComplex |
6 |
| -using LinearAlgebra.LAPACK: liblapack, chklapackerror |
7 | 5 |
|
8 | 6 | set_num_blas_threads(n::Integer) = LinearAlgebra.BLAS.set_num_threads(n)
|
9 |
| - |
10 |
| -# TODO: Remove once we drop support for Julia 1.5 or less |
11 |
| -function get_num_blas_threads() |
12 |
| - blas = LinearAlgebra.BLAS.vendor() |
13 |
| - # Wrap in a try to catch unsupported blas versions |
14 |
| - if blas == :openblas |
15 |
| - return ccall((:openblas_get_num_threads, Base.libblas_name), Cint, ()) |
16 |
| - elseif blas == :openblas64 |
17 |
| - return ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ()) |
18 |
| - elseif blas == :mkl |
19 |
| - return ccall((:mkl_get_max_threads, Base.libblas_name), Cint, ()) |
20 |
| - end |
21 |
| - |
22 |
| - # OSX BLAS looks at an environment variable |
23 |
| - @static if Sys.isapple() |
24 |
| - return tryparse(Cint, get(ENV, "VECLIB_MAXIMUM_THREADS", "1")) |
25 |
| - end |
26 |
| - error("unknown BLAS") |
27 |
| -end |
| 7 | +get_num_blas_threads(n::Integer) = LinearAlgebra.BLAS.get_num_threads(n) |
28 | 8 |
|
29 | 9 | # TODO: define for CuMatrix if we support this
|
30 | 10 | function _one!(A::DenseMatrix)
|
@@ -108,16 +88,6 @@ function _leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{QR, QRpos}, atol::
|
108 | 88 | end
|
109 | 89 | return Q, R
|
110 | 90 | end
|
111 |
| -# TODO: reconsider the following implementation |
112 |
| -# Unfortunately, geqrfp seems a bit slower than geqrt in the intermediate region |
113 |
| -# around matrix size 100, which is the interesting region. => Investigate and maybe fix |
114 |
| -# function _leftorth!(A::StridedMatrix{<:BlasFloat}) |
115 |
| -# m, n = size(A) |
116 |
| -# A, τ = geqrfp!(A) |
117 |
| -# Q = LAPACK.ormqr!('L', 'N', A, τ, eye(eltype(A), m, min(m, n))) |
118 |
| -# R = triu!(A[1:min(m, n), :]) |
119 |
| -# return Q, R |
120 |
| -# end |
121 | 91 |
|
122 | 92 | function _leftorth!(A::StridedMatrix{<:BlasFloat}, alg::Union{QL, QLpos}, atol::Real)
|
123 | 93 | iszero(atol) || throw(ArgumentError("nonzero atol not supported by $alg"))
|
@@ -282,6 +252,9 @@ function _svd!(A::StridedMatrix{<:BlasFloat}, alg::Union{SVD, SDD})
|
282 | 252 | end
|
283 | 253 |
|
284 | 254 | # TODO: override Julia's eig interface
|
| 255 | +# using LinearAlgebra.BLAS: @blasfunc, libblas, BlasReal, BlasComplex |
| 256 | +# using LinearAlgebra.LAPACK: liblapack, chklapackerror |
| 257 | +# |
285 | 258 | # eig!(A::StridedMatrix{<:BlasFloat}) = LinAlg.LAPACK.gees!('V', A)
|
286 | 259 | # function eig!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) where T<:BlasReal
|
287 | 260 | # n = size(A, 2)
|
@@ -312,43 +285,49 @@ end
|
312 | 285 | # return Eigen(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]...)
|
313 | 286 | # end
|
314 | 287 | #
|
315 |
| -# |
316 | 288 | # eigfact!(A::RealHermSymComplexHerm{<:BlasReal, <:StridedMatrix}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
|
317 | 289 | #
|
318 |
| -# |
319 |
| -# |
320 |
| -# Modified / missing Lapack wrappers |
321 |
| -#------------------------------------ |
322 |
| -# geqrfp!: computes qrpos factorization, missing in Base |
323 |
| -geqrfp!(A::StridedMatrix{<:BlasFloat}) = |
324 |
| - ((m, n) = size(A); geqrfp!(A, similar(A, min(m, n)))) |
| 290 | +# TODO: reconsider the following implementation |
| 291 | +# Unfortunately, geqrfp seems a bit slower than geqrt in the intermediate region |
| 292 | +# around matrix size 100, which is the interesting region. => Investigate and maybe fix |
| 293 | +# function _leftorth!(A::StridedMatrix{<:BlasFloat}) |
| 294 | +# m, n = size(A) |
| 295 | +# A, τ = geqrfp!(A) |
| 296 | +# Q = LAPACK.ormqr!('L', 'N', A, τ, eye(eltype(A), m, min(m, n))) |
| 297 | +# R = triu!(A[1:min(m, n), :]) |
| 298 | +# return Q, R |
| 299 | +# end |
325 | 300 |
|
326 |
| -for (geqrfp, elty, relty) in |
327 |
| - ((:dgeqrfp_, :Float64, :Float64), (:sgeqrfp_, :Float32, :Float32), |
328 |
| - (:zgeqrfp_, :ComplexF64, :Float64), (:cgeqrfp_, :ComplexF32, :Float32)) |
329 |
| - @eval begin |
330 |
| - function geqrfp!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}) |
331 |
| - chkstride1(A, tau) |
332 |
| - m, n = size(A) |
333 |
| - if length(tau) != min(m, n) |
334 |
| - throw(DimensionMismatch("tau has length $(length(tau)), but needs length $(min(m, n))")) |
335 |
| - end |
336 |
| - work = Vector{$elty}(1) |
337 |
| - lwork = BlasInt(-1) |
338 |
| - info = Ref{BlasInt}() |
339 |
| - for i = 1:2 # first call returns lwork as work[1] |
340 |
| - ccall((@blasfunc($geqrfp), liblapack), Nothing, |
341 |
| - (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, |
342 |
| - Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), |
343 |
| - Ref(m), Ref(n), A, Ref(max(1, stride(A, 2))), |
344 |
| - tau, work, Ref(lwork), info) |
345 |
| - chklapackerror(info[]) |
346 |
| - if i == 1 |
347 |
| - lwork = BlasInt(real(work[1])) |
348 |
| - resize!(work, lwork) |
349 |
| - end |
350 |
| - end |
351 |
| - A, tau |
352 |
| - end |
353 |
| - end |
354 |
| -end |
| 301 | +# geqrfp!: computes qrpos factorization, missing in Base |
| 302 | +# geqrfp!(A::StridedMatrix{<:BlasFloat}) = |
| 303 | +# ((m, n) = size(A); geqrfp!(A, similar(A, min(m, n)))) |
| 304 | +# |
| 305 | +# for (geqrfp, elty, relty) in |
| 306 | +# ((:dgeqrfp_, :Float64, :Float64), (:sgeqrfp_, :Float32, :Float32), |
| 307 | +# (:zgeqrfp_, :ComplexF64, :Float64), (:cgeqrfp_, :ComplexF32, :Float32)) |
| 308 | +# @eval begin |
| 309 | +# function geqrfp!(A::StridedMatrix{$elty}, tau::StridedVector{$elty}) |
| 310 | +# chkstride1(A, tau) |
| 311 | +# m, n = size(A) |
| 312 | +# if length(tau) != min(m, n) |
| 313 | +# throw(DimensionMismatch("tau has length $(length(tau)), but needs length $(min(m, n))")) |
| 314 | +# end |
| 315 | +# work = Vector{$elty}(1) |
| 316 | +# lwork = BlasInt(-1) |
| 317 | +# info = Ref{BlasInt}() |
| 318 | +# for i = 1:2 # first call returns lwork as work[1] |
| 319 | +# ccall((@blasfunc($geqrfp), liblapack), Nothing, |
| 320 | +# (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, |
| 321 | +# Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), |
| 322 | +# Ref(m), Ref(n), A, Ref(max(1, stride(A, 2))), |
| 323 | +# tau, work, Ref(lwork), info) |
| 324 | +# chklapackerror(info[]) |
| 325 | +# if i == 1 |
| 326 | +# lwork = BlasInt(real(work[1])) |
| 327 | +# resize!(work, lwork) |
| 328 | +# end |
| 329 | +# end |
| 330 | +# A, tau |
| 331 | +# end |
| 332 | +# end |
| 333 | +# end |
0 commit comments