reset GLOBAL_RNG state in/out at-testset for reproducibility

This is a follow-up on the recently introduced `guardsrand` functionality,
first suggested at
#16940 (comment).
Each `testset` block is now implicitly wrapped in a `guardsrand` block, which
is more user-friendly and more systematic (only few users know about
`guardsrand`).

These are essentially two new features:

1) "in": in base, tests are run with the global RNG randomly seeded,
   but the seed is printed in case of failure to allow reproducing the
   failure; but even if the failure occurs after many tests, when they
   alter the global RNG state, one has to re-run the whole file, which
   can be time-consuming; with this change, at the beginning of each
   `testset`, the global RNG is re-seeded with its own seed: this
   allows to re-run only the failing `testset`, which can be done
   easily in the REPL (the seeding occurs for each loop in a `testset for`;
   this also allows to re-arrange `testset`s in arbitrary order
   w.r.t. the global RNG;

2) "out": a `testset` leaves no tracks of its use of `srand` or `rand`
   (this "feature" should be less and less needed/useful with the
   generalization of the use of `testset` blocks).

Example:
```
@testset begin
   srand(123)
   rand()
   @testset for T in (Int, Float64)
       # here is an implicit srand(123), by 1)
       rand()
   end
   rand() # this value will not be affected if the sub-`testset` block
          # above is removed, or if another rand() call is added therein, by 2)
end
```

Note that guardsrand can't be used directly, as then the testset's body is
wrapped in a function, which causes problem with overwriting loop variable
("outer"), using `using`, defining new methods, etc. So we need to duplicate
guardsrand's logic.

Author: rfourquet

stdlib/IterativeEigensolvers/test/runtests.jl

@@ -4,118 +4,118 @@ using IterativeEigensolvers
 using Test
 
 @testset "eigs" begin
-    guardsrand(1234) do
-        n = 10
-        areal  = sprandn(n,n,0.4)
-        breal  = sprandn(n,n,0.4)
-        acmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4))
-        bcmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4))
+    srand(1234)
+    n = 10
+    areal  = sprandn(n,n,0.4)
+    breal  = sprandn(n,n,0.4)
+    acmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4))
+    bcmplx = complex.(sprandn(n,n,0.4), sprandn(n,n,0.4))
 
-        testtol = 1e-6
+    testtol = 1e-6
 
-        @testset for elty in (Float64, ComplexF64)
-            if elty == ComplexF32 || elty == ComplexF64
-                a = acmplx
-                b = bcmplx
-            else
-                a = areal
-                b = breal
-            end
-            a_evs = eigvals(Array(a))
-            a     = convert(SparseMatrixCSC{elty}, a)
-            asym  = a' + a                  # symmetric indefinite
-            apd   = a'*a                    # symmetric positive-definite
+    @testset for elty in (Float64, ComplexF64)
+        if elty == ComplexF32 || elty == ComplexF64
+            a = acmplx
+            b = bcmplx
+        else
+            a = areal
+            b = breal
+        end
+        a_evs = eigvals(Array(a))
+        a     = convert(SparseMatrixCSC{elty}, a)
+        asym  = a' + a                  # symmetric indefinite
+        apd   = a'*a                    # symmetric positive-definite
 
-            b     = convert(SparseMatrixCSC{elty}, b)
-            bsym  = b' + b
-            bpd   = b'*b
+        b     = convert(SparseMatrixCSC{elty}, b)
+        bsym  = b' + b
+        bpd   = b'*b
 
