Skip to content

Commit

Permalink
Use $relty as a macro parameter for the type of RWORK in LAPACK inter…
Browse files Browse the repository at this point in the history 
…face.
  • Loading branch information
@ViralBShah
ViralBShah committed Jan 17, 2013
1 parent 4f71aa3 commit d3a8dd1
Showing 1 changed file with 24 additions and 29 deletions.
53 changes: 24 additions & 29 deletions base/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -143,11 +143,11 @@ end
# gegp3 - pivoted QR decomposition
# gerqf - unpivoted RQ decomposition
# getrf - LU decomposition
for (gebrd, gelqf, geqlf, geqrf, geqp3, gerqf, getrf, elty) in
((:dgebrd_,:dgelqf_,:dgeqlf_,:dgeqrf_,:dgeqp3_,:dgerqf_,:dgetrf_,:Float64),
(:sgebrd_,:sgelqf_,:sgeqlf_,:sgeqrf_,:sgeqp3_,:sgerqf_,:sgetrf_,:Float32),
(:zgebrd_,:zgelqf_,:zgeqlf_,:zgeqrf_,:zgeqp3_,:zgerqf_,:zgetrf_,:Complex128),
(:cgebrd_,:cgelqf_,:cgeqlf_,:cgeqrf_,:cgeqp3_,:cgerqf_,:cgetrf_,:Complex64))
for (gebrd, gelqf, geqlf, geqrf, geqp3, gerqf, getrf, elty, relty) in
((:dgebrd_,:dgelqf_,:dgeqlf_,:dgeqrf_,:dgeqp3_,:dgerqf_,:dgetrf_,:Float64, Float64),
(:sgebrd_,:sgelqf_,:sgeqlf_,:sgeqrf_,:sgeqp3_,:sgerqf_,:sgetrf_,:Float32, Float32),
(:zgebrd_,:zgelqf_,:zgeqlf_,:zgeqrf_,:zgeqp3_,:zgerqf_,:zgetrf_,:Complex128, Float64),
(:cgebrd_,:cgelqf_,:cgeqlf_,:cgeqrf_,:cgeqp3_,:cgerqf_,:cgetrf_,:Complex64, Float32))

This comment has been minimized.

Sign in to view

Copy link
@kmsquire

kmsquire Jan 17, 2013

Member

Curious why relty is defined as a type, while elty is defined as a symbol. At

julia/base/lapack.jl

Line 90 in d3a8dd1

