← Back to branch summary

~ubuntu-branches/ubuntu/wily/julia/wily

« back to all changes in this revision

Viewing changes to base/lapack.jl

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-02-06 17:54:29 UTC
  • mfrom: (1.1.3)
  • Revision ID: package-import@ubuntu.com-20130206175429-13br5kqpkfjqdmre
Tags: 0.0.0+20130206.git32ff5759-1
* New upstream snapshot.
* debian/copyright: reflect upstream changes
* debian/rules: update get-orig-source to reflect upstream changes
   + Don't ship nginx
   + Adapt for new configure-random target in deps/Makefile
* Enable build of Tk wrapper.
   + debian/control: add build dependency on tk-dev
   + debian/rules: add tk rule to build-arch
* debian/julia.install: install VERSION and COMMIT files
* no-webrepl.patch: new patch
* Refresh other patches
* Add source override for config.status file under deps/random/

Show diffs side-by-side

added added

removed removed

Lines of Context:
base/lapack.jl
143
143
# gegp3 - pivoted QR decomposition
144
144
# gerqf - unpivoted RQ decomposition
145
145
# getrf - LU decomposition
146
 
for (gebrd, gelqf, geqlf, geqrf, geqp3, gerqf, getrf, elty) in
147
 
    ((:dgebrd_,:dgelqf_,:dgeqlf_,:dgeqrf_,:dgeqp3_,:dgerqf_,:dgetrf_,:Float64),
148
 
     (:sgebrd_,:sgelqf_,:sgeqlf_,:sgeqrf_,:sgeqp3_,:sgerqf_,:sgetrf_,:Float32),
149
 
     (:zgebrd_,:zgelqf_,:zgeqlf_,:zgeqrf_,:zgeqp3_,:zgerqf_,:zgetrf_,:Complex128),
150
 
     (:cgebrd_,:cgelqf_,:cgeqlf_,:cgeqrf_,:cgeqp3_,:cgerqf_,:cgetrf_,:Complex64))
 
146
for (gebrd, gelqf, geqlf, geqrf, geqp3, gerqf, getrf, elty, relty) in
 
147
    ((:dgebrd_,:dgelqf_,:dgeqlf_,:dgeqrf_,:dgeqp3_,:dgerqf_,:dgetrf_,:Float64,:Float64),
 
148
     (:sgebrd_,:sgelqf_,:sgeqlf_,:sgeqrf_,:sgeqp3_,:sgerqf_,:sgetrf_,:Float32,:Float32),
 
149
     (:zgebrd_,:zgelqf_,:zgeqlf_,:zgeqrf_,:zgeqp3_,:zgerqf_,:zgetrf_,:Complex128,:Float64),
 
150
     (:cgebrd_,:cgelqf_,:cgeqlf_,:cgeqrf_,:cgeqp3_,:cgerqf_,:cgetrf_,:Complex64,:Float32))
151
151
    @eval begin
152
152
        # SUBROUTINE DGEBRD( M, N, A, LDA, D, E, TAUQ, TAUP, WORK, LWORK,
153
153
        #                    INFO )
248
248
            work  = Array($elty, 1)
249
249
            lwork = blas_int(-1)
250
250
            info  = Array(BlasInt, 1)
251
 
            Rtyp  = typeof(real(A[1]))
252
251
            cmplx = iscomplex(A)
253
 
            if cmplx rwork = Array(Rtyp, 2n) end
 
252
            if cmplx; rwork = Array($relty, 2n); end
254
253
            for i in 1:2
255
254
                if cmplx
256
255
                    ccall(($(string(geqp3)),liblapack), Void,
257
256
                          (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
258
257
                           Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
259
 
                           Ptr{Rtyp}, Ptr{BlasInt}),
 
258
                           Ptr{$relty}, Ptr{BlasInt}),
260
259
                          &m, &n, A, &stride(A,2), jpvt, tau, work, &lwork, rwork, info)
261
260
                else
