Merge pull request #12592 from marcusps/normest

Fix complex matrix support in normestinv
JuliaLang · Aug 20, 2015 · 4dbab5f · 4dbab5f
2 parents 7420be9 + 6b7cb08
commit 4dbab5f
Show file tree

Hide file tree

Showing 2 changed files with 44 additions and 29 deletions.
diff --git 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
@@ -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))