Merge pull request #5427 from JuliaLang/anj/symmetrize

Deprecate symmetrize! in favor of copytri!

Author: jiahao

NEWS.md

@@ -136,6 +136,8 @@ Deprecated or removed
   * cholpfact and qrpfact are deprecated in favor of keyword arguments in 
     cholfact and qrfact.
 
+  * symmetrize! is deprecated in favor copytri! but the new function is not exported
+
 [#4042]: https://github.com/JuliaLang/julia/issues/4042
 [#5164]: https://github.com/JuliaLang/julia/issues/5164
 [#4026]: https://github.com/JuliaLang/julia/issues/4026

base/deprecated.jl

@@ -176,6 +176,8 @@ export PipeString
 @deprecate cholpfact!(A,tol=tol)        cholfact!(A, :U, pivot=true, tol=tol)
 @deprecate cholpfact!(A,uplo,tol=tol)   cholfact!(A, uplo, pivot=true, tol=tol)
 @deprecate cholpfact(A)             cholfact(A, :U, pivot=true)
+@deprecate symmetrize!(A)      Base.LinAlg.copytri!(A, 'U')
+@deprecate symmetrize!(A, uplo)      Base.LinAlg.copytri!(A, uplo)
 
 deprecated_ls() = run(`ls -l`)
 deprecated_ls(args::Cmd) = run(`ls -l $args`)

base/exports.jl

@@ -649,7 +649,6 @@ export
     svdfact,
     svdvals!,
     svdvals,
-    symmetrize!,
     trace,
     transpose,
     tril!,

base/linalg.jl

@@ -110,8 +110,6 @@ export
     svdfact,
     svdvals!,
     svdvals,
-    symmetrize!,
-    symmetrize_conj!,
     trace,
     transpose,
     tril,

base/linalg/bunchkaufman.jl

@@ -27,13 +27,13 @@ size(B::BunchKaufman,d::Integer) = size(B.LD,d)
 issym(B::BunchKaufman) = B.symmetric
 ishermitian(B::BunchKaufman) = !B.symmetric
 
-inv{T<:BlasReal}(B::BunchKaufman{T})=symmetrize_conj!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
+inv{T<:BlasReal}(B::BunchKaufman{T})=copytri!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
 
 function inv{T<:BlasComplex}(B::BunchKaufman{T})
     if issym(B)
-        symmetrize!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
+        copytri!(LAPACK.sytri!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
     else
-        symmetrize_conj!(LAPACK.hetri!(B.uplo, copy(B.LD), B.ipiv), B.uplo)
+        copytri!(LAPACK.hetri!(B.uplo, copy(B.LD), B.ipiv), B.uplo, true)
     end
 end
 

base/linalg/dense.jl

@@ -141,8 +141,6 @@ kron(a::Vector, b::Vector)=vec(kron(reshape(a,length(a),1),reshape(b,length(b),1
 kron(a::Matrix, b::Vector)=kron(a,reshape(b,length(b),1))
 kron(a::Vector, b::Matrix)=kron(reshape(a,length(a),1),b)
 
-randsym(n) = symmetrize!(randn(n,n))
-
 ^(A::Matrix, p::Integer) = p < 0 ? inv(A^-p) : Base.power_by_squaring(A,p)
 
 function ^(A::Matrix, p::Number)

base/linalg/factorization.jl

@@ -104,12 +104,12 @@ function logdet{T}(C::Cholesky{T})
     dd + dd # instead of 2.0dd which can change the type
 end
 
-inv(C::Cholesky)=symmetrize_conj!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo)
+inv(C::Cholesky) = copytri!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo, true)
 
 function inv(C::CholeskyPivoted)
     chkfullrank(C)
     ipiv = invperm(C.piv)
-    (symmetrize!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo))[ipiv, ipiv]
+    copytri!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo, true)[ipiv, ipiv]
 end
 
 chkfullrank(C::CholeskyPivoted) = C.rank<size(C.UL, 1) && throw(RankDeficientException(C.info))

base/linalg/hermitian.jl

@@ -17,7 +17,7 @@ function copy!(A::Hermitian, B::Hermitian)
 end
 size(A::Hermitian, args...) = size(A.S, args...)
 getindex(A::Hermitian, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.S, i, j) : conj(getindex(A.S, j, i))
-full(A::Hermitian) = symmetrize_conj!(A.S, A.uplo)
+full(A::Hermitian) = copytri!(A.S, A.uplo, true)
 ishermitian(A::Hermitian) = true
 issym{T<:Real}(A::Hermitian{T}) = true
 issym{T<:Complex}(A::Hermitian{T}) = all(imag(A.S) .== 0)
@@ -57,10 +57,10 @@ function sqrtm(A::Hermitian, cond::Bool=false)
     length(F[:values]) == 0 && return A
     vsqrt = sqrt(complex(F[:values]))
     if all(imag(vsqrt) .== 0)
-        retmat = symmetrize!(scale(F[:vectors], real(vsqrt)) * F[:vectors]')
+        retmat = copytri!(scale(F[:vectors], real(vsqrt)) * F[:vectors]', 'U')
     else
         zc = complex(F[:vectors])
-        retmat = symmetrize!(scale(zc, vsqrt) * zc')
+        retmat = copytri!(scale(zc, vsqrt) * zc', 'U')
     end
     return cond ? (retmat, norm(vsqrt, Inf)^2/norm(F[:values], Inf)) : retmat
 end

base/linalg/matmul.jl

@@ -141,32 +141,19 @@ Ac_mul_Bt!{T,S,R}(C::StridedMatrix{R}, A::StridedMatrix{T}, B::StridedMatrix{S})
 
 # Supporting functions for matrix multiplication
 
-function symmetrize!(A::StridedMatrix, uplo::Char='U')
+function copytri!(A::StridedMatrix, uplo::Char, conjugate::Bool=false)
     n = chksquare(A)
     @chkuplo
     if uplo == 'U'
         for i = 1:(n-1), j = (i+1):n
-            A[j,i] = A[i,j]
+            A[j,i] = conjugate ? conj(A[i,j]) : A[i,j]
         end
-    else #uplo == 'L'
+    elseif uplo == 'L'
         for i = 1:(n-1), j = (i+1):n
-            A[i,j] = A[j,i]
-        end
-    end
-    A
-end
-
-function symmetrize_conj!(A::StridedMatrix, uplo::Char='U')
-    n = chksquare(A)
-    @chkuplo
-    if uplo == 'U'
-        for i = 1:(n-1), j = (i+1):n
-            A[j,i] = conj(A[i,j])
-        end
-    else #uplo == 'L'
-        for i = 1:(n-1), j = (i+1):n
-            A[i,j] = conj(A[j,i])
+            A[i,j] = conjugate ? conj(A[j,i]) : A[j,i]
         end
+    else
+        throw(ArgumentError("second argument must be 'U' or 'L'"))
     end
     A
 end
@@ -200,7 +187,7 @@ function syrk_wrapper{T<:BlasFloat}(tA::Char, A::StridedMatrix{T})
     if mA == 3 && nA == 3; return matmul3x3(tA,tAt,A,A); end
 
     stride(A, 1) == 1 || (return generic_matmatmul(tA, tAt, A, A))
-    symmetrize!(BLAS.syrk('U', tA, one(T), A))
+    copytri!(BLAS.syrk('U', tA, one(T), A), 'U')
 end
 
 function herk_wrapper{T<:BlasFloat}(tA::Char, A::StridedMatrix{T})
@@ -220,7 +207,7 @@ function herk_wrapper{T<:BlasFloat}(tA::Char, A::StridedMatrix{T})
     # Result array does not need to be initialized as long as beta==0
     #    C = Array(T, mA, mA)
 
-    symmetrize_conj!(BLAS.herk('U', tA, one(T), A))
+    copytri!(BLAS.herk('U', tA, one(T), A), 'U', true)
 end
 
 function gemm_wrapper{T<:BlasFloat}(tA::Char, tB::Char,

base/linalg/symmetric.jl

@@ -17,7 +17,7 @@ function copy!(A::Symmetric, B::Symmetric)
 end
 size(A::Symmetric, args...) = size(A.S, args...)
 getindex(A::Symmetric, i::Integer, j::Integer) = (A.uplo == 'U') == (i < j) ? getindex(A.S, i, j) : getindex(A.S, j, i)
-full(A::Symmetric) = symmetrize!(A.S, A.uplo)
+full(A::Symmetric) = copytri!(A.S, A.uplo)
 ishermitian{T<:Real}(A::Symmetric{T}) = true
 ishermitian{T<:Complex}(A::Symmetric{T}) = all(imag(A.S) .== 0)
 issym(A::Symmetric) = true
@@ -57,10 +57,10 @@ function sqrtm{T<:Real}(A::Symmetric{T}, cond::Bool)
     F = eigfact(A)
     vsqrt = sqrt(complex(F[:values]))
     if all(imag(vsqrt) .== 0)
-        retmat = symmetrize!(scale(F[:vectors], real(vsqrt)) * F[:vectors]')
+        retmat = copytri!(scale(F[:vectors], real(vsqrt)) * F[:vectors]', 'U')
     else
         zc = complex(F[:vectors])
-        retmat = symmetrize!(scale(zc, vsqrt) * zc')
+        retmat = copytri!(scale(zc, vsqrt) * zc', 'U')
     end
     cond ? (retmat, norm(vsqrt, Inf)^2/norm(F[:values], Inf)) : retmat
 end

base/linalg/triangular.jl

@@ -89,8 +89,8 @@ getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? (A.unitdiag ==
 istril(A::Triangular) = A.uplo == 'L' || istriu(A.UL)
 istriu(A::Triangular) = A.uplo == 'U' || istril(A.UL)
 
-transpose(A::Triangular) = Triangular(symmetrize!(A.UL, A.uplo), A.uplo=='U'?'L':'U', A.unitdiag)
-ctranspose(A::Triangular) = Triangular(symmetrize_conj!(A.UL, A.uplo), A.uplo=='U'?'L':'U', A.unitdiag)
+transpose(A::Triangular) = Triangular(copytri!(A.UL, A.uplo), A.uplo=='U'?'L':'U', A.unitdiag)
+ctranspose(A::Triangular) = Triangular(copytri!(A.UL, A.uplo, true), A.uplo=='U'?'L':'U', A.unitdiag)
 diag(A::Triangular) = diag(A.UL)
 big(A::Triangular) = Triangular(big(A.UL), A.uplo, A.unitdiag)
 

doc/stdlib/linalg.rst

@@ -246,12 +246,6 @@ Linear algebra functions in Julia are largely implemented by calling functions f
    ``A`` by ``b[i]`` (similar to ``diagm(b)*A``), again operating in-place
    on ``A``.
 
-.. function:: symmetrize!(A[, UL::Char])
-
-   ``symmetrize!(A)`` converts from the BLAS/LAPACK symmetric storage
-   format, in which only the ``UL`` ('U'pper or 'L'ower, default 'U')
-   triangle is used, to a full symmetric matrix.
-
 .. function:: Tridiagonal(dl, d, du)
 
    Construct a tridiagonal matrix from the lower diagonal, diagonal, and upper diagonal, respectively.  The result is of type ``Tridiagonal`` and provides efficient specialized linear solvers, but may be converted into a regular matrix with ``full``.

test/blas.jl

@@ -77,7 +77,7 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @test all(tril(BLAS.herk('L', 'C', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
         ans = similar(L4)
         @test all(tril(BLAS.herk('L','C', L4)) .== tril(BLAS.herk!('L', 'C', one(elty), L4, zero(elty), ans)))
-        @test all(symmetrize!(ans, 'L') .== LinAlg.BLAS.gemm('T', 'N', L4, L4))
+        @test all(Base.LinAlg.copytri!(ans, 'L') .== LinAlg.BLAS.gemm('T', 'N', L4, L4))
     else
         @test all(triu(BLAS.syrk('U', 'N', U4)) .== triu(BLAS.gemm('N', 'T', U4, U4)))
         @test all(tril(BLAS.syrk('L', 'N', U4)) .== tril(BLAS.gemm('N', 'T', U4, U4)))
@@ -89,6 +89,6 @@ for elty in (Float32, Float64, Complex64, Complex128)
         @test all(tril(BLAS.syrk('L', 'T', L4)) .== tril(BLAS.gemm('T', 'N', L4, L4)))
         ans = similar(L4)
         @test all(tril(BLAS.syrk('L','T', L4)) .== tril(BLAS.syrk!('L', 'T', one(elty), L4, zero(elty), ans)))
-        @test all(symmetrize!(ans, 'L') .== BLAS.gemm('T', 'N', L4, L4))
+        @test all(Base.LinAlg.copytri!(ans, 'L') .== BLAS.gemm('T', 'N', L4, L4))
     end
 end