Merge pull request #1867 from andreasnoackjensen/cond

RFC: Use LAPACK to estimate 1 and Inf norm. Change norm to return Float.

base/lapack.jl

@@ -1519,4 +1519,75 @@ for (syevr, elty, relty) in
 end
 syevr!(A::StridedMatrix, Z::StridedMatrix) = syevr!('V', 'A', 'U', A, 0.0, 0.0, 0, 0, Z, -1.0)
 syevr!{T}(A::StridedMatrix{T}) = syevr!('N', 'A', 'U', A, 0.0, 0.0, 0, 0, zeros(T,0,0), -1.0)
+
+# Estimate condition number
+for (gecon, elty) in
+ ((:dgecon_,:Float64),
+ (:sgecon_,:Float32))
+ @eval begin
+ function gecon!(normtype::BlasChar, A::StridedMatrix{$elty}, anorm::$elty)
+# SUBROUTINE DGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, IWORK,
+# $ INFO )
+# * .. Scalar Arguments ..
+# CHARACTER NORM
+# INTEGER INFO, LDA, N
+# DOUBLE PRECISION ANORM, RCOND
+# * ..
+# * .. Array Arguments ..
+# INTEGER IWORK( * )
+# DOUBLE PRECISION A( LDA, * ), WORK( * )
+ chkstride1(A)
+ n = size(A, 2)
+ lda = max(1, size(A, 1))
+ rcond = Array($elty, 1)
+ work = Array($elty, 4n)
+ iwork = Array(BlasInt, n)
+ info = Array(BlasInt, 1)
+ ccall(($(string(gecon)),liblapack), Void,
+ (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,
+ info)
+ if info[1] != 0 throw(LapackException(info[1])) end
+ return rcond[1]
+ end
+ end
+end
+for (gecon, elty, relty) in
+ ((:zgecon_,:Complex128,:Float64),
+ (:cgecon_,:Complex64,:Float32))
+ @eval begin
+ function gecon!(normtype::BlasChar, A::StridedMatrix{$elty}, anorm::$relty)
+ chkstride1(A)
+# SUBROUTINE ZGECON( NORM, N, A, LDA, ANORM, RCOND, WORK, RWORK,
+# $ INFO )
+# * .. Scalar Arguments ..
+# CHARACTER NORM
+# INTEGER INFO, LDA, N
+# DOUBLE PRECISION ANORM, RCOND
+# * ..
+# * .. Array Arguments ..
+# DOUBLE PRECISION RWORK( * )
+# COMPLEX*16 A( LDA, * ), WORK( * )
+ chkstride1(A)
+ n = size(A, 2)
+ lda = max(1, size(A, 1))
+ rcond = Array($relty, 1)
+ work = Array($elty, 2n)
+ rwork = Array(BlasInt, 2n)
+ info = Array(BlasInt, 1)
+ ccall(($(string(gecon)),liblapack), Void,
+ (Ptr{Uint8}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, 
+ Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{$relty},
+ Ptr{BlasInt}),
+ &normtype, &n, A, &lda, 
+ &anorm, rcond, work, rwork,
+ info)
+ if info[1] < 0 throw(LapackException(info[1])) end
+ return rcond[1]
+ end
+ end
+end
 end # module

base/linalg.jl

@@ -35,42 +35,44 @@ diag(A::AbstractVector) = error("Perhaps you meant to use diagm().")
 
 function norm{T}(x::AbstractVector{T}, p::Number)
  if length(x) == 0
- return zero(eltype(x))
+ a = zero(eltype(x))
  elseif p == Inf
- return max(abs(x))
+ a = max(abs(x))
  elseif p == -Inf
- return min(abs(x))
+ a = min(abs(x))
  else
  absx = abs(x)
  dx = max(absx)
  if dx != zero(T)
  scale!(absx, 1/dx)
- return dx * (sum(absx.^p).^(1/p))
+ a = dx * (sum(absx.^p).^(1/p))
  else
- return sum(absx.^p).^(1/p)
+ a = sum(absx.^p).^(1/p)
  end
  end
+ return float(a)
 end
 norm{T<:Integer}(x::AbstractVector{T}, p::Number) = norm(float(x), p)
 norm(x::AbstractVector) = norm(x, 2)
 
 function norm(A::AbstractMatrix, p)
  m, n = size(A)
  if m == 0 || n == 0
- return zero(eltype(A))
+ a = zero(eltype(A))
  elseif m == 1 || n == 1
- return norm(reshape(A, numel(A)), p)
+ a = norm(reshape(A, numel(A)), p)
  elseif p == 1
- return max(sum(abs(A),1))
+ a = max(sum(abs(A),1))
  elseif p == 2
- return max(svdvals(A))
+ a = max(svdvals(A))
  elseif p == Inf
- return max(sum(abs(A),2))
+ a = max(sum(abs(A),2))
  elseif p == :fro
- return norm(reshape(A, numel(A)))
+ a = norm(reshape(A, numel(A)))
  else
  error("invalid parameter to matrix norm")
  end
+ return float(a)
 end
 
 norm(A::AbstractMatrix) = norm(A, 2)
@@ -96,26 +98,7 @@ trace(x::Number) = x
 #det(a::AbstractMatrix)
 inv(a::AbstractMatrix) = a \ one(a)
 
-function cond(a::AbstractMatrix)
- s = svdvals(a)
- condno = max(s) / min(s)
- # Return Inf if condno is NaN (input is all zeros)
- isnan(condno) ? Inf : condno
-end
-
-function cond(a::AbstractMatrix, p) 
- if p == 2 
- return cond(a)
- else
- try
- return norm(a, p) * norm(inv(a), p)
- catch e
- isa(e,LAPACK.SingularException) ? (return Inf) : rethrow(e)
- end
- end
-end
-
-cond(x::Number) = x == 0 ? Inf : 1
+cond(x::Number) = x == 0 ? Inf : 1.0
 cond(x::Number, p) = cond(x)
 
 function issym(A::AbstractMatrix)

base/linalg_dense.jl

@@ -809,6 +809,22 @@ end
 null{T<:Integer}(A::StridedMatrix{T}) = null(float(A))
 null(a::StridedVector) = null(reshape(a, length(a), 1))