← Back to branch summary

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

« back to all changes in this revision

Viewing changes to base/linalg_dense.jl

  • Committer: Package Import Robot
  • Author(s): Sébastien Villemot
  • Date: 2013-02-11 03:51:26 UTC
  • mfrom: (1.1.4)
  • Revision ID: package-import@ubuntu.com-20130211035126-hap464pbhd97wjbl
Tags: 0.1~20130211.git86fbe98d-1
* New upstream snapshot.
* debian/control:
   + add git to Recommends (for Julia package manager)
   + remove dependencies on libglpk-dev (it moved to its own package)
   + add explicit dependency on libgmp10 (there is no more a wrapper)
* fix-clean-rules.patch: remove patch, applied upstream
* gsvddense_blasint.patch: new patch
* Refresh other patches

Show diffs side-by-side

added added

removed removed

Lines of Context:
base/linalg_dense.jl
1
1
# Linear algebra functions for dense matrices in column major format
2
2
 
3
3
scale!(X::Array{Float32}, s::Real) = BLAS.scal!(length(X), float32(s), X, 1)
4
 
 
5
4
scale!(X::Array{Float64}, s::Real) = BLAS.scal!(length(X), float64(s), X, 1)
6
 
 
7
5
scale!(X::Array{Complex64}, s::Real) = (ccall(("sscal_",Base.libblas_name), Void, (Ptr{BlasInt}, Ptr{Float32}, Ptr{Complex64}, Ptr{BlasInt}), &(2*length(X)), &s, X, &1); X)
8
 
 
9
6
scale!(X::Array{Complex128}, s::Real) = (ccall(("dscal_",Base.libblas_name), Void, (Ptr{BlasInt}, Ptr{Float64}, Ptr{Complex128}, Ptr{BlasInt}), &(2*length(X)), &s, X, &1); X)
10
7
 
11
8
#Test whether a matrix is positive-definite
389
386
 
390
387
\{T<:BlasFloat}(B::BunchKaufman{T}, R::StridedVecOrMat{T}) =
391
388
    LAPACK.sytrs!(B.UL, B.LD, B.ipiv, copy(R))
392
 
    
 
389
 
393
390
type CholeskyDense{T<:BlasFloat} <: Factorization{T}
394
391
    LR::Matrix{T}
395
392
    UL::BlasChar
396
393
    function CholeskyDense(A::Matrix{T}, UL::BlasChar)
397
394
        A, info = LAPACK.potrf!(UL, A)
398
 
        if info != 0 error("Matrix A not positive-definite") end
 
395
        if info != 0; throw(LAPACK.PosDefException(info)); end
399
396
        if UL == 'U'
400
397
            new(triu!(A), UL)
401
398
        elseif UL == 'L'
423
420
    
424
421
function inv(C::CholeskyDense)
425
422
    Ci, info = LAPACK.potri!(C.UL, copy(C.LR))
426
 
    if info != 0 error("Matrix singular") end 
 
423
    if info != 0; throw(LAPACK.SingularException(info)); end 
427
424
    symmetrize!(Ci, C.UL)
428
425
end
429
426
 
430
427
## Should these functions check that the matrix is Hermitian?
431
 
chold!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = CholeskyDense{T}(A, UL)
432
 
chold!{T<:BlasFloat}(A::Matrix{T}) = chold!(A, 'U')
433
 