-            (d,v) = eigs(a, nev=3)
-            @test a*v[:,2] ≈ d[2]*v[:,2]
-            @test norm(v) > testtol # eigenvectors cannot be null vectors
-            (d,v) = eigs(a, I, nev=3) # test eigs(A, B; kwargs...)
-            @test a*v[:,2] ≈ d[2]*v[:,2]
-            @test norm(v) > testtol # eigenvectors cannot be null vectors
-            @test_logs (:warn,"Use symbols instead of strings for specifying which eigenvalues to compute") eigs(a, which="LM")
-            @test_logs (:warn,"Adjusting ncv from 1 to 4") eigs(a, ncv=1, nev=2)
-            @test_logs (:warn,"Adjusting nev from $n to $(n-2)") eigs(a, nev=n)
-            # (d,v) = eigs(a, b, nev=3, tol=1e-8) # not handled yet
-            # @test a*v[:,2] ≈ d[2]*b*v[:,2] atol=testtol
-            # @test norm(v) > testtol # eigenvectors cannot be null vectors
-            if elty <: Base.LinAlg.BlasComplex
-                sr_ind = indmin(real.(a_evs))
-                (d, v) = eigs(a, nev=1, which=:SR)
-                @test d[1] ≈ a_evs[sr_ind]
-                si_ind = indmin(imag.(a_evs))
-                (d, v) = eigs(a, nev=1, which=:SI)
-                @test d[1] ≈ a_evs[si_ind]
-                lr_ind = indmax(real.(a_evs))
-                (d, v) = eigs(a, nev=1, which=:LR)
-                @test d[1] ≈ a_evs[lr_ind]
-                li_ind = indmax(imag.(a_evs))
-                (d, v) = eigs(a, nev=1, which=:LI)
-                @test d[1] ≈ a_evs[li_ind]
-            end
+        (d,v) = eigs(a, nev=3)
+        @test a*v[:,2] ≈ d[2]*v[:,2]
+        @test norm(v) > testtol # eigenvectors cannot be null vectors
+        (d,v) = eigs(a, I, nev=3) # test eigs(A, B; kwargs...)
+        @test a*v[:,2] ≈ d[2]*v[:,2]
+        @test norm(v) > testtol # eigenvectors cannot be null vectors
+        @test_logs (:warn,"Use symbols instead of strings for specifying which eigenvalues to compute") eigs(a, which="LM")
+        @test_logs (:warn,"Adjusting ncv from 1 to 4") eigs(a, ncv=1, nev=2)
+        @test_logs (:warn,"Adjusting nev from $n to $(n-2)") eigs(a, nev=n)
+        # (d,v) = eigs(a, b, nev=3, tol=1e-8) # not handled yet
+        # @test a*v[:,2] ≈ d[2]*b*v[:,2] atol=testtol
+        # @test norm(v) > testtol # eigenvectors cannot be null vectors
+        if elty <: Base.LinAlg.BlasComplex
+            sr_ind = indmin(real.(a_evs))
+            (d, v) = eigs(a, nev=1, which=:SR)
+            @test d[1] ≈ a_evs[sr_ind]
+            si_ind = indmin(imag.(a_evs))
+            (d, v) = eigs(a, nev=1, which=:SI)
+            @test d[1] ≈ a_evs[si_ind]
+            lr_ind = indmax(real.(a_evs))
+            (d, v) = eigs(a, nev=1, which=:LR)
+            @test d[1] ≈ a_evs[lr_ind]
+            li_ind = indmax(imag.(a_evs))
+            (d, v) = eigs(a, nev=1, which=:LI)
+            @test d[1] ≈ a_evs[li_ind]
+        end
 
-            (d,v) = eigs(asym, nev=3)
-            @test asym*v[:,1] ≈ d[1]*v[:,1]
-            @test eigs(asym; nev=1, sigma=d[3])[1][1] ≈ d[3]
-            @test norm(v) > testtol # eigenvectors cannot be null vectors
+        (d,v) = eigs(asym, nev=3)
+        @test asym*v[:,1] ≈ d[1]*v[:,1]
+        @test eigs(asym; nev=1, sigma=d[3])[1][1] ≈ d[3]
+        @test norm(v) > testtol # eigenvectors cannot be null vectors
 