262
261
                    ccall(($(string(geqp3)),liblapack), Void,
435
434
        function getri!(A::StridedMatrix{$elty}, ipiv::Vector{BlasInt})
436
435
            chkstride1(A)
437
436
            m, n    = size(A)
438
 
            if m != n || n != numel(ipiv) error("getri!: dimension mismatch") end
 
437
            if m != n || n != length(ipiv) error("getri!: dimension mismatch") end
439
438
            lda     = stride(A, 2)
440
439
            info    = Array(BlasInt, 1)
441
440
            lwork   = -1
455
454
        end
456
455
    end
457
456
end
458
 
for (gelsd, elty) in ((:dgelsd_, Float64),
459
 
                      (:sgelsd_, Float32))
 
457
for (gelsd, elty) in ((:dgelsd_, :Float64),
 
458
                      (:sgelsd_, :Float32))
460
459
    @eval begin
461
460
        # SUBROUTINE DGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
462
461
        #      $                   WORK, LWORK, IWORK, INFO )
503
502
        gelsd!(A::StridedMatrix{$elty}, B::StridedVecOrMat{$elty}) = gelsd!(A, B, -1.)
504
503
    end
505
504
end
506
 
for (gelsd, elty, relty) in ((:zgelsd_, Complex128, Float64),
507
 
                             (:cgelsd_, Complex64, Float32))
 
505
for (gelsd, elty, relty) in ((:zgelsd_, :Complex128, :Float64),
 
506
                             (:cgelsd_, :Complex64, :Float32))
508
507
    @eval begin
509
508
        # SUBROUTINE ZGELSD( M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK,
510
509
        #      $                   WORK, LWORK, RWORK, IWORK, INFO )
556
555
end
557
556
 
558
557
# (GE) general matrices eigenvalue-eigenvector and singular value decompositions
559
 
for (geev, gesvd, gesdd, elty) in
560
 
    ((:dgeev_,:dgesvd_,:dgesdd_,:Float64),
561
 
     (:sgeev_,:sgesvd_,:sgesdd_,:Float32),
562
 
     (:zgeev_,:zgesvd_,:zgesdd_,:Complex128),
563
 
     (:cgeev_,:cgesvd_,:cgesdd_,:Complex64))
 
558
for (geev, gesvd, gesdd, elty, relty) in
 
559
    ((:dgeev_,:dgesvd_,:dgesdd_,:Float64,Float64),
 
560
     (:sgeev_,:sgesvd_,:sgesdd_,:Float32,:Float32),
 
561
     (:zgeev_,:zgesvd_,:zgesdd_,:Complex128,:Float64),
 
562
     (:cgeev_,:cgesvd_,:cgesdd_,:Complex64,:Float32))
564
563
    @eval begin
565
564
        #      SUBROUTINE DGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR,
566
565
        #      $                  LDVR, WORK, LWORK, INFO )
578
577
            rvecs = jobvr == 'V'
579
578
            VL    = Array($elty, (n, lvecs ? n : 0))
580
579
            VR    = Array($elty, (n, rvecs ? n : 0))
581
 
            Rtyp  = typeof(real(A[1]))
582
580
            cmplx = iscomplex(A)
583
581
            if cmplx
584
582
                W     = Array($elty, n)
585
 
                rwork = Array(Rtyp, 2n)
 
583
                rwork = Array($relty, 2n)
586
584
            else
587
585
                WR    = Array($elty, n)
588
586
                WI    = Array($elty, n)
596
594
                          (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty},
597
595
                           Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, 
598
596
                           Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
599
 
                           Ptr{Rtyp}, Ptr{BlasInt}),
 
597
                           Ptr{$relty}, Ptr{BlasInt}),
600
598
                          &jobvl, &jobvr, &n, A, &stride(A,2), W, VL, &n, VR, &n,
601
599
                          work, &lwork, rwork, info)
602
600
                else
645
643
            end
646
644
            work   = Array($elty, 1)
647
645
            lwork  = blas_int(-1)
648
 
            Rtyp   = typeof(real(A[1]))
649
 
            S      = Array(Rtyp, minmn)
 
646
            S      = Array($relty, minmn)
650
647
            cmplx  = iscomplex(A)
651
648
            if cmplx
652
 
                rwork = Array(Rtyp, job == 'N' ? 5*minmn : minmn*max(5*minmn+7,2*max(m,n)+2*minmn+1))
 
649
                rwork = Array($relty, job == 'N' ? 7*minmn : 5*minmn*minmn + 5*minmn)
653
650
            end
654
651
            iwork  = Array(BlasInt, 8*minmn)
655
652
            info   = Array(BlasInt, 1)
