diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 450349b5631bb..f9af90061ead0 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -111,7 +111,7 @@ immutable QRPackedQ{T} <: AbstractMatrix{T} factors::Matrix{T} τ::Vector{T} end -immutable QRCompactWYQ{S} <: AbstractMatrix{S} +immutable QRCompactWYQ{S} <: AbstractMatrix{S} factors::Matrix{S} T::Triangular{S} end @@ -445,7 +445,7 @@ end function eigfact!{T<:BlasComplex}(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true) n = size(A, 2) n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0)) - ishermitian(A) && return eigfact!(Hermitian(A)) + ishermitian(A) && return eigfact!(Hermitian(A)) return Eigen(LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]...) end eigfact{T}(A::AbstractMatrix{T}, args...; kwargs...) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigfact!(convert(AbstractMatrix{S}, A), args...; kwargs...) : eigfact!(copy(A), args...; kwargs...)) @@ -668,9 +668,9 @@ immutable GeneralizedSchur{Ty<:BlasFloat} <: Factorization{Ty} Z::Matrix{Ty} end -schurfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}) = GeneralizedSchur(LinAlg.LAPACK.gges!('V', 'V', A, B)...) -schurfact{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}) = schurfact!(copy(A),copy(B)) -schurfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); schurfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B))) +schurfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}, sort::Char='N', selctg::Ptr{None}=C_NULL) = GeneralizedSchur(LinAlg.LAPACK.gges!('V', 'V', A, B)...) +schurfact{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}, sort::Char='N', selctg::Ptr{None}=C_NULL) = schurfact!(copy(A),copy(B)) +schurfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}, sort::Char='N', selctg::Ptr{None}=C_NULL) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); schurfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B))) function getindex(F::GeneralizedSchur, d::Symbol) d == :S && return F.S @@ -692,16 +692,16 @@ end convert{T}(::Type{Factorization{T}}, F::Factorization{T}) = F inv{T}(F::Factorization{T}) = A_ldiv_B!(F, eye(T, size(F,1))) function \{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N}) - TFB = typeof(one(TF)/one(TB)) + TFB = typeof(one(TF)/one(TB)) A_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B)) end function Ac_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N}) - TFB = typeof(one(TF)/one(TB)) + TFB = typeof(one(TF)/one(TB)) Ac_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B)) end function At_ldiv_B{TF<:Number,TB<:Number,N}(F::Factorization{TF}, B::AbstractArray{TB,N}) - TFB = typeof(one(TF)/one(TB)) + TFB = typeof(one(TF)/one(TB)) At_ldiv_B!(convert(Factorization{TFB}, F), TB == TFB ? copy(B) : convert(AbstractArray{TFB,N}, B)) end diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 15be3b058c5d9..6be6c3963d750 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -260,16 +260,16 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}), - &m, &n, A, &lda, - jpvt, tau, work, &lwork, + &m, &n, A, &lda, + jpvt, tau, work, &lwork, rwork, info) else ccall(($(string(geqp3)),liblapack), Void, (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &m, &n, A, &lda, - jpvt, tau, work, + &m, &n, A, &lda, + jpvt, tau, work, &lwork, info) end @lapackerror @@ -290,11 +290,11 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty work = Array($elty, nb*n) if n > 0 info = Array(BlasInt, 1) - ccall(($(string(geqrt)), liblapack), Void, - (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, + ccall(($(string(geqrt)), liblapack), Void, + (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - &m, &n, &nb, A, + &m, &n, &nb, A, &lda, T, &max(1,stride(T,2)), work, info) @lapackerror @@ -308,7 +308,7 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, geqrt, geqrt3, gerqf, getrf, elty, relty if p < n || q < n throw(DimensionMismatch("block reflector")) end if n > 0 info = Array(BlasInt, 1) - ccall(($(string(geqrt3)), liblapack), Void, + ccall(($(string(geqrt3)), liblapack), Void, (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), &m, &n, A, &max(1, stride(A, 2)), @@ -407,14 +407,14 @@ function geqp3!{T<:BlasFloat}(A::StridedMatrix{T}) end ## Complete orthogonaliztion tools -for (tzrzf, ormrz, elty) in +for (tzrzf, ormrz, elty) in ((:dtzrzf_,:dormrz_,:Float64), (:stzrzf_,:sormrz_,:Float32), (:ztzrzf_,:zunmrz_,:Complex128), (:ctzrzf_,:cunmrz_,:Complex64)) @eval begin # * SUBROUTINE ZTZRZF( M, N, A, LDA, TAU, WORK, LWORK, INFO ) - # 22 * + # 22 * # 23 * .. Scalar Arguments .. # 24 * INTEGER INFO, LDA, LWORK, M, N # 25 * .. @@ -432,7 +432,7 @@ for (tzrzf, ormrz, elty) in ccall(($(string(tzrzf)), liblapack), Void, (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &m, &n, A, &lda, + &m, &n, A, &lda, tau, work, &lwork, info) if i == 1 lwork = blas_int(real(work[1])) @@ -444,7 +444,7 @@ for (tzrzf, ormrz, elty) in end # 21 * SUBROUTINE ZUNMRZ( SIDE, TRANS, M, N, K, L, A, LDA, TAU, C, LDC, # 22 * WORK, LWORK, INFO ) - # 23 * + # 23 * # 24 * .. Scalar Arguments .. # 25 * CHARACTER SIDE, TRANS # 26 * INTEGER INFO, K, L, LDA, LDC, LWORK, M, N @@ -466,9 +466,9 @@ for (tzrzf, ormrz, elty) in Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &side, &trans, &m, &n, - &k, &l, A, &lda, - tau, C, &ldc, work, + &side, &trans, &m, &n, + &k, &l, A, &lda, + tau, C, &ldc, work, &lwork, info) if i == 1 lwork = blas_int(real(work[1])) @@ -601,7 +601,7 @@ for (gesvx, elty) in # SUBROUTINE DGESVX( FACT, TRANS, N, NRHS, A, LDA, AF, LDAF, IPIV, # EQUED, R, C, B, LDB, X, LDX, RCOND, FERR, BERR, # WORK, IWORK, INFO ) - # + # # .. Scalar Arguments .. # CHARACTER EQUED, FACT, TRANS # INTEGER INFO, LDA, LDAF, LDB, LDX, N, NRHS @@ -611,7 +611,7 @@ for (gesvx, elty) in # INTEGER IPIV( * ), IWORK( * ) # DOUBLE PRECISION A( LDA, * ), AF( LDAF, * ), B( LDB, * ), # $ BERR( * ), C( * ), FERR( * ), R( * ), - # $ WORK( * ), X( LDX, * + # $ WORK( * ), X( LDX, * # function gesvx!(fact::BlasChar, trans::BlasChar, A::StridedMatrix{$elty}, AF::StridedMatrix{$elty}, ipiv::Vector{BlasInt}, equed::BlasChar, @@ -635,7 +635,7 @@ for (gesvx, elty) in &fact, &trans, &n, &nrhs, A, &lda, AF, &ldaf, ipiv, &equed, R, C, B, &ldb, X, &n, rcond, ferr, berr, work, iwork, info) @lapackerror - if info[1] == n+1 warn("Matrix is singular to working precision.") + if info[1] == n+1 warn("Matrix is singular to working precision.") else @assertnonsingular end #WORK(1) contains the reciprocal pivot growth factor norm(A)/norm(U) @@ -704,7 +704,7 @@ for (gesvx, elty, relty) in end end -for (gelsd, gelsy, elty) in +for (gelsd, gelsy, elty) in ((:dgelsd_,:dgelsy_,:Float64), (:sgelsd_,:sgelsy_,:Float32)) @eval begin @@ -772,14 +772,14 @@ for (gelsd, gelsy, elty) in lwork = -1 info = Array(BlasInt, 1) for i = 1:2 - ccall(($(string(gelsy)), liblapack), Void, + ccall(($(string(gelsy)), liblapack), Void, (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &m, &n, &nrhs, A, - &lda, newB, &ldb, jpvt, - &rcond, rnk, work, &lwork, + &m, &n, &nrhs, A, + &lda, newB, &ldb, jpvt, + &rcond, rnk, work, &lwork, info) if i == 1 lwork = blas_int(work[1]) @@ -791,7 +791,7 @@ for (gelsd, gelsy, elty) in end end end -for (gelsd, gelsy, elty, relty) in +for (gelsd, gelsy, elty, relty) in ((:zgelsd_,:zgelsy_,:Complex128,:Float64), (:cgelsd_,:cgelsy_,:Complex64,:Float32)) @eval begin @@ -863,14 +863,14 @@ for (gelsd, gelsy, elty, relty) in rwork = Array($relty, 2n) info = Array(BlasInt, 1) for i = 1:2 - ccall(($(string(gelsy)), liblapack), Void, + ccall(($(string(gelsy)), liblapack), Void, (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}), - &m, &n, &nrhs, A, - &lda, newB, &ldb, jpvt, - &rcond, rnk, work, &lwork, + &m, &n, &nrhs, A, + &lda, newB, &ldb, jpvt, + &rcond, rnk, work, &lwork, rwork, info) if i == 1 lwork = blas_int(real(work[1])) @@ -964,7 +964,7 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in if cmplx ccall(($(string(geev)),liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, + Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}), &jobvl, &jobvr, &n, A, &max(1,stride(A,2)), W, VL, &n, VR, &n, @@ -1136,11 +1136,11 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}), - &jobu, &jobv, &jobq, &m, - &n, &p, k, l, - A, &lda, B, &ldb, - alpha, beta, U, &ldu, - V, &ldv, Q, &ldq, + &jobu, &jobv, &jobq, &m, + &n, &p, k, l, + A, &lda, B, &ldb, + alpha, beta, U, &ldu, + V, &ldv, Q, &ldq, work, rwork, iwork, info) else ccall(($(string(ggsvd)),liblapack), Void, @@ -1150,11 +1150,11 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &jobu, &jobv, &jobq, &m, - &n, &p, k, l, - A, &lda, B, &ldb, - alpha, beta, U, &ldu, - V, &ldv, Q, &ldq, + &jobu, &jobv, &jobq, &m, + &n, &p, k, l, + A, &lda, B, &ldb, + alpha, beta, U, &ldu, + V, &ldv, Q, &ldq, work, iwork, info) end @lapackerror @@ -1168,14 +1168,14 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in end end ## Expert driver and generalized eigenvlua problem -for (geevx, ggev, elty) in +for (geevx, ggev, elty) in ((:dgeevx_,:dggev_,:Float64), (:sgeevx_,:sggev_,:Float32)) @eval begin # SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, # VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, # RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) - # + # # .. Scalar Arguments .. # CHARACTER BALANC, JOBVL, JOBVR, SENSE # INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N @@ -1213,11 +1213,11 @@ for (geevx, ggev, elty) in Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), - &balanc, &jobvl, &jobvr, &sense, - &n, A, &lda, wr, - wi, VL, &max(1,ldvl), VR, - &max(1,ldvr), ilo, ihi, scale, - abnrm, rconde, rcondv, work, + &balanc, &jobvl, &jobvr, &sense, + &n, A, &lda, wr, + wi, VL, &max(1,ldvl), VR, + &max(1,ldvr), ilo, ihi, scale, + abnrm, rconde, rcondv, work, &lwork, iwork, info) lwork = convert(BlasInt, work[1]) work = Array($elty, lwork) @@ -1259,10 +1259,10 @@ for (geevx, ggev, elty) in Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &jobvl, &jobvr, &n, A, - &lda, B, &ldb, alphar, - alphai, beta, vl, &ldvl, - vr, &ldvr, work, &lwork, + &jobvl, &jobvr, &n, A, + &lda, B, &ldb, alphar, + alphai, beta, vl, &ldvl, + vr, &ldvr, work, &lwork, info) if i == 1 lwork = blas_int(work[1]) @@ -1274,14 +1274,14 @@ for (geevx, ggev, elty) in end end end -for (geevx, ggev, elty, relty) in +for (geevx, ggev, elty, relty) in ((:zgeevx_,:zggev_,:Complex128,:Float64), (:cgeevx_,:cggev_,:Complex64,:Float32)) @eval begin # SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL, # LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, # RCONDV, WORK, LWORK, RWORK, INFO ) - # + # # .. Scalar Arguments .. # CHARACTER BALANC, JOBVL, JOBVR, SENSE # INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N @@ -1314,15 +1314,15 @@ for (geevx, ggev, elty, relty) in ccall(($(string(geevx)),Base.liblapack_name), Void, (Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, - Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty}, - Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, + Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty}, + Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}), - &balanc, &jobvl, &jobvr, &sense, - &n, A, &lda, w, - VL, &max(1,ldvl), VR, &max(1,ldvr), - ilo, ihi, scale, abnrm, - rconde, rcondv, work, &lwork, + &balanc, &jobvl, &jobvr, &sense, + &n, A, &lda, w, + VL, &max(1,ldvl), VR, &max(1,ldvr), + ilo, ihi, scale, abnrm, + rconde, rcondv, work, &lwork, rwork, info) lwork = convert(BlasInt, work[1]) work = Array($elty, lwork) @@ -1361,12 +1361,12 @@ for (geevx, ggev, elty, relty) in ccall(($(string(ggev)), liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, + Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}), - &jobvl, &jobvr, &n, A, - &lda, B, &ldb, alpha, - beta, vl, &ldvl, vr, + &jobvl, &jobvr, &n, A, + &lda, B, &ldb, alpha, + beta, vl, &ldvl, vr, &ldvr, work, &lwork, rwork, info) if i == 1 @@ -1386,7 +1386,7 @@ for (laic1, elty) in (:slaic1_,:Float32)) @eval begin # 21 * SUBROUTINE DLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) - # 22 * + # 22 * # 23 * .. Scalar Arguments .. # 24 * INTEGER J, JOB # 25 * DOUBLE PRECISION C, GAMMA, S, SEST, SESTPR @@ -1403,8 +1403,8 @@ for (laic1, elty) in (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}), - &job, &j, x, &sest, - w, &gamma, sestpr, s, + &job, &j, x, &sest, + w, &gamma, sestpr, s, c) sestpr[1], s[1], c[1] end @@ -1415,7 +1415,7 @@ for (laic1, elty, relty) in (:claic1_,:Complex64,:Float32)) @eval begin # 21 * SUBROUTINE ZLAIC1( JOB, J, X, SEST, W, GAMMA, SESTPR, S, C ) - # 22 * + # 22 * # 23 * .. Scalar Arguments .. # 24 * INTEGER J, JOB # 25 * DOUBLE PRECISION SEST, SESTPR @@ -1433,8 +1433,8 @@ for (laic1, elty, relty) in (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$relty}, Ptr{$elty}, Ptr{$elty}, Ptr{$relty}, Ptr{$elty}, Ptr{$elty}), - &job, &j, x, &sest, - w, &gamma, sestpr, s, + &job, &j, x, &sest, + w, &gamma, sestpr, s, c) sestpr[1], s[1], c[1] end @@ -1447,7 +1447,7 @@ for (gtsv, gttrf, gttrs, elty) in ((:dgtsv_,:dgttrf_,:dgttrs_,:Float64), (:sgtsv_,:sgttrf_,:sgttrs_,:Float32), (:zgtsv_,:zgttrf_,:zgttrs_,:Complex128), - (:cgtsv_,:cgttrf_,:cgttrs_,:Complex64)) + (:cgtsv_,:cgttrf_,:cgttrs_,:Complex64)) @eval begin # SUBROUTINE DGTSV( N, NRHS, DL, D, DU, B, LDB, INFO ) # .. Scalar Arguments .. @@ -1543,7 +1543,7 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), &m, &n, &k, A, &max(1,stride(A,2)), tau, work, &lwork, info) @lapackerror - if lwork < 0 + if lwork < 0 lwork = blas_int(real(work[1])) work = Array($elty, lwork) end @@ -1571,11 +1571,11 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in ccall(($(string(orgqr)),liblapack), Void, (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &m, &n, &k, A, - &max(1,stride(A,2)), tau, work, &lwork, + &m, &n, &k, A, + &max(1,stride(A,2)), tau, work, &lwork, info) @lapackerror - if lwork < 0 + if lwork < 0 lwork = blas_int(real(work[1])) work = Array($elty, lwork) end @@ -1600,7 +1600,7 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in nA = size(A, 2) k = length(tau) if side == 'L' && m != nA throw(DimensionMismatch("")) end - if side == 'R' && n != nA throw(DimensionMismatch("")) end + if side == 'R' && n != nA throw(DimensionMismatch("")) end if (side == 'L' && k > m) || (side == 'R' && k > n) throw(DimensionMismatch("invalid number of reflectors")) end work = Array($elty, 1) lwork = blas_int(-1) @@ -1613,7 +1613,7 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in &side, &trans, &m, &n, &k, A, &max(1,stride(A,2)), tau, C, &max(1,stride(C,2)), work, &lwork, info) @lapackerror - if lwork < 0 + if lwork < 0 lwork = blas_int(real(work[1])) work = Array($elty, lwork) end @@ -1634,23 +1634,23 @@ for (orglq, orgqr, ormlq, ormqr, gemqrt, elty) in mA = size(A, 1) k = length(tau) if side == 'L' && m != mA throw(DimensionMismatch("")) end - if side == 'R' && n != mA throw(DimensionMismatch("")) end + if side == 'R' && n != mA throw(DimensionMismatch("")) end if (side == 'L' && k > m) || (side == 'R' && k > n) throw(DimensionMismatch("invalid number of reflectors")) end work = Array($elty, 1) lwork = blas_int(-1) info = Array(BlasInt, 1) for i in 1:2 ccall(($(string(ormqr)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, - Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, + Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, + Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &side, &trans, &m, &n, + &side, &trans, &m, &n, &k, A, &max(1,stride(A,2)), tau, - C, &max(1, stride(C,2)), work, &lwork, + C, &max(1, stride(C,2)), work, &lwork, info) @lapackerror - if lwork < 0 + if lwork < 0 lwork = blas_int(real(work[1])) work = Array($elty, lwork) end @@ -1778,9 +1778,9 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in ldb = max(1,stride(B,2)) info = Array(BlasInt, 1) ccall(($(string(potrs)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, + (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &uplo, &n, &nrhs, A, + &uplo, &n, &nrhs, A, &lda, B, &ldb, info) @lapackerror return B @@ -1815,8 +1815,8 @@ end ## Direct solvers for general tridiagonal and symmetric positive-definite tridiagonal for (ptsv, pttrf, pttrs, elty, relty) in ((:dptsv_,:dpttrf_,:dpttrs_,:Float64,:Float64), - (:sptsv_,:spttrf_,:spttrs_,:Float32,:Float32), - (:zptsv_,:zpttrf_,:zpttrs_,:Complex128,:Float64), + (:sptsv_,:spttrf_,:spttrs_,:Float32,:Float32), + (:zptsv_,:zpttrf_,:zpttrs_,:Complex128,:Float64), (:cptsv_,:cpttrf_,:cpttrs_,:Complex64,:Float32)) @eval begin # SUBROUTINE DPTSV( N, NRHS, D, E, B, LDB, INFO ) @@ -2054,7 +2054,7 @@ for (trcon, trevc, trrfs, elty) in iwork=Array(BlasInt, n) info=Array(BlasInt, 1) ccall(($(string(trrfs)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), &uplo, &trans, &diag, &n, @@ -2168,8 +2168,8 @@ for (trcon, trevc, trrfs, elty, relty) in rwork=Array($relty, n) info=Array(BlasInt, 1) ccall(($(string(trrfs)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, - Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, + Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}), &uplo, &trans, &diag, &n, &nrhs, A, &max(1,stride(A,2)), B, &max(1,stride(B,2)), X, &max(1,stride(X,2)), @@ -2200,7 +2200,7 @@ for (stev, stebz, stegr, stein, elty) in (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), &job, &n, dv, ev, Zmat, &n, work, info) - @lapackerror + @lapackerror dv, Zmat end #* DSTEBZ computes the eigenvalues of a symmetric tridiagonal @@ -2225,10 +2225,10 @@ for (stev, stebz, stegr, stein, elty) in Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &range, &order, &n, &vl, - &vu, &il, &iu, &abstol, + &range, &order, &n, &vl, + &vu, &il, &iu, &abstol, dv, ev, m, nsplit, - w, iblock, isplit, work, + w, iblock, isplit, work, iwork, info) @lapackerror w[1:m[1]], iblock[1:m[1]], isplit[1:nsplit[1]] @@ -2257,7 +2257,7 @@ for (stev, stebz, stegr, stein, elty) in liwork = -one(BlasInt) info = Array(BlasInt, 1) for i = 1:2 - ccall(($(string(stegr)), liblapack), Void, + ccall(($(string(stegr)), liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, @@ -2322,7 +2322,7 @@ for (stev, stebz, stegr, stein, elty) in Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), &n, dv, ev, &m, w, iblock, isplit, z, &ldz, work, iwork, ifail, info) - + @lapackerror all(ifail.==0) || error("failed to converge eigenvectors:\n$(nonzeros(ifail))") z @@ -2330,7 +2330,7 @@ for (stev, stebz, stegr, stein, elty) in end end stegr!(jobz::BlasChar, dv::Vector, ev::Vector) = stegr!(jobz, 'A', dv, ev, 0.0, 0.0, 0, 0) - + # Allow user to skip specification of iblock and isplit stein!(dv::Vector, ev::Vector, w_in::Vector)=stein!(dv, ev, w_in, zeros(BlasInt,0), zeros(BlasInt,0)) # Allow user to specify just one eigenvector to get in stein! @@ -2474,7 +2474,7 @@ for (syconv, sysv, sytrf, sytri, sytrs, elty) in A end # SUBROUTINE DSYTRS( UPLO, N, NRHS, A, LDA, IPIV, B, LDB, INFO ) - # + # # .. Scalar Arguments .. # CHARACTER UPLO # INTEGER INFO, LDA, LDB, N, NRHS @@ -2504,7 +2504,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in (:csyconv_,:chesv_,:chetrf_,:chetri_,:chetrs_,:Complex64, :Float32)) @eval begin # SUBROUTINE ZSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO ) - # 22 * + # 22 * # 23 * .. Scalar Arguments .. # 24 * CHARACTER UPLO, WAY # 25 * INTEGER INFO, LDA, N @@ -2653,7 +2653,7 @@ for (syconv, hesv, hetrf, hetri, hetrs, elty, relty) in (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), &uplo, &n, &size(B,2), A, &max(1,stride(A,2)), ipiv, B, &max(1,stride(B,2)), info) - @lapackerror + @lapackerror B end end @@ -2840,14 +2840,14 @@ for (syev, syevr, sygvd, elty) in # * .. # * .. Array Arguments .. # INTEGER ISUPPZ( * ), IWORK( * ) - # DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * ) + # DOUBLE PRECISION A( LDA, * ), W( * ), WORK( * ), Z( LDZ, * ) function syevr!(jobz::BlasChar, range::BlasChar, uplo::BlasChar, A::StridedMatrix{$elty}, vl::FloatingPoint, vu::FloatingPoint, il::Integer, iu::Integer, abstol::FloatingPoint) chkstride1(A) n = chksquare(A) if range == 'I' 1 <= il <= iu <= n || throw(ArgumentError("illegal choice of eigenvalue indices")) end - if range == 'V' + if range == 'V' vl < vu || throw(ArgumentError("lower boundary must be less than upper boundary")) end lda = max(1,stride(A,2)) @@ -2868,17 +2868,17 @@ for (syev, syevr, sygvd, elty) in info = Array(BlasInt, 1) for i in 1:2 ccall(($(string(syevr)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, - Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, + Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), - &jobz, &range, &uplo, &n, - A, &lda, &vl, &vu, + &jobz, &range, &uplo, &n, + A, &lda, &vl, &vu, &il, &iu, &abstol, m, w, Z, &ldz, isuppz, - work, &lwork, iwork, &liwork, + work, &lwork, iwork, &liwork, info) @lapackerror if lwork < 0 @@ -2889,8 +2889,8 @@ for (syev, syevr, sygvd, elty) in end end w[1:m[1]], Z[:,1:(jobz == 'V' ? m[1] : 0)] - end - syevr!(jobz::BlasChar, A::StridedMatrix{$elty}) = syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0) + end + syevr!(jobz::BlasChar, A::StridedMatrix{$elty}) = syevr!(jobz, 'A', 'U', A, 0.0, 0.0, 0, 0, -1.0) # Generalized eigenproblem # SUBROUTINE DSYGVD( ITYPE, JOBZ, UPLO, N, A, LDA, B, LDB, W, WORK, # $ LWORK, IWORK, LIWORK, INFO ) @@ -2918,9 +2918,9 @@ for (syev, syevr, sygvd, elty) in Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), - &itype, &jobz, &uplo, &n, - A, &lda, B, &ldb, - w, work, &lwork, iwork, + &itype, &jobz, &uplo, &n, + A, &lda, B, &ldb, + w, work, &lwork, iwork, &liwork, info) if i == 1 lwork = blas_int(work[1]) @@ -2936,11 +2936,11 @@ for (syev, syevr, sygvd, elty) in end end # Hermitian eigensolvers -for (syev, syevr, sygvd, elty, relty) in +for (syev, syevr, sygvd, elty, relty) in ((:zheev_,:zheevr_,:zhegvd_,:Complex128,:Float64), (:cheev_,:cheevr_,:chegvd_,:Complex64,:Float32)) @eval begin -# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO ) +# SUBROUTINE ZHEEV( JOBZ, UPLO, N, A, LDA, W, WORK, LWORK, RWORK, INFO ) # * .. Scalar Arguments .. # CHARACTER JOBZ, UPLO # INTEGER INFO, LDA, LWORK, N @@ -2981,14 +2981,14 @@ for (syev, syevr, sygvd, elty, relty) in # * .. Array Arguments .. # INTEGER ISUPPZ( * ), IWORK( * ) # DOUBLE PRECISION RWORK( * ), W( * ) -# COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * ) +# COMPLEX*16 A( LDA, * ), WORK( * ), Z( LDZ, * ) function syevr!(jobz::BlasChar, range::BlasChar, uplo::BlasChar, A::StridedMatrix{$elty}, vl::FloatingPoint, vu::FloatingPoint, il::Integer, iu::Integer, abstol::FloatingPoint) chkstride1(A) n = chksquare(A) if range == 'I' 1 <= il <= iu <= n || throw(ArgumentError("illegal choice of eigenvalue indices")) end - if range == 'V' + if range == 'V' vl < vu || throw(ArgumentError("lower boundary must be less than upper boundary")) end lda = max(1,stride(A,2)) @@ -3011,14 +3011,14 @@ for (syev, syevr, sygvd, elty, relty) in info = Array(BlasInt, 1) for i in 1:2 ccall(($(string(syevr)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, - Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, + Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), - &jobz, &range, &uplo, &n, - A, &lda, &vl, &vu, + &jobz, &range, &uplo, &n, + A, &lda, &vl, &vu, &il, &iu, &abstol, m, w, Z, &ldz, isuppz, work, &lwork, rwork, &lrwork, @@ -3071,9 +3071,9 @@ for (syev, syevr, sygvd, elty, relty) in Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}), - &itype, &jobz, &uplo, &n, - A, &lda, B, &ldb, - w, work, &lwork, rwork, + &itype, &jobz, &uplo, &n, + A, &lda, B, &ldb, + w, work, &lwork, rwork, &lrwork, iwork, &liwork, info) if i == 1 lwork = blas_int(real(work[1])) @@ -3115,7 +3115,7 @@ for (bdsqr, relty, elty) in ccall(($(string(bdsqr)),liblapack), Void, (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), + Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), &uplo, &n, ncvt, &nru, &ncc, d, e_, vt, &ldvt, u, &ldu, c, &ldc, work, info) @@ -3177,7 +3177,7 @@ for (bdsdc, elty) in q, iq, work, iwork, info) @lapackerror - compq=='N' ? d : (compq=='P' ? (d, q, iq) : (u, d, vt')) + compq=='N' ? d : (compq=='P' ? (d, q, iq) : (u, d, vt')) end end end @@ -3206,7 +3206,7 @@ for (gecon, elty) in iwork = Array(BlasInt, n) info = Array(BlasInt, 1) ccall(($(string(gecon)),liblapack), Void, - (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), &normtype, &n, A, &lda, &anorm, rcond, work, iwork, @@ -3240,7 +3240,7 @@ for (gecon, elty, relty) in rwork = Array($relty, 2n) info = Array(BlasInt, 1) ccall(($(string(gecon)),liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty}, Ptr{BlasInt}), &normtype, &n, A, &lda, &anorm, rcond, work, rwork, @@ -3275,8 +3275,8 @@ for (gehrd, elty) in (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &n, &ilo, &ihi, A, - &max(1,n), tau, work, &lwork, + &n, &ilo, &ihi, A, + &max(1,n), tau, work, &lwork, info) @lapackerror if lwork < 0 @@ -3314,8 +3314,8 @@ for (orghr, elty) in (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}), - &n, &ilo, &ihi, A, - &max(1,n), tau, work, &lwork, + &n, &ilo, &ihi, A, + &max(1,n), tau, work, &lwork, info) @lapackerror if lwork < 0 @@ -3357,9 +3357,9 @@ for (gees, gges, elty) in Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, Ptr{BlasInt}), - &jobvs, &'N', C_NULL, &n, + &jobvs, &'N', C_NULL, &n, A, &max(1, n), sdim, wr, - wi, vs, &ldvs, work, + wi, vs, &ldvs, work, &lwork, C_NULL, info) @lapackerror if lwork < 0 @@ -3369,7 +3369,9 @@ for (gees, gges, elty) in end A, vs, all(wi .== 0) ? wr : complex(wr, wi) end - function gges!(jobvsl::Char, jobvsr::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function gges!(jobvsl::Char, jobvsr::Char, + A::StridedMatrix{$elty}, B::StridedMatrix{$elty}), + sort::Char, selctg::Ptr{None}) # * .. Scalar Arguments .. # CHARACTER JOBVSL, JOBVSR, SORT # INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM @@ -3401,11 +3403,11 @@ for (gees, gges, elty) in Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, Ptr{BlasInt}), - &jobvsl, &jobvsr, &'N', C_NULL, - &n, A, &max(1,n), B, - &max(1,n), &sdim, alphar, alphai, - beta, vsl, &ldvsl, vsr, - &ldvsr, work, &lwork, C_NULL, + &jobvsl, &jobvsr, &sort, selctg, + &n, A, &max(1,n), B, + &max(1,n), &sdim, alphar, alphai, + beta, vsl, &ldvsl, vsr, + &ldvsr, work, &lwork, C_NULL, info) if i == 1 lwork = blas_int(real(work[1])) @@ -3445,11 +3447,11 @@ for (gees, gges, elty, relty) in ccall(($(string(gees)),liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Void}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, - Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, + Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{Void}, Ptr{BlasInt}), - &jobvs, &sort, C_NULL, &n, + &jobvs, &sort, C_NULL, &n, A, &max(1, n), &sdim, w, - vs, &ldvs, work, &lwork, + vs, &ldvs, work, &lwork, rwork, C_NULL, info) @lapackerror if lwork < 0 @@ -3459,7 +3461,9 @@ for (gees, gges, elty, relty) in end A, vs, w end - function gges!(jobvsl::Char, jobvsr::Char, A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + function gges!(jobvsl::Char, jobvsr::Char, + A::StridedMatrix{$elty}, B::StridedMatrix{$elty}) + sort::Char, selctg::Ptr{None} # * .. Scalar Arguments .. # CHARACTER JOBVSL, JOBVSR, SORT # INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM @@ -3492,11 +3496,11 @@ for (gees, gges, elty, relty) in Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{Void}, Ptr{BlasInt}), - &jobvsl, &jobvsr, &'N', C_NULL, - &n, A, &max(1,n), B, - &max(1,n), &sdim, alpha, beta, - vsl, &ldvsl, vsr, &ldvsr, - work, &lwork, rwork, C_NULL, + &jobvsl, &jobvsr, &sort, selctg, + &n, A, &max(1,n), B, + &max(1,n), &sdim, alpha, beta, + vsl, &ldvsl, vsr, &ldvsr, + work, &lwork, rwork, C_NULL, info) if i == 1 lwork = blas_int(real(work[1])) @@ -3570,7 +3574,7 @@ for (fn, elty) in ((:dpftri_, :Float64), ccall(($(string(fn)), liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - &transr, &uplo, &n, A, + &transr, &uplo, &n, A, info) @assertargsok @assertnonsingular @@ -3637,9 +3641,9 @@ for (fn, elty) in ((:dtftri_, :Float64), n = int(div(sqrt(8length(A)), 2)) info = Array(BlasInt, 1) ccall(($(string(fn)), liblapack), Void, - (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - &transr, &uplo, &diag, &n, + &transr, &uplo, &diag, &n, A, info) @assertargsok @assertnonsingular @@ -3709,15 +3713,15 @@ for (fn, elty, relty) in ((:dtrsyl_, :Float64, :Float64), scale = Array($relty, 1) info = Array(BlasInt, 1) - - ccall(($(string(fn)), liblapack), Void, + + ccall(($(string(fn)), liblapack), Void, (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}), &transa, &transb, &isgn, &m, &n, A, &lda, B, &ldb, C, &ldc, scale, info) - @lapackerror + @lapackerror C, scale[1] end end diff --git a/base/linalg/test_sort.jl b/base/linalg/test_sort.jl new file mode 100644 index 0000000000000..b9c5150b006f6 --- /dev/null +++ b/base/linalg/test_sort.jl @@ -0,0 +1,116 @@ +const liblapack = Base.liblapack_name + +using Base, Base.LinAlg + +import Base.LinAlg: BlasFloat, BlasChar, BlasInt, blas_int, LAPACKException, + DimensionMismatch, SingularException, PosDefException, chkstride1, chksquare + +#Generic LAPACK error handlers +macro assertargsok() #Handle only negative info codes - use only if positive info code is useful! + :(info[1]<0 && throw(ArgumentError("invalid argument #$(-info[1]) to LAPACK call"))) +end +macro lapackerror() #Handle all nonzero info codes + :(info[1]>0 ? throw(LAPACKException(info[1])) : @assertargsok ) +end + +macro assertnonsingular() + :(info[1]>0 && throw(SingularException(info[1]))) +end +macro assertposdef() + :(info[1]>0 && throw(PosDefException(info[1]))) +end + +#Check that upper/lower (for special matrices) is correctly specified +macro chkuplo() + :((uplo=='U' || uplo=='L') || throw(ArgumentError("""invalid uplo = $uplo + +Valid choices are 'U' (upper) or 'L' (lower)."""))) +end + + +(gges, elty) = (:dgges_, :Float64) +@eval begin + function test_gges!(jobvsl::Char, jobvsr::Char, + A::StridedMatrix{$elty}, B::StridedMatrix{$elty}, + sort::Char, selctg::Ptr{None}) + + # * .. Scalar Arguments .. + # CHARACTER JOBVSL, JOBVSR, SORT + # INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM + # * .. + # * .. Array Arguments .. + # LOGICAL BWORK( * ) + # DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), + # $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ), + # $ VSR( LDVSR, * ), WORK( * ) + + chkstride1(A, B) + n, m = chksquare(A, B) + n==m || throw(DimensionMismatch("Matrices are not of the same size")) + sdim = blas_int(0) + alphar = similar(A, $elty, n) + alphai = similar(A, $elty, n) + beta = similar(A, $elty, n) + ldvsl = jobvsl == 'V' ? n : 1 + vsl = similar(A, $elty, ldvsl, n) + ldvsr = jobvsr == 'V' ? n : 1 + vsr = similar(A, $elty, ldvsr, n) + work = Array($elty, 1) + lwork = blas_int(-1) + info = Array(BlasInt, 1) + for i = 1:2 + ccall(($(string(gges)), liblapack), Void, + (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Void}, + Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, + Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, + Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, + Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Void}, + Ptr{BlasInt}), + &jobvsl, &jobvsr, &sort, selctg, + &n, A, &max(1,n), B, + &max(1,n), &sdim, alphar, alphai, + beta, vsl, &ldvsl, vsr, + &ldvsr, work, &lwork, C_NULL, + info) + if i == 1 + lwork = blas_int(real(work[1])) + work = Array($elty, lwork) + end + end + @lapackerror + A, B, complex(alphar, alphai), beta, vsl[1:(jobvsl == 'V' ? n : 0),:], vsr[1:(jobvsr == 'V' ? n : 0),:] + end +end + + +test_schurfact{T<:BlasFloat}(A::StridedMatrix{T},B::StridedMatrix{T}, sort::Char='V', selctg::Ptr{None}=C_NULL) = schurfact!(copy(A),copy(B), sort, selctg) +test_schurfact!{T<:BlasFloat}(A::StridedMatrix{T}, B::StridedMatrix{T}, sort::Char='V', selctg::Ptr{None}=C_NULL) = GeneralizedSchur(test_gges!('V', 'V', A, B, sort, selctg)...) +test_schurfact{TA,TB}(A::StridedMatrix{TA}, B::StridedMatrix{TB}, sort::Char='V', selctg::Ptr{None}=C_NULL) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); test_schurfact!(S != TA ? convert(AbstractMatrix{S},A) : copy(A), S != TB ? convert(AbstractMatrix{S},B) : copy(B), sort, selctg)) + + +function order_eigs(alphar, alphai, beta) + if abs(alphar + alphai)/beta > 1 + return true + else + return false + end +end + +const order_eigs_c = cfunction(order_eigs, Bool, (Cdouble, Cdouble, Cdouble)) + +# Test call +A = randn(3, 3) +B = reshape(1:9, 3, 3) + +their_schur = schurfact(A, B) +our_schur = test_schurfact(A, B, 'N') + + +println(their_schur[:S] - our_schur[:S]) +println(their_schur[:T] - our_schur[:T]) +println(their_schur[:alpha] - our_schur[:alpha]) +println(their_schur[:beta] - our_schur[:beta]) +println(their_schur[:Q] - our_schur[:Q]) +println(their_schur[:Z] - our_schur[:Z]) + +test_schurfact(A, B, 'S', order_eigs_c) \ No newline at end of file