-            (d,v) = eigs(apd, nev=3)
-            @test apd*v[:,3] ≈ d[3]*v[:,3]
-            @test eigs(apd; nev=1, sigma=d[3])[1][1] ≈ d[3]
+        (d,v) = eigs(apd, nev=3)
+        @test apd*v[:,3] ≈ d[3]*v[:,3]
+        @test eigs(apd; nev=1, sigma=d[3])[1][1] ≈ d[3]
 
-            (d,v) = eigs(apd, bpd, nev=3, tol=1e-8)
-            @test apd*v[:,2] ≈ d[2]*bpd*v[:,2] atol=testtol
-            @test norm(v) > testtol # eigenvectors cannot be null vectors
+        (d,v) = eigs(apd, bpd, nev=3, tol=1e-8)
+        @test apd*v[:,2] ≈ d[2]*bpd*v[:,2] atol=testtol
+        @test norm(v) > testtol # eigenvectors cannot be null vectors
 
-            @testset "(shift-and-)invert mode" begin
-                (d,v) = eigs(apd, nev=3, sigma=0)
-                @test apd*v[:,3] ≈ d[3]*v[:,3]
-                @test norm(v) > testtol # eigenvectors cannot be null vectors
+        @testset "(shift-and-)invert mode" begin
+            (d,v) = eigs(apd, nev=3, sigma=0)
+            @test apd*v[:,3] ≈ d[3]*v[:,3]
+            @test norm(v) > testtol # eigenvectors cannot be null vectors
 
-                (d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8)
-                @test apd*v[:,1] ≈ d[1]*bpd*v[:,1] atol=testtol
-                @test norm(v) > testtol # eigenvectors cannot be null vectors
-            end
+            (d,v) = eigs(apd, bpd, nev=3, sigma=0, tol=1e-8)
+            @test apd*v[:,1] ≈ d[1]*bpd*v[:,1] atol=testtol
+            @test norm(v) > testtol # eigenvectors cannot be null vectors
+        end
 