((:dgebal_, :dgebak_, :Float64, :Float64),
, both are defined as symbols.

This comment has been minimized.

Sign in to view

Copy link
@ViralBShah

ViralBShah Jan 17, 2013

Author Member

I see that it is not consistent. I am just converting them all to symbols.
@eval begin
# SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
# INFO )
Expand Down Expand Up @@ -248,15 +248,14 @@ for (gebrd, gelqf, geqlf, geqrf, geqp3, gerqf, getrf, elty) in
work = Array($elty, 1)
lwork = blas_int(-1)
info = Array(BlasInt, 1)
Rtyp = typeof(real(A[1]))
cmplx = iscomplex(A)
if cmplx rwork = Array(Rtyp, 2n) end
if cmplx; rwork = Array($relty, 2n); end
for i in 1:2
if cmplx
ccall(($(string(geqp3)),liblapack), Void,
(Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
Ptr{Rtyp}, Ptr{BlasInt}),
Ptr{$relty}, Ptr{BlasInt}),
&m, &n, A, &stride(A,2), jpvt, tau, work, &lwork, rwork, info)
else
ccall(($(string(geqp3)),liblapack), Void,
Expand Down Expand Up @@ -556,11 +555,11 @@ for (gelsd, elty, relty) in ((:zgelsd_, Complex128, Float64),
end

# (GE) general matrices eigenvalue-eigenvector and singular value decompositions
for (geev, gesvd, gesdd, elty) in
((:dgeev_,:dgesvd_,:dgesdd_,:Float64),
(:sgeev_,:sgesvd_,:sgesdd_,:Float32),
(:zgeev_,:zgesvd_,:zgesdd_,:Complex128),
(:cgeev_,:cgesvd_,:cgesdd_,:Complex64))
for (geev, gesvd, gesdd, elty, relty) in
((:dgeev_,:dgesvd_,:dgesdd_,:Float64, Float64),
(:sgeev_,:sgesvd_,:sgesdd_,:Float32, Float32),
(:zgeev_,:zgesvd_,:zgesdd_,:Complex128, Float64),
(:cgeev_,:cgesvd_,:cgesdd_,:Complex64, Float32))
@eval begin
# SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
# $ LDVR, WORK, LWORK, INFO )
Expand All @@ -578,11 +577,10 @@ for (geev, gesvd, gesdd, elty) in
rvecs = jobvr == 'V'
VL = Array($elty, (n, lvecs ? n : 0))
VR = Array($elty, (n, rvecs ? n : 0))
Rtyp = typeof(real(A[1]))
cmplx = iscomplex(A)
if cmplx
W = Array($elty, n)
rwork = Array(Rtyp, 2n)
rwork = Array($relty, 2n)
else
WR = Array($elty, n)
WI = Array($elty, n)
Expand All @@ -596,7 +594,7 @@ for (geev, gesvd, gesdd, elty) in
(Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{Rtyp}, Ptr{BlasInt}),
Ptr{$relty}, Ptr{BlasInt}),
&jobvl, &jobvr, &n, A, &stride(A,2), W, VL, &n, VR, &n,
work, &lwork, rwork, info)
else
Expand Down Expand Up @@ -645,11 +643,10 @@ for (geev, gesvd, gesdd, elty) in
end
work = Array($elty, 1)
lwork = blas_int(-1)
Rtyp = typeof(real(A[1]))
S = Array(Rtyp, minmn)
S = Array($relty, minmn)
cmplx = iscomplex(A)
if cmplx
rwork = Array(Rtyp, job == 'N' ? 7*minmn : 5*minmn*minmn + 5*minmn)
rwork = Array($relty, job == 'N' ? 7*minmn : 5*minmn*minmn + 5*minmn)
end
iwork = Array(BlasInt, 8*minmn)
info = Array(BlasInt, 1)
Expand All @@ -659,7 +656,7 @@ for (geev, gesvd, gesdd, elty) in
(Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{Rtyp}, Ptr{BlasInt}, Ptr{BlasInt}),
Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}),
&job, &m, &n, A, &stride(A,2), S, U, &max(1,stride(U,2)), VT, &max(1,stride(VT,2)),
work, &lwork, rwork, iwork, info)
else
Expand Down Expand Up @@ -693,15 +690,14 @@ for (geev, gesvd, gesdd, elty) in
# $ VT( LDVT, * ), WORK( * )
function gesvd!(jobu::BlasChar, jobvt::BlasChar, A::StridedMatrix{$elty})
chkstride1(A)
Rtyp = typeof(real(A[1]))
m, n = size(A)
minmn = min(m, n)
S = Array(Rtyp, minmn)
S = Array($relty, minmn)
U = Array($elty, jobu == 'A'? (m, m):(jobu == 'S'? (m, minmn) : (m, 0)))
VT = Array($elty, jobvt == 'A'? (n, n):(jobvt == 'S'? (minmn, n) : (n, 0)))
work = Array($elty, 1)
cmplx = iscomplex(A)
if cmplx rwork = Array(Rtyp, 5minmn) end
if cmplx; rwork = Array($relty, 5minmn); end
lwork = blas_int(-1)
info = Array(BlasInt, 1)
for i in 1:2
Expand All @@ -710,7 +706,7 @@ for (geev, gesvd, gesdd, elty) in
(Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{Rtyp}, Ptr{BlasInt}),
Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}),
&jobu, &jobvt, &m, &n, A, &stride(A,2), S, U, &max(1,stride(U,2)), VT, &max(1,stride(VT,2)),
work, &lwork, rwork, info)
else
Expand Down Expand Up @@ -1239,25 +1235,24 @@ for (syconv, syev, sysv, sytrf, sytri, sytrs, elty) in
chkstride1(A)
chksquare(A)
cmplx = iscomplex(A)
Rtyp = typeof(real(A[1]))
n = size(A, 1)
W = Array(Rtyp, n)
W = Array($relty, n)

This comment has been minimized.

Sign in to view

Copy link
@kmsquire

kmsquire Jan 17, 2013

Member

relty doesn't seem to be defined here.

This comment has been minimized.

Sign in to view

Copy link
@ViralBShah

ViralBShah Jan 17, 2013

Author Member

Good catch. Adding it.
work = Array($elty, 1)
lwork = blas_int(-1)
if cmplx
rwork = Array(Rtyp, max(1, 3n-2))
rwork = Array($relty, max(1, 3n-2))
end
info = Array(BlasInt, 1)
for i in 1:2
if cmplx
ccall(($(string(syev)),liblapack), Void,
(Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{Rtyp}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Rtyp}, Ptr{BlasInt}),
Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}),
&jobz, &uplo, &n, A, &stride(A,2), W, work, &lwork, rwork, info)
else
ccall(($(string(syev)),liblapack), Void,
(Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{Rtyp}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
&jobz, &uplo, &n, A, &stride(A,2), W, work, &lwork, info)
end
if info[1] != 0 throw(LapackException(info[1])) end
Expand Down

0 comments on commit d3a8dd1

Please sign in to comment.