659
656
                          (Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
660
657
                           Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
661
658
                           Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
662
 
                           Ptr{Rtyp}, Ptr{BlasInt}, Ptr{BlasInt}),
 
659
                           Ptr{$relty}, Ptr{BlasInt}, Ptr{BlasInt}),
663
660
                          &job, &m, &n, A, &stride(A,2), S, U, &max(1,stride(U,2)), VT, &max(1,stride(VT,2)),
664
661
                          work, &lwork, rwork, iwork, info)
665
662
                else
693
690
        #      $                   VT( LDVT, * ), WORK( * )
694
691
        function gesvd!(jobu::BlasChar, jobvt::BlasChar, A::StridedMatrix{$elty})
695
692
            chkstride1(A)
696
 
            Rtyp   = typeof(real(A[1]))
697
693
            m, n   = size(A)
698
694
            minmn  = min(m, n)
699
 
            S      = Array(Rtyp, minmn)
 
695
            S      = Array($relty, minmn)
700
696
            U      = Array($elty, jobu  == 'A'? (m, m):(jobu  == 'S'? (m, minmn) : (m, 0)))
701
697
            VT     = Array($elty, jobvt == 'A'? (n, n):(jobvt == 'S'? (minmn, n) : (n, 0)))
702
698
            work   = Array($elty, 1)
703
699
            cmplx  = iscomplex(A)
704
 
            if cmplx rwork = Array(Rtyp, 5minmn) end
 
700
            if cmplx; rwork = Array($relty, 5minmn); end
705
701
            lwork  = blas_int(-1)
706
702
            info   = Array(BlasInt, 1)
707
703
            for i in 1:2
710
706
                          (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{BlasInt},
711
707
                           Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty},
712
708
                           Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
713
 
                           Ptr{BlasInt}, Ptr{Rtyp}, Ptr{BlasInt}),
 
709
                           Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}),
714
710
                          &jobu, &jobvt, &m, &n, A, &stride(A,2), S, U, &max(1,stride(U,2)), VT, &max(1,stride(VT,2)),
715
711
                          work, &lwork, rwork, info)
716
712
                else
1130
1126
## (TR) triangular matrices: solver and inverse
1131
1127
for (trtri, trtrs, elty) in
1132
1128
    ((:dtrtri_,:dtrtrs_,:Float64),
1133
 
     (:strtri_,:strtrs_,Float32),
 
1129
     (:strtri_,:strtrs_,:Float32),
1134
1130
     (:ztrtri_,:ztrtrs_,:Complex128),
1135
1131
     (:ctrtri_,:ctrtrs_,:Complex64))
1136
1132
    @eval begin
1203
1199
 
1204
1200
## (SY) symmetric matrices - eigendecomposition, Bunch-Kaufman decomposition,
1205
1201
## solvers (direct and factored) and inverse.
1206
 
for (syconv, syev, sysv, sytrf, sytri, sytrs, elty) in
1207
 
    ((:dsyconv_,:dsyev_,:dsysv_,:dsytrf_,:dsytri_,:dsytrs_,:Float64),
1208
 
     (:ssyconv_,:ssyev_,:ssysv_,:ssytrf_,:ssytri_,:ssytrs_,:Float32),
1209
 
     (:zheconv_,:zheev_,:zhesv_,:zhetrf_,:zhetri_,:zhetrs_,:Complex128),
1210
 
     (:checonv_,:cheev_,:chesv_,:chetrf_,:chetri_,:chetrs_,:Complex64))
 
1202
for (syconv, syev, sysv, sytrf, sytri, sytrs, elty, relty) in
 
1203
    ((:dsyconv_,:dsyev_,:dsysv_,:dsytrf_,:dsytri_,:dsytrs_,:Float64, :Float64),
 
1204
     (:ssyconv_,:ssyev_,:ssysv_,:ssytrf_,:ssytri_,:ssytrs_,:Float32, :Float32),
 
1205
     (:zheconv_,:zheev_,:zhesv_,:zhetrf_,:zhetri_,:zhetrs_,:Complex128, :Float64),
 
1206
     (:checonv_,:cheev_,:chesv_,:chetrf_,:chetri_,:chetrs_,:Complex64, :Float32))
1211
1207
    @eval begin
1212
1208
        #       SUBROUTINE DSYCONV( UPLO, WAY, N, A, LDA, IPIV, WORK, INFO )
1213
1209
        # *     .. Scalar Arguments ..
