Add Stewart's fast orthogonal/unitary matrix generation algorithm

Author: andreasnoack

src/HaarMeasure.jl

@@ -70,3 +70,44 @@ end
 #Zyczkowski and Kus, Random unitary matrices, J. Phys. A: Math. Gen. 27,
 #4235–4245 (1994).
 
+### Stewarts algorithm for n^2 orthogonal random matrix
+using Base.LinAlg: BlasInt
+for (s, elty) in (("dlarfg_", Float64),
+                  ("zlarfg_", Complex128))
+    @eval begin
+        function larfg!(n::Int, α::Ptr{$elty}, x::Ptr{$elty}, incx::Int, τ::Ptr{$elty})
+            ccall(($(Base.blasfunc(s)), Base.liblapack_name), Void,
+                (Ptr{BlasInt}, Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}),
+                &n, α, x, &incx, τ)
+        end
+    end
+end
+
+function Stewart(::Type{Float64}, n)
+    τ = Array(Float64, n)
+    H = randn(n, n)
+
+    pτ = pointer(τ)
+    pβ = pointer(H)
+    pH = pointer(H, 2)
+
+    for i = 0:n-2
+        larfg!(n - i, pβ + (n + 1)*8i, pH + (n + 1)*8i, 1, pτ + 8i)
+    end
+    LinAlg.QRPackedQ(H,τ)
+end
+function Stewart(::Type{Complex128}, n)
+    τ = Array(Complex128, n)
+    H = complex(randn(n, n), randn(n,n))
+
+    pτ = pointer(τ)
+    pβ = pointer(H)
+    pH = pointer(H, 2)
+
+    for i = 0:n-2
+        larfg!(n - i, pβ + (n + 1)*16i, pH + (n + 1)*16i, 1, pτ + 16i)
+    end
+    LinAlg.QRPackedQ(H,τ)
+end
+
+

test/Haar.jl

@@ -12,4 +12,9 @@ Q=UniformHaar(2, N)
 println("Case 3")
 println("E(A*Q*B*Q'*A*Q*B*Q') = ", eval(expectation(N, :Q, :(A*Q*B*Q'*A*Q*B*Q'))))
 
+for elty in (Float64, Complex128)
+	A = Stewart(elty, N)
+	@test_approx_eq A'A eye(N)
+end
+
 end #_HAVE_GSL