-            @testset "ArgumentErrors" begin
-                @test_throws ArgumentError eigs(rand(elty,2,2))
-                @test_throws ArgumentError eigs(a, nev=-1)
-                @test_throws ArgumentError eigs(a, which=:Z)
-                @test_throws ArgumentError eigs(a, which=:BE)
-                @test_throws DimensionMismatch eigs(a, v0=zeros(elty,n+2))
-                @test_throws ArgumentError eigs(a, v0=zeros(Int,n))
-                if elty == Float64
-                    @test_throws ArgumentError eigs(a+a.',which=:SI)
-                    @test_throws ArgumentError eigs(a+a.',which=:LI)
-                    @test_throws ArgumentError eigs(a,sigma=rand(ComplexF32))
-                end
+        @testset "ArgumentErrors" begin
+            @test_throws ArgumentError eigs(rand(elty,2,2))
+            @test_throws ArgumentError eigs(a, nev=-1)
+            @test_throws ArgumentError eigs(a, which=:Z)
+            @test_throws ArgumentError eigs(a, which=:BE)
+            @test_throws DimensionMismatch eigs(a, v0=zeros(elty,n+2))
+            @test_throws ArgumentError eigs(a, v0=zeros(Int,n))
+            if elty == Float64
+                @test_throws ArgumentError eigs(a+a.',which=:SI)
+                @test_throws ArgumentError eigs(a+a.',which=:LI)
+                @test_throws ArgumentError eigs(a,sigma=rand(ComplexF32))
             end
         end
+    end
 
-        @testset "Symmetric generalized with singular B" begin
-            n = 10
-            k = 3
-            A = randn(n,n); A = A'A
-            B = randn(n,k); B = B*B'
-            @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k])
-        end
+    @testset "Symmetric generalized with singular B" begin
+        srand(124)
+        n = 10
+        k = 3
+        A = randn(n,n); A = A'A
+        B = randn(n,k); B = B*B'
+        @test sort(eigs(A, B, nev = k, sigma = 1.0)[1]) ≈ sort(eigvals(A, B)[1:k])
     end
 end
 
 # Problematic example from #6965A
 let A6965 = [
-         1.0   1.0   1.0   1.0   1.0   1.0   1.0  1.0
-        -1.0   2.0   0.0   0.0   0.0   0.0   0.0  1.0
-        -1.0   0.0   3.0   0.0   0.0   0.0   0.0  1.0
-        -1.0   0.0   0.0   4.0   0.0   0.0   0.0  1.0
-        -1.0   0.0   0.0   0.0   5.0   0.0   0.0  1.0
-        -1.0   0.0   0.0   0.0   0.0   6.0   0.0  1.0
-        -1.0   0.0   0.0   0.0   0.0   0.0   7.0  1.0
-        -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  8.0
-       ]
+    1.0   1.0   1.0   1.0   1.0   1.0   1.0  1.0
+    -1.0   2.0   0.0   0.0   0.0   0.0   0.0  1.0
+    -1.0   0.0   3.0   0.0   0.0   0.0   0.0  1.0
+    -1.0   0.0   0.0   4.0   0.0   0.0   0.0  1.0
+    -1.0   0.0   0.0   0.0   5.0   0.0   0.0  1.0
+    -1.0   0.0   0.0   0.0   0.0   6.0   0.0  1.0
+    -1.0   0.0   0.0   0.0   0.0   0.0   7.0  1.0
+    -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  -1.0  8.0
+]
     d, = eigs(A6965,which=:SM,nev=2,ncv=4,tol=eps())
     @test d[1] ≈ 2.5346936860350002
     @test real(d[2]) ≈ 2.6159972444834976

stdlib/Test/src/Test.jl

@@ -964,12 +964,20 @@ function testset_beginend(args, tests, source)
         # which is needed for backtrace scrubbing to work correctly.
         while false; end
         push_testset(ts)
+        # we reproduce the logic of guardsrand, but this function
+        # cannot be used as it changes slightly the semantic of @testset,
+        # by wrapping the body in a function
+        oldrng = copy(Base.GLOBAL_RNG)
         try
+            # GLOBAL_RNG is re-seeded with its own seed to ease reproduce a failed test
+            srand(Base.GLOBAL_RNG.seed)
             $(esc(tests))
         catch err
             # something in the test block threw an error. Count that as an
             # error in this test set
             record(ts, Error(:nontest_error, :(), err, catch_backtrace(), $(QuoteNode(source))))
+        finally
+            copy!(Base.GLOBAL_RNG, oldrng)
         end
         pop_testset()
         finish(ts)
@@ -1030,12 +1038,16 @@ function testset_forloop(args, testloop, source)
         ts = $(testsettype)($desc; $options...)
         push_testset(ts)
         first_iteration = false
+        oldrng = copy(Base.GLOBAL_RNG)
         try
+            srand(Base.GLOBAL_RNG.seed)
             $(esc(tests))
         catch err
             # Something in the test block threw an error. Count that as an
             # error in this test set
             record(ts, Error(:nontest_error, :(), err, catch_backtrace(), $(QuoteNode(source))))
+        finally
+            copy!(Base.GLOBAL_RNG, oldrng)
         end
     end
     quote
@@ -1477,7 +1489,7 @@ end
 
 "`guardsrand(f, seed)` is equivalent to running `srand(seed); f()` and
 then restoring the state of the global RNG as it was before."
-guardsrand(f::Function, seed::Integer) = guardsrand() do
+guardsrand(f::Function, seed::Union{Vector{UInt32},Integer}) = guardsrand() do
     srand(seed)
     f()
 end

stdlib/Test/test/runtests.jl

@@ -287,10 +287,8 @@ end
                 end
             end
             @testset "some loops fail" begin
-                guardsrand(123) do
-                    @testset for i in 1:5
-                        @test i <= rand(1:10)
-                    end
+                @testset for i in 1:5
+                    @test i <= 4
                 end
                 # should add 3 errors and 3 passing tests
                 @testset for i in 1:6
@@ -739,3 +737,22 @@ end
                  be tested in --depwarn=error mode"""
     end
 end
+
+@testset "@testset preserves GLOBAL_RNG's state, and re-seeds it" begin
+    # i.e. it behaves as if it was wrapped in a `guardsrand(GLOBAL_RNG.seed)` block
+    seed = rand(UInt128)
+    srand(seed)
+    a = rand()
+    @testset begin
+        # global RNG must re-seeded at the beginning of @testset
+        @test a == rand()
+    end
+    @testset for i=1:3
+        @test a == rand()
+    end
+    # the @testset's above must have no consequence for rand() below
+    b = rand()
+    srand(seed)
+    @test a == rand()
+    @test b == rand()
+end

test/linalg/lapack.jl

@@ -10,25 +10,24 @@ import Base.LinAlg.BlasInt
 @test_throws ArgumentError Base.LinAlg.LAPACK.chktrans('Z')
 
 @testset "syevr" begin
-    guardsrand(123) do
-        Ainit = randn(5,5)
-        @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
-            if elty == ComplexF32 || elty == ComplexF64
-                A = complex.(Ainit, Ainit)
-            else
-                A = Ainit
-            end
-            A = convert(Array{elty, 2}, A)
-            Asym = A'A
-            vals, Z = LAPACK.syevr!('V', copy(Asym))
-            @test Z*(Diagonal(vals)*Z') ≈ Asym
-            @test all(vals .> 0.0)
-            @test LAPACK.syevr!('N','V','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[vals .< 1.0]
-            @test LAPACK.syevr!('N','I','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[4:5]
-            @test vals ≈ LAPACK.syev!('N','U',copy(Asym))
-
-            @test_throws DimensionMismatch LAPACK.sygvd!(1,'V','U',copy(Asym),ones(elty,6,6))
+    srand(123)
+    Ainit = randn(5,5)
+    @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
+        if elty == ComplexF32 || elty == ComplexF64
+            A = complex.(Ainit, Ainit)
+        else
+            A = Ainit
         end
+        A = convert(Array{elty, 2}, A)
+        Asym = A'A
+        vals, Z = LAPACK.syevr!('V', copy(Asym))
+        @test Z*(Diagonal(vals)*Z') ≈ Asym
+        @test all(vals .> 0.0)
+        @test LAPACK.syevr!('N','V','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[vals .< 1.0]
+        @test LAPACK.syevr!('N','I','U',copy(Asym),0.0,1.0,4,5,-1.0)[1] ≈ vals[4:5]
+        @test vals ≈ LAPACK.syev!('N','U',copy(Asym))
+
+        @test_throws DimensionMismatch LAPACK.sygvd!(1,'V','U',copy(Asym),ones(elty,6,6))
     end
 end
 
@@ -207,12 +206,11 @@ end
 
 @testset "gels" begin
     @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
-        guardsrand(913) do
-            A = rand(elty,10,10)
-            X = rand(elty,10)
-            B,Y,z = LAPACK.gels!('N',copy(A),copy(X))
-            @test A\X ≈ Y
-        end
+        srand(913)
+        A = rand(elty,10,10)
+        X = rand(elty,10)
+        B,Y,z = LAPACK.gels!('N',copy(A),copy(X))
+        @test A\X ≈ Y
     end
 end
 
@@ -435,36 +433,34 @@ end
 
 @testset "sysv" begin
     @testset for elty in (Float32, Float64, ComplexF32, ComplexF64)
-        guardsrand(123) do
-            A = rand(elty,10,10)
-            A = A + A.' #symmetric!
-            b = rand(elty,10)
-            c = A \ b
-            b,A = LAPACK.sysv!('U',A,b)
-            @test b ≈ c
-            @test_throws DimensionMismatch LAPACK.sysv!('U',A,rand(elty,11))
-        end
+        srand(123)
+        A = rand(elty,10,10)
+        A = A + A.' #symmetric!
+        b = rand(elty,10)
+        c = A \ b
+        b,A = LAPACK.sysv!('U',A,b)
+        @test b ≈ c
+        @test_throws DimensionMismatch LAPACK.sysv!('U',A,rand(elty,11))
     end
 end
 
 @testset "hesv" begin
     @testset for elty in (ComplexF32, ComplexF64)
-        guardsrand(935) do
-            A = rand(elty,10,10)
-            A = A + A' #hermitian!
-            b = rand(elty,10)
-            c = A \ b
-            b,A = LAPACK.hesv!('U',A,b)
-            @test b ≈ c
-            @test_throws DimensionMismatch LAPACK.hesv!('U',A,rand(elty,11))
-            A = rand(elty,10,10)
-            A = A + A' #hermitian!
-            b = rand(elty,10)
-            c = A \ b
-            b,A = LAPACK.hesv_rook!('U',A,b)
-            @test b ≈ c
-            @test_throws DimensionMismatch LAPACK.hesv_rook!('U',A,rand(elty,11))
-        end
+        srand(935)
+        A = rand(elty,10,10)
+        A = A + A' #hermitian!
+        b = rand(elty,10)
+        c = A \ b
+        b,A = LAPACK.hesv!('U',A,b)
+        @test b ≈ c
+        @test_throws DimensionMismatch LAPACK.hesv!('U',A,rand(elty,11))
+        A = rand(elty,10,10)
+        A = A + A' #hermitian!
+        b = rand(elty,10)
+        c = A \ b
+        b,A = LAPACK.hesv_rook!('U',A,b)
+        @test b ≈ c
+        @test_throws DimensionMismatch LAPACK.hesv_rook!('U',A,rand(elty,11))
     end
 end
 

test/linalg/lu.jl

@@ -195,6 +195,7 @@ dimg  = randn(n)/2
 end
 
 @testset "conversion" begin
+    srand(3)
     a = Tridiagonal(rand(9),rand(10),rand(9))
     fa = Array(a)
     falu = lufact(fa)

test/linalg/tridiag.jl

@@ -5,8 +5,8 @@ using Main.TestHelpers.OAs
 
 using Base.Random: Sampler, SamplerRangeFast, SamplerRangeInt
 
-# Issue #6573
-guardsrand(0) do
+@testset "Issue #6573" begin
+    srand(0)
     rand()
     x = rand(384)
     @test find(x .== rand()) == []

test/random.jl

@@ -226,73 +226,73 @@ function randn_with_nans(n,p)
 end
 
 @testset "advanced sorting" begin
-    guardsrand(0xdeadbeef) do
-        for n in [0:10; 100; 101; 1000; 1001]
-            local r
-            r = -5:5
-            v = rand(r,n)
-            h = [sum(v .== x) for x in r]
-
-            for rev in [false,true]
-                # insertion sort (stable) as reference
-                pi = sortperm(v, alg=InsertionSort, rev=rev)
-                @test pi == sortperm(float(v), alg=InsertionSort, rev=rev)
-                @test isperm(pi)
-                si = v[pi]
-                @test [sum(si .== x) for x in r] == h
-                @test issorted(si, rev=rev)
-                @test all(issorted,[pi[si.==x] for x in r])
-                c = copy(v)
-                permute!(c, pi)
-                @test c == si
-                ipermute!(c, pi)
-                @test c == v
-
-                # stable algorithms
-                for alg in [MergeSort]
-                    p = sortperm(v, alg=alg, rev=rev)
-                    @test p == sortperm(float(v), alg=alg, rev=rev)
-                    @test p == pi
-                    s = copy(v)
-                    permute!(s, p)
-                    @test s == si
-                    ipermute!(s, p)
-                    @test s == v
-                end
-
-                # unstable algorithms
-                for alg in [QuickSort, PartialQuickSort(n)]
-                    p = sortperm(v, alg=alg, rev=rev)
-                    @test p == sortperm(float(v), alg=alg, rev=rev)
-                    @test isperm(p)
-                    @test v[p] == si
-                    s = copy(v)
-                    permute!(s, p)
-                    @test s == si
-                    ipermute!(s, p)
-                    @test s == v
-                end
+    srand(0xdeadbeef)
+    for n in [0:10; 100; 101; 1000; 1001]
+        local r
+        r = -5:5
+        v = rand(r,n)
+        h = [sum(v .== x) for x in r]
+
+        for rev in [false,true]
+            # insertion sort (stable) as reference
+            pi = sortperm(v, alg=InsertionSort, rev=rev)
+            @test pi == sortperm(float(v), alg=InsertionSort, rev=rev)
+            @test isperm(pi)
+            si = v[pi]
+            @test [sum(si .== x) for x in r] == h
+            @test issorted(si, rev=rev)
+            @test all(issorted,[pi[si.==x] for x in r])
+            c = copy(v)
+            permute!(c, pi)
+            @test c == si
+            ipermute!(c, pi)
+            @test c == v
+
+            # stable algorithms
+            for alg in [MergeSort]
+                p = sortperm(v, alg=alg, rev=rev)
+                @test p == sortperm(float(v), alg=alg, rev=rev)
+                @test p == pi
+                s = copy(v)
+                permute!(s, p)
+                @test s == si
+                ipermute!(s, p)
+                @test s == v
             end
 
-            v = randn_with_nans(n,0.1)
-            # TODO: alg = PartialQuickSort(n) fails here
-            for alg in [InsertionSort, QuickSort, MergeSort],
-                rev in [false,true]
-                # test float sorting with NaNs
-                s = sort(v, alg=alg, rev=rev)
-                @test issorted(s, rev=rev)
-                @test reinterpret(UInt64,v[isnan.(v)]) == reinterpret(UInt64,s[isnan.(s)])
-
-                # test float permutation with NaNs
+            # unstable algorithms
+            for alg in [QuickSort, PartialQuickSort(n)]
                 p = sortperm(v, alg=alg, rev=rev)
+                @test p == sortperm(float(v), alg=alg, rev=rev)
                 @test isperm(p)
-                vp = v[p]
-                @test isequal(vp,s)
-                @test reinterpret(UInt64,vp) == reinterpret(UInt64,s)
+                @test v[p] == si
+                s = copy(v)
+                permute!(s, p)
+                @test s == si
+                ipermute!(s, p)
+                @test s == v
             end
         end
+
+        v = randn_with_nans(n,0.1)
+        # TODO: alg = PartialQuickSort(n) fails here
+        for alg in [InsertionSort, QuickSort, MergeSort],
+            rev in [false,true]
+            # test float sorting with NaNs
+            s = sort(v, alg=alg, rev=rev)
+            @test issorted(s, rev=rev)
+            @test reinterpret(UInt64,v[isnan.(v)]) == reinterpret(UInt64,s[isnan.(s)])
+
+            # test float permutation with NaNs
+            p = sortperm(v, alg=alg, rev=rev)
+            @test isperm(p)
+            vp = v[p]
+            @test isequal(vp,s)
+            @test reinterpret(UInt64,vp) == reinterpret(UInt64,s)
+        end
     end
 end
+
 @testset "sortperm" begin
     inds = [
         1,1,1,2,2,2,3,3,3,4,4,4,5,5,5,6,6,6,7,7,7,8,8,8,9,9,9,10,

test/sorting.jl

@@ -1691,21 +1691,20 @@ end
 end
 
 @testset "sparse matrix normestinv" begin
-    guardsrand(1234) do
-        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)
-        if Base.USE_GPL_LIBS
-            # NOTE: normestinv is probabilistic, so must be included in the guardsrand block
-            @test Base.SparseArrays.normestinv(Ac,3) ≈ norm(inv(Array(Ac)),1) atol=1e-4
-            @test Base.SparseArrays.normestinv(Aci,3) ≈ norm(inv(Array(Aci)),1) atol=1e-4
-            @test Base.SparseArrays.normestinv(Ar) ≈ norm(inv(Array(Ar)),1) atol=1e-4
-            @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,0)
-            @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,21)
-        end
-        @test_throws DimensionMismatch Base.SparseArrays.normestinv(sprand(3,5,.9))
+    srand(1234)
+    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)
+    if Base.USE_GPL_LIBS
+        # NOTE: normestinv is probabilistic, so requires a fixed seed (set above in srand(1234))
+        @test Base.SparseArrays.normestinv(Ac,3) ≈ norm(inv(Array(Ac)),1) atol=1e-4
+        @test Base.SparseArrays.normestinv(Aci,3) ≈ norm(inv(Array(Aci)),1) atol=1e-4
+        @test Base.SparseArrays.normestinv(Ar) ≈ norm(inv(Array(Ar)),1) atol=1e-4
+        @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,0)
+        @test_throws ArgumentError Base.SparseArrays.normestinv(Ac,21)
     end
+    @test_throws DimensionMismatch Base.SparseArrays.normestinv(sprand(3,5,.9))
 end
 
 @testset "issue #13008" begin
@@ -1745,30 +1744,29 @@ end
 end
 
 @testset "factorization" begin
-    guardsrand(123) do
-        local A
-        A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2)
-        A = A + A'
-        @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A))))
-        A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2)
-        A = A*A'
-        @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A))))
-        A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2)
-        A = A + A.'
-        @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A))))
-        A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2)
-        A = A*A.'
-        @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A))))
-        @test factorize(triu(A)) == triu(A)
-        @test isa(factorize(triu(A)), UpperTriangular{Float64, SparseMatrixCSC{Float64, Int}})
-        @test factorize(tril(A)) == tril(A)
-        @test isa(factorize(tril(A)), LowerTriangular{Float64, SparseMatrixCSC{Float64, Int}})
-        @test !Base.USE_GPL_LIBS || factorize(A[:, 1:4])\ones(size(A, 1)) ≈ Array(A[:, 1:4])\ones(size(A, 1))
-        @test_throws ErrorException chol(A)
-        @test_throws ErrorException lu(A)
-        @test_throws ErrorException eig(A)
-        @test_throws ErrorException inv(A)
-    end
+    srand(123)
+    local A
+    A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2)
+    A = A + A'
+    @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A))))
+    A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2) + im*sprandn(5, 5, 0.2)
+    A = A*A'
+    @test !Base.USE_GPL_LIBS || abs(det(factorize(Hermitian(A)))) ≈ abs(det(factorize(Array(A))))
+    A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2)
+    A = A + A.'
+    @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A))))
+    A = sparse(Diagonal(rand(5))) + sprandn(5, 5, 0.2)
+    A = A*A.'
+    @test !Base.USE_GPL_LIBS || abs(det(factorize(Symmetric(A)))) ≈ abs(det(factorize(Array(A))))
+    @test factorize(triu(A)) == triu(A)
+    @test isa(factorize(triu(A)), UpperTriangular{Float64, SparseMatrixCSC{Float64, Int}})
+    @test factorize(tril(A)) == tril(A)
+    @test isa(factorize(tril(A)), LowerTriangular{Float64, SparseMatrixCSC{Float64, Int}})
+    @test !Base.USE_GPL_LIBS || factorize(A[:, 1:4])\ones(size(A, 1)) ≈ Array(A[:, 1:4])\ones(size(A, 1))
+    @test_throws ErrorException chol(A)
+    @test_throws ErrorException lu(A)
+    @test_throws ErrorException eig(A)
+    @test_throws ErrorException inv(A)
 end
 
 @testset "issue #13792, use sparse triangular solvers for sparse triangular solves" begin