1239
1235
            chkstride1(A)
1240
1236
            chksquare(A)
1241
1237
            cmplx = iscomplex(A)
1242
 
            Rtyp  = typeof(real(A[1]))
1243
1238
            n     = size(A, 1)
1244
 
            W     = Array(Rtyp, n)
 
1239
            W     = Array($relty, n)
1245
1240
            work  = Array($elty, 1)
1246
1241
            lwork = blas_int(-1)
1247
1242
            if cmplx
1248
 
                rwork = Array(Rtyp, max(1, 3n-2))
 
1243
                rwork = Array($relty, max(1, 3n-2))
1249
1244
            end
1250
1245
            info  = Array(BlasInt, 1)
1251
1246
            for i in 1:2
1252
1247
                if cmplx
1253
1248
                    ccall(($(string(syev)),liblapack), Void,
1254
1249
                          (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
1255
 
                          Ptr{Rtyp}, Ptr{$elty}, Ptr{BlasInt}, Ptr{Rtyp}, Ptr{BlasInt}),
 
1250
                          Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$relty}, Ptr{BlasInt}),
1256
1251
                          &jobz, &uplo, &n, A, &stride(A,2), W, work, &lwork, rwork, info)
1257
1252
                else
1258
1253
                    ccall(($(string(syev)),liblapack), Void,
1259
1254
                          (Ptr{Uint8}, Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
1260
 
                          Ptr{Rtyp}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
 
1255
                          Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}),
1261
1256
                          &jobz, &uplo, &n, A, &stride(A,2), W, work, &lwork, info)
1262
1257
                end
1263
1258
                if info[1] != 0 throw(LapackException(info[1])) end
1398
1393
        end
1399
1394
    end
1400
1395
end
 
1396
 
 
1397
# New symmetric eigen solver
1401
1398
for (syevr, elty) in
1402
1399
    ((:dsyevr_,:Float64),
1403
1400
     (:ssyevr_,:Float32))
1521
1518
                    lwork = blas_int(real(work[1]))
1522
1519
                    work = Array($elty, lwork)
1523
1520
                    lrwork = blas_int(rwork[1])
1524
 
                    rwork = Array($elty, lrwork)
 
1521
                    rwork = Array($relty, lrwork)
1525
1522
                    liwork = iwork[1]
1526
1523
                    iwork = Array(BlasInt, liwork)
1527
1524
                end
1557
1554
            iwork = Array(BlasInt, n)
1558
1555
            info = Array(BlasInt, 1)
1559
1556
            ccall(($(string(gecon)),liblapack), Void,
1560
 
                (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, 
1561
 
                    Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
1562
 
                    Ptr{BlasInt}),
1563
 
                &normtype, &n, A, &lda,
1564
 
                &anorm, rcond, work, iwork,
1565
 
                info)
 
1557
                  (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, 
 
1558
                   Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
 
1559
                   Ptr{BlasInt}),
 
1560
                  &normtype, &n, A, &lda, &anorm, rcond, work, iwork,
 
1561
                  info)
1566
1562
            if info[1] != 0 throw(LapackException(info[1])) end
1567
1563
            return rcond[1]
1568
1564
        end
1570
1566
end
1571
1567
for (gecon, elty, relty) in
1572
1568
    ((:zgecon_,:Complex128,:Float64),
1573
 
     (:cgecon_,:Complex64,:Float32))
 
1569
     (:cgecon_,:Complex64, :Float32))
1574
1570
    @eval begin
1575
1571
        function gecon!(normtype::BlasChar, A::StridedMatrix{$elty}, anorm::$relty)
1576
1572
            chkstride1(A)
1589
1585
            lda = max(1, size(A, 1))
1590
1586
            rcond = Array($relty, 1)
1591
1587
            work = Array($elty, 2n)
1592
 
            rwork = Array(BlasInt, 2n)
 
1588
            rwork = Array($relty, 2n)
1593
1589
            info = Array(BlasInt, 1)
1594
1590
            ccall(($(string(gecon)),liblapack), Void,
1595
 
                (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, 
1596
 
                    Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty},
1597
 
                    Ptr{BlasInt}),
1598
 
                &normtype, &n, A, &lda, 
1599
 
                &anorm, rcond, work, rwork,
1600
 
                info)
 
1591
                  (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, 
 
1592
                   Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty},
 
