fix nprocs acting weird in MKL when libpardiso is loaded (#67)

src/Pardiso.jl

@@ -123,9 +123,13 @@ function __init__()
         @warn "MKLROOT not set, MKL Pardiso solver will not be functional"
     end
 
+    # This is apparently needed for MKL to not get stuck on 1 thread when
+    # libpardiso is loaded in the block below...
+    get_nprocs_mkl()
+
     if PARDISO_LIB_FOUND
         try
-            libpardiso = Libdl.dlopen(PARDISO_PATH, Libdl.RTLD_GLOBAL)
+            libpardiso = Libdl.dlopen(PARDISO_PATH)
             init[] = Libdl.dlsym(libpardiso, "pardisoinit")
             pardiso_f[] = Libdl.dlsym(libpardiso, "pardiso")
             pardiso_chkmatrix[] = Libdl.dlsym(libpardiso, "pardiso_chkmatrix")

src/mkl_pardiso.jl

@@ -34,9 +34,11 @@ show(io::IO, ps::MKLPardisoSolver) = print(io, string("$MKLPardisoSolver:\n",
                                   "\tMatrix type: $(MATRIX_STRING[get_matrixtype(ps)])\n",
                                   "\tPhase: $(PHASE_STRING[get_phase(ps)])"))
 
-set_nprocs!(ps::MKLPardisoSolver, n::Integer) =
+set_nprocs!(ps::MKLPardisoSolver, n::Integer) = set_nprocs_mkl!(n)
+set_nprocs_mkl!(n::Integer) =
     ccall((:mkl_domain_set_num_threads, libmkl_rt), Cvoid, (Ptr{Int32}, Ptr{Int32}), Ref((Int32(n))), Ref(MKL_DOMAIN_PARDISO))
-get_nprocs(ps::MKLPardisoSolver) =
+get_nprocs(ps::MKLPardisoSolver) = get_nprocs_mkl()
+get_nprocs_mkl() =
     ccall((:mkl_domain_get_max_threads, libmkl_rt), Int32, (Ptr{Int32},), Ref(MKL_DOMAIN_PARDISO))
 
 valid_phases(ps::MKLPardisoSolver) = keys(MKL_PHASES)

test/runtests.jl

@@ -82,40 +82,39 @@ if Pardiso.PARDISO_LOADED[]
 
     # try some random matrices
     m = 50; n = 15; p = .1
-    for pardiso_type in psolvers
-        for T in (Float64, )#ComplexF64)
-            ps = pardiso_type()
-            pardisoinit(ps)
-            if T == Float64
-                set_matrixtype!(ps, 11)
-            else
-                set_matrixtype!(ps, 13)
-            end
-            for j ∈ 1:100
-                A = 5I + sprand(T,m,m,p)
-                A⁻¹ = inv(Matrix(A))
-                B = sprand(T,m,n,p)
-                C = sprand(T,n,m,p)
-                D = 5I + sprand(T,n,n,p)
-                M = [A B; C D]
-
-                # test integer block specification
-                S = schur_complement(ps, M, n);
-                @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
-
-                # test sparse vector block specification
-                x = spzeros(T,m+n)
-                x[(m+1):(m+n)] .= 1
-                S = schur_complement(ps, M, x);
-                @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
-
-                # test sparse matrix block specification
-                x = spzeros(T,m+n,2)
-                x[(m+1):(m+n-1),1] .= 1
-                x[end,2] = 1
-                S = schur_complement(ps, M, x);
-                @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
-            end
+    ps = PardisoSolver()
+    for T in (Float64, )#ComplexF64)
+        ps = PardisoSolver()
+        pardisoinit(ps)
+        if T == Float64
+            set_matrixtype!(ps, 11)
+        else
+            set_matrixtype!(ps, 13)
+        end
+        for j ∈ 1:100
+            A = 5I + sprand(T,m,m,p)
+            A⁻¹ = inv(Matrix(A))
+            B = sprand(T,m,n,p)
+            C = sprand(T,n,m,p)
+            D = 5I + sprand(T,n,n,p)
+            M = [A B; C D]
+
+            # test integer block specification
+            S = schur_complement(ps, M, n);
+            @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
+
+            # test sparse vector block specification
+            x = spzeros(T,m+n)
+            x[(m+1):(m+n)] .= 1
+            S = schur_complement(ps, M, x);
+            @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
+
+            # test sparse matrix block specification
+            x = spzeros(T,m+n,2)
+            x[(m+1):(m+n-1),1] .= 1
+            x[end,2] = 1
+            S = schur_complement(ps, M, x);
+            @test norm(D - C*A⁻¹*B - S) < 1e-10*(m+n)^2
         end
     end
 end # testset