chold{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = chold!(copy(A), UL)
434
 
chold{T<:Number}(A::Matrix{T}, UL::BlasChar) = chold(float64(A), UL)
435
 
chold{T<:Number}(A::Matrix{T}) = chold(A, 'U')
 
428
cholfact!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = CholeskyDense{T}(A, UL)
 
429
cholfact!{T<:BlasFloat}(A::Matrix{T}) = cholfact!(A, 'U')
 
430
cholfact{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholfact!(copy(A), UL)
 
431
cholfact{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholfact(float64(A), UL)
 
432
cholfact{T<:Number}(A::Matrix{T}) = cholfact(A, 'U')
436
433
 
437
434
## Matlab (and R) compatible
438
 
chol(A::Matrix, UL::BlasChar) = factors(chold(A, UL))
 
435
chol(A::Matrix, UL::BlasChar) = factors(cholfact(A, UL))
439
436
chol(A::Matrix) = chol(A, 'U')
440
437
chol(x::Number) = imag(x) == 0 && real(x) > 0 ? sqrt(x) : error("Argument not positive-definite")
441
438
 
442
 
type CholeskyDensePivoted{T<:BlasFloat} <: Factorization{T}
 
439
type CholeskyPivotedDense{T<:BlasFloat} <: Factorization{T}
443
440
    LR::Matrix{T}
444
441
    UL::BlasChar
445
442
    piv::Vector{BlasInt}
446
443
    rank::BlasInt
447
444
    tol::Real
448
 
    function CholeskyDensePivoted(A::Matrix{T}, UL::BlasChar, tol::Real)
 
445
    function CholeskyPivotedDense(A::Matrix{T}, UL::BlasChar, tol::Real)
449
446
        A, piv, rank, info = LAPACK.pstrf!(UL, A, tol)
 
447
        if info != 0; throw(LAPACK.RankDeficientException(info)); end
450
448
        if UL == 'U'
451
449
            new(triu!(A), UL, piv, rank, tol)
452
450
        elseif UL == 'L'
457
455
    end
458
456
end
459
457
 
460
 
size(C::CholeskyDensePivoted) = size(C.LR)
461
 
size(C::CholeskyDensePivoted,d::Integer) = size(C.LR,d)
462
 
 
463
 
factors(C::CholeskyDensePivoted) = C.LR, C.piv
464
 
 
465
 
function \{T<:BlasFloat}(C::CholeskyDensePivoted{T}, B::StridedVector{T})
466
 
    if C.rank < size(C.LR, 1) error("Matrix is not positive-definite") end
 
458
size(C::CholeskyPivotedDense) = size(C.LR)
 
459
size(C::CholeskyPivotedDense,d::Integer) = size(C.LR,d)
 
460
 
 
461
factors(C::CholeskyPivotedDense) = C.LR, C.piv
 
462
 
 
463
function \{T<:BlasFloat}(C::CholeskyPivotedDense{T}, B::StridedVector{T})
 
464
    if C.rank < size(C.LR, 1); throw(LAPACK.RankDeficientException(info)); end
467
465
    LAPACK.potrs!(C.UL, C.LR, copy(B)[C.piv])[invperm(C.piv)]
468
466
end
469
 
function \{T<:BlasFloat}(C::CholeskyDensePivoted{T}, B::StridedMatrix{T})
470
 
    if C.rank < size(C.LR, 1) error("Matrix is not positive-definite") end
 
467
 
 
468
function \{T<:BlasFloat}(C::CholeskyPivotedDense{T}, B::StridedMatrix{T})
 
469
    if C.rank < size(C.LR, 1); throw(LAPACK.RankDeficientException(info)); end
471
470
    LAPACK.potrs!(C.UL, C.LR, copy(B)[C.piv,:])[invperm(C.piv),:]
472
471
end
473
472
 
474
 
rank(C::CholeskyDensePivoted) = C.rank
 
473
rank(C::CholeskyPivotedDense) = C.rank
475
474
 
476
 
function det{T}(C::CholeskyDensePivoted{T})
 
475
function det{T}(C::CholeskyPivotedDense{T})
477
476
    if C.rank < size(C.LR, 1) 
478
477
        return real(zero(T))
479
478
    else 
481
480
    end
482
481
end
483
482
    
484
 
function inv(C::CholeskyDensePivoted)
 
483
function inv(C::CholeskyPivotedDense)
485
484
    if C.rank < size(C.LR, 1) error("Matrix singular") end
486
485
    Ci, info = LAPACK.potri!(C.UL, copy(C.LR))
487
486
    if info != 0 error("Matrix is singular") end
490
489
end
491
490
 
492
491
## Should these functions check that the matrix is Hermitian?
493
 
cholpd!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = CholeskyDensePivoted{T}(A, UL, tol)
494
 
cholpd!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpd!(A, UL, -1.)
495
 
cholpd!{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpd!(A, 'U', tol)
496
 
cholpd!{T<:BlasFloat}(A::Matrix{T}) = cholpd!(A, 'U', -1.)
497
 
cholpd{T<:Number}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpd(float64(A), UL, tol)
498
 
cholpd{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholpd(float64(A), UL, -1.)
499
 
cholpd{T<:Number}(A::Matrix{T}, tol::Real) = cholpd(float64(A), true, tol)
500
 
cholpd{T<:Number}(A::Matrix{T}) = cholpd(float64(A), 'U', -1.)
501
 
cholpd{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpd!(copy(A), UL, tol)
502
 
cholpd{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpd!(copy(A), UL, -1.)
503
 
cholpd{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpd!(copy(A), 'U', tol)
504
 
cholpd{T<:BlasFloat}(A::Matrix{T}) = cholpd!(copy(A), 'U', -1.)
 
492
cholpfact!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = CholeskyPivotedDense{T}(A, UL, tol)
 
493
cholpfact!{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpfact!(A, UL, -1.)
 
494
cholpfact!{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpfact!(A, 'U', tol)
 
495
cholpfact!{T<:BlasFloat}(A::Matrix{T}) = cholpfact!(A, 'U', -1.)
 
496
cholpfact{T<:Number}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpfact(float64(A), UL, tol)
 
497
cholpfact{T<:Number}(A::Matrix{T}, UL::BlasChar) = cholpfact(float64(A), UL, -1.)
 
498
cholpfact{T<:Number}(A::Matrix{T}, tol::Real) = cholpfact(float64(A), true, tol)
 
499
cholpfact{T<:Number}(A::Matrix{T}) = cholpfact(float64(A), 'U', -1.)
 
500
cholpfact{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar, tol::Real) = cholpfact!(copy(A), UL, tol)
 
501
cholpfact{T<:BlasFloat}(A::Matrix{T}, UL::BlasChar) = cholpfact!(copy(A), UL, -1.)
 
502
cholpfact{T<:BlasFloat}(A::Matrix{T}, tol::Real) = cholpfact!(copy(A), 'U', tol)
 
503
cholpfact{T<:BlasFloat}(A::Matrix{T}) = cholpfact!(copy(A), 'U', -1.)
505
504
 
506
505
type LUDense{T} <: Factorization{T}
507
506
    lu::Matrix{T}
531
530
    L, U, P
532
531
end
533
532
 
534
 
function lud!{T<:BlasFloat}(A::Matrix{T})
 
533
function lufact!{T<:BlasFloat}(A::Matrix{T})
535
534
    lu, ipiv, info = LAPACK.getrf!(A)
536
535
    LUDense{T}(lu, ipiv, info)
537
536
end
538
537
 
539
 
lud{T<:BlasFloat}(A::Matrix{T}) = lud!(copy(A))
540
 
lud{T<:Number}(A::Matrix{T}) = lud(float64(A))
 
538
lufact{T<:BlasFloat}(A::Matrix{T}) = lufact!(copy(A))
 
539
lufact{T<:Number}(A::Matrix{T}) = lufact(float64(A))
541
540
 
542
541
## Matlab-compatible
543
 
lu{T<:Number}(A::Matrix{T}) = factors(lud(A))
 
542
lu{T<:Number}(A::Matrix{T}) = factors(lufact(A))
544
543
lu(x::Number) = (one(x), x, [1])
545
544
 
546
545
function det{T}(lu::LUDense{T})
572
571
size(A::QRDense) = size(A.hh)
573
572
size(A::QRDense,n) = size(A.hh,n)
574
573
 
575
 
qrd!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDense{T}(LAPACK.geqrf!(A)...)
576
 
qrd{T<:BlasFloat}(A::StridedMatrix{T}) = qrd!(copy(A))
577
 
qrd{T<:Real}(A::StridedMatrix{T}) = qrd(float64(A))
 
574
qrfact!{T<:BlasFloat}(A::StridedMatrix{T}) = QRDense{T}(LAPACK.geqrf!(A)...)
 
575
qrfact{T<:BlasFloat}(A::StridedMatrix{T}) = qrfact!(copy(A))
 
576
qrfact{T<:Real}(A::StridedMatrix{T}) = qrfact(float64(A))
578
577
 
579
 
function factors{T<:BlasFloat}(qrd::QRDense{T})
580
 
    aa  = copy(qrd.hh)
 
578
function factors{T<:BlasFloat}(qrfact::QRDense{T})
 
579
    aa  = copy(qrfact.hh)
581
580
    R   = triu(aa[1:min(size(aa)),:])   # must be *before* call to orgqr!
582
 
    LAPACK.orgqr!(aa, qrd.tau, size(aa,2)), R
 
581
    LAPACK.orgqr!(aa, qrfact.tau, size(aa,2)), R
583
582
end
584
583
 
585
 
qr{T<:Number}(x::StridedMatrix{T}) = factors(qrd(x))
 
584
qr{T<:Number}(x::StridedMatrix{T}) = factors(qrfact(x))
586
585
qr(x::Number) = (one(x), x)
587
586
 
588
587
## Multiplication by Q from the QR decomposition
589
 
(*){T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
 
588
qmulQR{T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
590
589
    LAPACK.ormqr!('L', 'N', A.hh, size(A.hh,2), A.tau, copy(B))
591
590
 
592
591
## Multiplication by Q' from the QR decomposition
593
 
Ac_mul_B{T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
 
592
qTmulQR{T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T}) =
594
593
    LAPACK.ormqr!('L', iscomplex(A.tau)?'C':'T', A.hh, size(A.hh,2), A.tau, copy(B))
595
594
 
596
595
## Least squares solution.  Should be more careful about cases with m < n
597
596
function (\){T<:BlasFloat}(A::QRDense{T}, B::StridedVecOrMat{T})
598
597
    n   = length(A.tau)
599
 
    ans, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(A'*B)[1:n,:])
 
598
    ans, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(qTmulQR(A,B))[1:n,:])
600
599
    if info > 0; throw(LAPACK.SingularException(info)); end
601
600
    return ans
602
601
end
603
602
 
604
 
type QRPDense{T} <: Factorization{T}
 
603
type QRPivotedDense{T} <: Factorization{T}
605
604
    hh::Matrix{T}
606
605
    tau::Vector{T}
607
606
    jpvt::Vector{BlasInt}
608
 
    function QRPDense(hh::Matrix{T}, tau::Vector{T}, jpvt::Vector{BlasInt})
 
607
    function QRPivotedDense(hh::Matrix{T}, tau::Vector{T}, jpvt::Vector{BlasInt})
609
608
        m, n = size(hh)
610
609
        if length(tau) != min(m,n) || length(jpvt) != n
611
 
            error("QRPDense: mismatched dimensions")
 
610
            error("QRPivotedDense: mismatched dimensions")
612
611
        end
613
612
        new(hh,tau,jpvt)
614
613
    end
615
614
end
616
 
size(x::QRPDense)   = size(x.hh)
617
 
size(x::QRPDense,d) = size(x.hh,d)
 
615
size(x::QRPivotedDense)   = size(x.hh)
 
616
size(x::QRPivotedDense,d) = size(x.hh,d)
618
617
## Multiplication by Q from the QR decomposition
619
 
(*){T<:BlasFloat}(A::QRPDense{T}, B::StridedVecOrMat{T}) =
 
618
qmulQR{T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T}) =
620
619
    LAPACK.ormqr!('L', 'N', A.hh, size(A,2), A.tau, copy(B))
621
620
## Multiplication by Q' from the QR decomposition
622
 
Ac_mul_B{T<:BlasFloat}(A::QRPDense{T}, B::StridedVecOrMat{T}) =
 
621
qTmulQR{T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T}) =
623
622
    LAPACK.ormqr!('L', iscomplex(A.tau)?'C':'T', A.hh, size(A,2), A.tau, copy(B))
624
623
 
625
 
qrpd!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPDense{T}(LAPACK.geqp3!(A)...)
626
 
qrpd{T<:BlasFloat}(A::StridedMatrix{T}) = qrpd!(copy(A))
627
 
qrpd{T<:Real}(x::StridedMatrix{T}) = qrpd(float64(x))
 
624
qrpfact!{T<:BlasFloat}(A::StridedMatrix{T}) = QRPivotedDense{T}(LAPACK.geqp3!(A)...)
 
625
qrpfact{T<:BlasFloat}(A::StridedMatrix{T}) = qrpfact!(copy(A))
 
626
qrpfact{T<:Real}(x::StridedMatrix{T}) = qrpfact(float64(x))
628
627
 
629
 
function factors{T<:BlasFloat}(x::QRPDense{T})
 
628
function factors{T<:BlasFloat}(x::QRPivotedDense{T})
630
629
    aa = copy(x.hh)
631
630
    R  = triu(aa[1:min(size(aa)),:])
632
631
    LAPACK.orgqr!(aa, x.tau, size(aa,2)), R, x.jpvt
633
632
end
634
633
 
635
 
qrp{T<:BlasFloat}(x::StridedMatrix{T}) = factors(qrpd(x))
 
634
qrp{T<:BlasFloat}(x::StridedMatrix{T}) = factors(qrpfact(x))
636
635
qrp{T<:Real}(x::StridedMatrix{T}) = qrp(float64(x))
637
636
 
638
 
function (\){T<:BlasFloat}(A::QRPDense{T}, B::StridedVecOrMat{T})
 
637
function (\){T<:BlasFloat}(A::QRPivotedDense{T}, B::StridedVecOrMat{T})
639
638
    n = length(A.tau)
640
 
    x, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(A'*B)[1:n,:])
 
639
    x, info = LAPACK.trtrs!('U','N','N',A.hh[1:n,:],(qTmulQR(A,B))[1:n,:])
641
640
    if info > 0; throw(LAPACK.SingularException(info)); end
642
641
    isa(B, Vector) ? x[invperm(A.jpvt)] : x[:,invperm(A.jpvt)]
643
642
end
675
674
    m, n = size(A)
676
675
    if m != n; error("det only defined for square matrices"); end
677
676
    if istriu(A) | istril(A); return prod(diag(A)); end
678
 
    return det(lud(A))
 
677
    return det(lufact(A))
679
678
end
680
679
det(x::Number) = x
681
680
 
728
727
eig(x) = eig(x, true)
729
728
eigvals(x) = eig(x, false)
730
729
 
731
 
# This is the svd based on the LAPACK GESVD, which is slower, but takes
732
 
# lesser memory. It should be made available through a keyword argument
733
 
# at a later date.
734
 
#
735
 
# function svd{T<:BlasFloat}(A::StridedMatrix{T},vecs::Bool,thin::Bool)
736
 
#     m,n = size(A)
737
 
#     if m == 0 || n == 0
738
 
#         if vecs; return (eye(m, thin ? n : m), zeros(0), eye(n,n)); end
739
 
#         return (zeros(T, 0, 0), zeros(T, 0), zeros(T, 0, 0))
740
 
#     end
741
 
#     if vecs; return LAPACK.gesvd!(thin ? 'S' : 'A', thin ? 'S' : 'A', copy(A)); end
742
 
#     LAPACK.gesvd!('N', 'N', copy(A))
743
 
# end
744
 
#
745
 
# svd{T<:Integer}(x::StridedMatrix{T},vecs,thin) = svd(float64(x),vecs,thin)
746
 
# svd(A::StridedMatrix) = svd(A,true,false)
747
 
# svd(A::StridedMatrix, thin::Bool) = svd(A,true,thin)
748
 
# svdvals(A) = svd(A,false,true)[2]
749
 
 
750
 
svdt(x::Number,vecs::Bool,thin::Bool) = vecs ? (x==0?one(x):x/abs(x),abs(x),one(x)) : ([],abs(x),[])
751
 
svdt(x::Number,vecs::Bool,thin::Bool) = svdt(x, vecs, thin)
752
 
 
753
 
function svdt{T<:BlasFloat}(A::StridedMatrix{T},vecs::Bool,thin::Bool)
754
 
    m,n = size(A)
755
 
    if m == 0 || n == 0
756
 
        if vecs; return (eye(m, thin ? n : m), zeros(0), eye(n,n)); end
 
730
# SVD
 
731
type SVDDense{T,Tr} <: Factorization{T}
 
732
    U::Matrix{T}
 
733
    S::Vector{Tr}
 
734
    V::Matrix{T}
 
735
end
 
736
 
 
737
factors(F::SVDDense) = (F.U, F.S, F.V)
 
738
 
 
739
function svdfact!{T<:BlasFloat}(A::StridedMatrix{T}, thin::Bool)
 
740
    m,n = size(A)
 
741
    if m == 0 || n == 0
 
742
        u,s,v = (eye(m, thin ? n : m), zeros(0), eye(n,n))
 
743
    else
 
744
        u,s,v = LAPACK.gesdd!(thin ? 'S' : 'A', A)
 
745
    end
 
746
    return SVDDense(u,s,v)
 
747
end
 
748
 
 
749
svdfact!(A::StridedMatrix) = svdfact(A, false)
 
750
 
 
751
svdfact(A::StridedMatrix, thin::Bool) = svdfact!(copy(A), thin)
 
752
svdfact(A::StridedMatrix) = svdfact(A, false)
 
753
 
 
754
function svdvals!(A::StridedMatrix)
 
755
    m,n = size(A)
 
756
    if m == 0 || n == 0
757
757
        return (zeros(T, 0, 0), zeros(T, 0), zeros(T, 0, 0))
758
758
    end
759
 
    if vecs; return LAPACK.gesdd!(thin ? 'S' : 'A', copy(A)); end
760
 
    LAPACK.gesdd!('N', copy(A))
761
 
end
762
 
 
763
 
svdt{T<:Integer}(x::StridedMatrix{T},vecs,thin) = svdt(float64(x),vecs,thin)
764
 
svdt(A) = svdt(A,true,false)
765
 
svdt(A, thin::Bool) = svdt(A,true,thin)
766
 
 
767
 
svdt(x::Number,vecs::Bool,thin::Bool) = vecs ? (x==0?one(x):x/abs(x),abs(x),one(x)) : ([],abs(x),[])
768
 
 
769
 
function svd(x::StridedMatrix,vecs::Bool,thin::Bool) 
770
 
    (u, s, vt) = svdt(x,vecs,thin)
771
 
    return (u, s, vt')
772
 
end
773
 
 
774
 
svd(A) = svd(A,true,false)
775
 
svd(A, thin::Bool) = svd(A,true,thin)
776
 
 
777
 
svdvals(A) = svdt(A,false,true)[2]
 
759
    U, S, V = LAPACK.gesdd!('N', A)
 
760
    return S
 
761
end
 
762
 
 
763
svdvals(A) = svdvals!(copy(A))
 
764
 
 
765
svdt(A::StridedMatrix, thin::Bool) = factors(svdfact(A, thin))
 
766
svdt(A::StridedMatrix) = svdt(A, false)
 
767
svdt(x::Number, thin::Bool) = (x==0?one(x):x/abs(x),abs(x),one(x))
 
768
 
 
769
function svd(A::StridedMatrix, thin::Bool)
 
770
    u,s,v = factors(svdfact(A, thin))
 
771
    return (u,s,v')
 
772
end
 
773
 
 
774
svd(A::StridedMatrix) = svd(A, false)
 
775
svd(x::Number, thin::Bool) = (x==0?one(x):x/abs(x),abs(x),one(x))
 
776
 
 
777
 
 
778
# Generalized svd
 
779
type GSVDDense{T} <: Factorization{T}
 
780
    U::Matrix{T}
 
781
    V::Matrix{T}
 
782
    Q::Matrix{T}
 
783
    a::Vector #{eltype(real(one(T)))}
 
784
    b::Vector #{eltype(real(one(T)))}
 
785
    k::Int
 
786
    l::Int
 
787
    R::Matrix{T}
 
788
end
 
789
 
 
790
function svdfact(A::StridedMatrix, B::StridedMatrix)
 
791
    U, V, Q, a, b, k, l, R = LAPACK.ggsvd!('U', 'V', 'Q', copy(A), copy(B))
 
792
    return GSVDDense(U, V, Q, a, b, int(k), int(l), R)
 
793
end
 
794
 
 
795
svd(A::StridedMatrix, B::StridedMatrix) = factors(svdfact(A, B))
 
796
 
 
797
function factors{T}(obj::GSVDDense{T})
 
798
    m = size(obj.U, 1)
 
799
    p = size(obj.V, 1)
 
800
    n = size(obj.Q, 1)
 
801
    if m - obj.k - obj.l >= 0
 
802
        D1 = [eye(T, obj.k) zeros(T, obj.k, obj.l); zeros(T, obj.l, obj.k) diagm(obj.a[obj.k + 1:obj.k + obj.l]); zeros(T, m - obj.k - obj.l, obj.k + obj.l)]
 
803
        D2 = [zeros(T, obj.l, obj.k) diagm(obj.b[obj.k + 1:obj.k + obj.l]); zeros(T, p - obj.l, obj.k + obj.l)]
 
804
        R0 = [zeros(T, obj.k + obj.l, n - obj.k - obj.l) obj.R]
 
805
    else
 
806
        D1 = [eye(T, m, obj.k) [zeros(T, obj.k, m - obj.k); diagm(a[obj.k + 1:m])] zeros(T, m, obj.k + obj.l - m)]
 
807
        D2 = [zeros(T, p, obj.k) [diagm(b[obj.k + 1:m]); zeros(T, obj.k + p - m, m - obj.k)] [zeros(T, m - obj.k, obj.k + obj.l - m); eye(T, obj.k + p - m, obj.k + obj.l - m)]]
 
808
        R0 = [zeros(T, obj.k + obj.l, n - obj.k - obj.l) obj.R]
 
809
    end
 
810
    return obj.U, obj.V, obj.Q, D1, D2, R0
 
811
end
 
812
 
 
813
svdvals(A::StridedMatrix, B::StridedMatrix) = LAPACK.ggsvd!('N', 'N', 'N', copy(A), copy(B))[4:5]
778
814
 
779
815
schur{T<:BlasFloat}(A::StridedMatrix{T}) = LAPACK.gees!('V', copy(A))
780
816
 
889
925
    elseif p == 1 || p == Inf
890
926
        m, n = size(A)
891
927
        if m != n; error("Use 2-norm for non-square matrices"); end
892
 
        cnd = 1 / LAPACK.gecon!(p == 1 ? '1' : 'I', lud(A).lu, norm(A, p))
 
928
        cnd = 1 / LAPACK.gecon!(p == 1 ? '1' : 'I', lufact(A).lu, norm(A, p))
893
929
    else
894
930
        error("Norm type must be 1, 2 or Inf")
895
931
    end
910
946
end
911
947
 
912
948
SymTridiagonal{T<:BlasFloat}(dv::Vector{T}, ev::Vector{T}) = SymTridiagonal{T}(copy(dv), copy(ev))
 
949
 
913
950
function SymTridiagonal{T<:Real}(dv::Vector{T}, ev::Vector{T})
914
951
    SymTridiagonal{Float64}(float64(dv),float64(ev))
915
952
end
 
953
 
916
954
function SymTridiagonal{Td<:Number,Te<:Number}(dv::Vector{Td}, ev::Vector{Te})
917
955
    T = promote(Td,Te)
918
956
    SymTridiagonal(convert(Vector{T}, dv), convert(Vector{T}, ev))
919
957
end
 
958
 
 
959
SymTridiagonal(A::AbstractMatrix) = SymTridiagonal(diag(A), diag(A,1))
 
960
 
920
961
copy(S::SymTridiagonal) = SymTridiagonal(S.dv,S.ev)
 
962
 
921
963
function full(S::SymTridiagonal)
922
964
    M = diagm(S.dv)
923
965
    for i in 1:length(S.ev)
940
982
 
941
983
eig(m::SymTridiagonal, vecs::Bool) = LAPACK.stev!(vecs ? 'V' : 'N', copy(m.dv), copy(m.ev))
942
984
eig(m::SymTridiagonal) = eig(m::SymTridiagonal, true)
943
 
## This function has been in Julia for some time.  Could probably be dropped.
944
 
trideig{T<:BlasFloat}(d::Vector{T}, e::Vector{T}) = LAPACK.stev!('N', copy(d), copy(e))[1]
 
985
eigvals(m::SymTridiagonal) = eig(m::SymTridiagonal, false)[1]
945
986
 
946
987
## Tridiagonal matrices ##
947
988
type Tridiagonal{T} <: AbstractMatrix{T}
1185
1226
 
1186
1227
#show(io, lu::LUTridiagonal) = print(io, "LU decomposition of ", summary(lu.lu))
1187
1228
 
1188
 
lud!{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
1189
 
lud{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(copy(A.dl),copy(A.d),copy(A.du))...)
1190
 
lu(A::Tridiagonal) = factors(lud(A))
 
1229
lufact!{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(A.dl,A.d,A.du)...)
 
1230
lufact{T}(A::Tridiagonal{T}) = LUTridiagonal{T}(LAPACK.gttrf!(copy(A.dl),copy(A.d),copy(A.du))...)
 
1231
lu(A::Tridiagonal) = factors(lufact(A))
1191
1232
 
1192
1233
function det{T}(lu::LUTridiagonal{T})
1193
1234
    n = length(lu.d)
1194
1235
    prod(lu.d) * (bool(sum(lu.ipiv .!= 1:n) % 2) ? -one(T) : one(T))
1195
1236
end
1196
1237
 
1197
 
det(A::Tridiagonal) = det(lud(A))
 
1238
det(A::Tridiagonal) = det(lufact(A))
1198
1239
 
1199
1240
(\){T<:BlasFloat}(lu::LUTridiagonal{T}, B::StridedVecOrMat{T}) =
1200
1241
    LAPACK.gttrs!('N', lu.dl, lu.d, lu.du, lu.du2, lu.ipiv, copy(B))