1593
                   Ptr{BlasInt}),
 
1594
                  &normtype, &n, A, &lda, &anorm, rcond, work, rwork,
 
1595
                  info)
1601
1596
            if info[1] < 0 throw(LapackException(info[1])) end
1602
1597
            return rcond[1]
1603
1598
        end
1604
1599
    end
1605
1600
end
 
1601
 
 
1602
# Hessenberg form
 
1603
for (gehrd, elty) in
 
1604
    ((:dgehrd_,:Float64),
 
1605
     (:sgehrd_,:Float32),
 
1606
     (:zgehrd_,:Complex128),
 
1607
     (:cgehrd_,:Complex64))
 
1608
    @eval begin
 
1609
        function gehrd!(ilo::BlasInt, ihi::BlasInt, A::StridedMatrix{$elty})
 
1610
#                 .. Scalar Arguments ..
 
1611
#       INTEGER            IHI, ILO, INFO, LDA, LWORK, N
 
1612
# *     ..
 
1613
# *     .. Array Arguments ..
 
1614
#       DOUBLE PRECISION  A( LDA, * ), TAU( * ), WORK( * )
 
1615
            chkstride1(A)
 
1616
            chksquare(A)
 
1617
            n = size(A, 1)
 
1618
            tau = Array($elty, n - 1)
 
1619
            work = Array($elty, 1)
 
1620
            lwork = blas_int(-1)
 
1621
            info = Array(BlasInt, 1)
 
1622
            for i = 1:2
 
1623
                ccall(($(string(gehrd)),liblapack), Void,
 
1624
                    (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
 
1625
                     Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
 
1626
                     Ptr{BlasInt}),
 
1627
                    &n, &ilo, &ihi, A, 
 
1628
                    &max(1,n), tau, work, &lwork, 
 
1629
                    info)
 
1630
                if info[1] < 0 throw(LapackException(info[1])) end
 
1631
                if lwork < 0
 
1632
                    lwork = blas_int(work[1])
 
1633
                    work = Array($elty, lwork)
 
1634
                end
 
1635
            end
 
1636
            return A, tau
 
1637
        end
 
1638
    end
 
1639
end
 
1640
gehrd!(A::StridedMatrix) = gehrd!(blas_int(1), blas_int(size(A, 1)), A)
 
1641
 
 
1642
# construct Q from Hessenberg
 
1643
for (orghr, elty) in
 
1644
    ((:dorghr_,:Float64),
 
1645
     (:sorghr_,:Float32),
 
1646
     (:zunghr_,:Complex128),
 
1647
     (:cunghr_,:Complex64))
 
1648
    @eval begin
 
1649
        function orghr!(ilo::BlasInt, ihi::BlasInt, A::StridedMatrix{$elty}, tau::StridedVector{$elty})
 
1650
# *     .. Scalar Arguments ..
 
1651
#       INTEGER            IHI, ILO, INFO, LDA, LWORK, N
 
1652
# *     ..
 
1653
# *     .. Array Arguments ..
 
1654
#       DOUBLE PRECISION   A( LDA, * ), TAU( * ), WORK( * )
 
1655
            chkstride1(A)
 
1656
            chksquare(A)
 
1657
            n = size(A, 1)
 
1658
            if n - length(tau) != 1 throw(LapackDimMismatch) end
 
1659
            work = Array($elty, 1)
 
1660
            lwork = blas_int(-1)
 
1661
            info = Array(BlasInt, 1)
 
1662
            for i = 1:2
 
1663
                ccall(($(string(orghr)),liblapack), Void,
 
1664
                    (Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
 
1665
                     Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt},
 
1666
                     Ptr{BlasInt}),
 
1667
                    &n, &ilo, &ihi, A, 
 
1668
                    &max(1,n), tau, work, &lwork, 
 
1669
                    info)
 
1670
                if info[1] < 0 throw(LapackException(info[1])) end
 
1671
                if lwork < 0
 
1672
                    lwork = blas_int(work[1])
 
1673
                    work = Array($elty, lwork)
 
1674
                end
 
1675
            end
 
1676
            return A
 
1677
        end
 
1678
    end
 
1679
end
 
1680
# Schur form
 
1681
for (gees, elty) in
 
1682
    ((:dgees_,:Float64),
 
1683
     (:sgees_,:Float32))
 
