diff --git a/base/sparse/linalg.jl b/base/sparse/linalg.jl index 421faeb179483..44e9524e02d6e 100644 --- a/base/sparse/linalg.jl +++ b/base/sparse/linalg.jl @@ -505,12 +505,11 @@ function cond(A::SparseMatrixCSC, p::Real=2) end end -function normestinv{T}(A::SparseMatrixCSC{T}, t::Integer = 2) +function normestinv{T}(A::SparseMatrixCSC{T}, t::Integer = min(2,maximum(size(A)))) maxiter = 5 # Check the input n = max(size(A)...) F = factorize(A) - t = min(t,n) if t <= 0 throw(ArgumentError("number of blocks must be a positive integer")) end @@ -519,10 +518,13 @@ function normestinv{T}(A::SparseMatrixCSC{T}, t::Integer = 2) end ind = Array(Int64, n) ind_hist = Array(Int64, maxiter * t) - S = zeros(T, n, t) + + Ti = typeof(float(zero(T))) + + S = zeros(T <: Real ? Int : Ti, n, t) # Generate the block matrix - X = Array(T, n, t) + X = Array(Ti, n, t) X[1:n,1] = 1 for j = 2:t repeated = true @@ -572,39 +574,41 @@ function normestinv{T}(A::SparseMatrixCSC{T}, t::Integer = 2) S_old = copy(S) for j = 1:t for i = 1:n - S[i,j] = Y[i,j]==0?1:sign(Y[i,j]) + S[i,j] = Y[i,j]==0?one(Y[i,j]):sign(Y[i,j]) end end - # Check wether cols of S are parallel to cols of S or S_old - for j = 1:t - done = false - while ~done - repeated = false - if j > 1 - saux = S[1:n,j]' * S[1:n,1:j-1] - for i = 1:j-1 - if abs(saux[i]) == n - repeated = true - break + if T <: Real + # Check wether cols of S are parallel to cols of S or S_old + for j = 1:t + done = false + while ~done + repeated = false + if j > 1 + saux = S[1:n,j]' * S[1:n,1:j-1] + for i = 1:j-1 + if abs(saux[i]) == n + repeated = true + break + end end end - end - if ~repeated - saux2 = S[1:n,j]' * S_old[1:n,1:t] - for i = 1:t - if abs(saux2[i]) == n - repeated = true - break + if ~repeated + saux2 = S[1:n,j]' * S_old[1:n,1:t] + for i = 1:t + if abs(saux2[i]) == n + repeated = true + break + end end end - end - if repeated - for i = 1:n - S[i,j] = rand()>=0.5?1:-1 + if repeated + for i = 1:n + S[i,j] = rand()>=0.5?1:-1 + end + else + done = true end - else - done = true end end end diff --git a/test/sparsedir/sparse.jl b/test/sparsedir/sparse.jl index 763a28e65a886..079cc797aacfb 100644 --- a/test/sparsedir/sparse.jl +++ b/test/sparsedir/sparse.jl @@ -1079,5 +1079,16 @@ Ar = sprandn(20,20,.5) @test_throws ArgumentError cond(A,2) @test_throws ArgumentError cond(A,3) +# test sparse matrix normestinv +Ac = sprandn(20,20,.5) + im* sprandn(20,20,.5) +Aci = ceil(Int64,100*sprand(20,20,.5))+ im*ceil(Int64,sprand(20,20,.5)) +Ar = sprandn(20,20,.5) +Ari = ceil(Int64,100*Ar) +@test_approx_eq_eps Base.SparseMatrix.normestinv(Ac,3) norm(inv(full(Ac)),1) 1e-4 +@test_approx_eq_eps Base.SparseMatrix.normestinv(Aci,3) norm(inv(full(Aci)),1) 1e-4 +@test_approx_eq_eps Base.SparseMatrix.normestinv(Ar) norm(inv(full(Ar)),1) 1e-4 +@test_throws ArgumentError Base.SparseMatrix.normestinv(Ac,0) +@test_throws ArgumentError Base.SparseMatrix.normestinv(Ac,21) + @test_throws ErrorException transpose(sub(sprandn(10, 10, 0.3), 1:4, 1:4)) @test_throws ErrorException ctranspose(sub(sprandn(10, 10, 0.3), 1:4, 1:4))