1684
    @eval begin
 
1685
        function gees!(jobvs::BlasChar, A::StridedMatrix{$elty})
 
1686
#     .. Scalar Arguments ..
 
1687
#     CHARACTER          JOBVS, SORT
 
1688
#     INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
 
1689
#     ..
 
1690
#     .. Array Arguments ..
 
1691
#     LOGICAL            BWORK( * )
 
1692
#     DOUBLE PRECISION   A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
 
1693
#    $                   WR( * )
 
1694
            chkstride1(A)
 
1695
            chksquare(A)
 
1696
            sort = 'N'
 
1697
            n = size(A, 1)
 
1698
            sdim = Array(BlasInt, 1)
 
1699
            wr = Array($elty, n)
 
1700
            wi = Array($elty, n)
 
1701
            ldvs = jobvs == 'V' ? n : 1
 
1702
            vs = Array($elty, ldvs, n)
 
1703
            work = Array($elty, 1)
 
1704
            lwork = blas_int(-1)
 
1705
            info = Array(BlasInt, 1)
 
1706
            for i = 1:2
 
1707
                ccall(($(string(gees)),liblapack), Void,
 
1708
                    (Ptr{Uint8}, Ptr{Uint8}, Ptr{Void}, Ptr{BlasInt},
 
1709
                        Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
 
1710
                        Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
 
1711
                        Ptr{BlasInt}, Ptr{Void}, Ptr{BlasInt}),
 
1712
                    &jobvs, &sort, [], &n, 
 
1713
                        A, &max(1, n), sdim, wr,
 
1714
                        wi, vs, &ldvs, work, 
 
1715
                        &lwork, [], info)
 
1716
                if info[1] != 0 throw(LapackException(info[1])) end
 
1717
                if lwork < 0
 
1718
                    lwork = blas_int(work[1])
 
1719
                    work = Array($elty, lwork)
 
1720
                end
 
1721
            end
 
1722
            if all(wi .== 0)
 
1723
                return A, vs, wr
 
1724
            else
 
1725
                return A, vs, complex(wr, wi)
 
1726
            end
 
1727
        end
 
1728
    end
 
1729
end
 
1730
for (gees, elty, relty) in
 
1731
    ((:zgees_,:Complex128,:Float64),
 
1732
     (:cgees_,:Complex64,:Float32))
 
1733
    @eval begin
 
1734
        function gees!(jobvs::BlasChar, A::StridedMatrix{$elty})
 
1735
# *     .. Scalar Arguments ..
 
1736
#       CHARACTER          JOBVS, SORT
 
1737
#       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
 
1738
# *     ..
 
1739
# *     .. Array Arguments ..
 
1740
#       LOGICAL            BWORK( * )
 
1741
#       DOUBLE PRECISION   RWORK( * )
 
1742
#       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
 
1743
            chkstride1(A)
 
1744
            chksquare(A)
 
1745
            sort = 'N'
 
1746
            n = size(A, 1)
 
1747
            sdim = Array(BlasInt, 1)
 
1748
            w = Array($elty, n)
 
1749
            ldvs = jobvs == 'V' ? n : 1
 
1750
            vs = Array($elty, ldvs, n)
 
1751
            work = Array($elty, 1)
 
1752
            lwork = blas_int(-1)
 
1753
            rwork = Array($relty, n)
 
1754
            info = Array(BlasInt, 1)
 
1755
            for i = 1:2
 
1756
                ccall(($(string(gees)),liblapack), Void,
 
1757
                    (Ptr{Uint8}, Ptr{Uint8}, Ptr{Void}, Ptr{BlasInt},
 
1758
                        Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
 
1759
                        Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, 
 
1760
                        Ptr{$relty}, Ptr{Void}, Ptr{BlasInt}),
 
1761
                    &jobvs, &sort, [], &n, 
 
1762
                        A, &max(1, n), sdim, w,
 
1763
                        vs, &ldvs, work, &lwork, 
 
1764
                        rwork, [], info)
 
1765
                if info[1] != 0 throw(LapackException(info[1])) end
 
1766
                if lwork < 0
 
1767
                    lwork = blas_int(work[1])
 
1768
                    work = Array($elty, lwork)
 
1769
                end
 
1770
            end
 
1771
            return A, vs, w
 
1772
        end
 
1773
    end
 
1774
end
1606
1775
end # module