Merge pull request #97 from mipals/fixing-macos-tests

Removing hard-coded MKL in tests

Project.toml

@@ -14,9 +14,10 @@ MKL_jll = "2021, 2022, 2023, 2024"
 julia = "1.6"
 
 [extras]
+Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "Random", "Printf"]
+test = ["Test", "Random", "Printf", "Pkg"]

README.md

@@ -19,7 +19,7 @@ library.
 ### MKL PARDISO
 
 By default Julia, will automatically install a suitable MKL for your platform.
-Note that if you use a mac you will need to pin `MKL_jll` to version 2022.
+Note that if you use a mac you will need to pin `MKL_jll` to version 2023.
 
 If you rather use a self installed MKL follow these instructions:
 

examples/exampleherm.jl

@@ -7,72 +7,72 @@ using SparseArrays
 using Random
 using Test
 
-# Script parameters.
-# -----------------
-verbose = false
-n       = 100
-lambda  = 3
-
-# Create the Hermitian positive definite matrix A and the vector b in the
-# linear system Ax = b.
-e = ones(n)
-e2 = ones(n-1)
-A = spdiagm(-1 => im*e2, 0 => lambda*e, 1 => -im*e2)
-b = rand(n) + im * zeros(n)
-
-# Initialize the PARDISO internal data structures.
-# ps = PardisoSolver()
-ps = MKLPardisoSolver()
-
-if verbose
-    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+function example_hermitian_psd(solver=MKLPardisoSolver)
+    # Script parameters.
+    # -----------------
+    verbose = false
+    n       = 100
+    lambda  = 3
+
+    # Create the Hermitian positive definite matrix A and the vector b in the
+    # linear system Ax = b.
+    e = ones(n)
+    e2 = ones(n-1)
+    A = spdiagm(-1 => im*e2, 0 => lambda*e, 1 => -im*e2)
+    b = rand(n) + im * zeros(n)
+
+    # Initialize the PARDISO internal data structures.
+    ps = solver()
+
+    if verbose
+        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    end
+
+    # If we want, we could just solve the system right now.
+    # Pardiso.jl will automatically detect the correct matrix type,
+    # solve the system and free the data
+    X1 = solve(ps, A, b)
+
+    # We also show how to do this in incremental steps.
+
+    ps = solver()
+
+    # First set the matrix type to handle general complex
+    # hermitian positive definite matrices
+    set_matrixtype!(ps, Pardiso.COMPLEX_HERM_POSDEF)
+
+    # Initialize the default settings with the current matrix type
+    pardisoinit(ps)
+
+    # Remember that we pass in a CSC matrix to Pardiso, so need
+    # to set the transpose iparm option.
+    fix_iparm!(ps, :N)
+
+    # Get the correct matrix to be sent into the pardiso function.
+    # :N for normal matrix, :T for transpose, :C for conjugate
+    A_pardiso = get_matrix(ps, A, :N)
+
+    # Analyze the matrix and compute a symbolic factorization.
+    set_phase!(ps, Pardiso.ANALYSIS)
+    pardiso(ps, A_pardiso, b)
+    @printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))
+
+    # Compute the numeric factorization.
+    set_phase!(ps, Pardiso.NUM_FACT)
+    pardiso(ps, A_pardiso, b)
+
+    # Compute the solutions X using the symbolic factorization.
+    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
+    x = similar(b) # Solution is stored in X
+    pardiso(ps, x, A_pardiso, b)
+    @printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))
+
+    # Compute the residuals.
+    r = abs.(A*x - b)
+    @printf("The maximum residual for the solution is %0.3g.\n",maximum(r))
+    @test norm(r) < 1e-10
+
+    # Free the PARDISO data structures.
+    set_phase!(ps, Pardiso.RELEASE_ALL)
+    pardiso(ps)
 end
-
-# If we want, we could just solve the system right now.
-# Pardiso.jl will automatically detect the correct matrix type,
-# solve the system and free the data
-X1 = solve(ps, A, b)
-
-# We also show how to do this in incremental steps.
-
-# ps = PardisoSolver()
-ps = MKLPardisoSolver()
-
-# First set the matrix type to handle general complex
-# hermitian positive definite matrices
-set_matrixtype!(ps, Pardiso.COMPLEX_HERM_POSDEF)
-
-# Initialize the default settings with the current matrix type
-pardisoinit(ps)
-
-# Remember that we pass in a CSC matrix to Pardiso, so need
-# to set the transpose iparm option.
-fix_iparm!(ps, :N)
-
-# Get the correct matrix to be sent into the pardiso function.
-# :N for normal matrix, :T for transpose, :C for conjugate
-A_pardiso = get_matrix(ps, A, :N)
-
-# Analyze the matrix and compute a symbolic factorization.
-set_phase!(ps, Pardiso.ANALYSIS)
-pardiso(ps, A_pardiso, b)
-@printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))
-
-# Compute the numeric factorization.
-set_phase!(ps, Pardiso.NUM_FACT)
-pardiso(ps, A_pardiso, b)
-
-# Compute the solutions X using the symbolic factorization.
-set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
-x = similar(b) # Solution is stored in X
-pardiso(ps, x, A_pardiso, b)
-@printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))
-
-# Compute the residuals.
-r = abs.(A*x - b)
-@printf("The maximum residual for the solution is %0.3g.\n",maximum(r))
-@test norm(r) < 1e-10
-
-# Free the PARDISO data structures.
-set_phase!(ps, Pardiso.RELEASE_ALL)
-pardiso(ps)

examples/examplesym.jl

@@ -12,75 +12,76 @@ using Random
 using Printf
 using Test
 
-# Script parameters.
-# -----------------
-verbose = false
-n       = 4  # The number of equations.
-m       = 3  # The number of right-hand sides.
-
-A = sparse([ 1. 0 -2  3
-             0  5  1  2
-            -2  1  4 -7
-             3  2 -7  5 ])
-
-# Generate a random collection of right-hand sides.
-B = rand(n,m)
-
-# Initialize the PARDISO internal data structures.
-# ps = PardisoSolver()
-ps = MKLPardisoSolver()
-
-if verbose
-    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+function example_symmetric(solver=MKLPardisoSolver)
+
+    # Script parameters.
+    # -----------------
+    verbose = false
+    n       = 4  # The number of equations.
+    m       = 3  # The number of right-hand sides.
+
+    A = sparse([ 1. 0 -2  3
+                0  5  1  2
+                -2  1  4 -7
+                3  2 -7  5 ])
+
+    # Generate a random collection of right-hand sides.
+    B = rand(n,m)
+
+    # Initialize the PARDISO internal data structures.
+    ps = solver()
+
+    if verbose
+        set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
+    end
+
+    # If we want, we could just solve the system right now.
+    # Pardiso.jl will automatically detect the correct matrix type,
+    # solve the system and free the data
+    X1 = solve(ps, A, B)
+
+    # We also show how to do this in incremental steps.
+
+    ps = solver()
+
+    # First set the matrix type to handle general real symmetric matrices
+    set_matrixtype!(ps, Pardiso.REAL_SYM_INDEF)
+
+    # Initialize the default settings with the current matrix type
+    pardisoinit(ps)
+
+    # Remember that we pass in a CSC matrix to Pardiso, so need
+    # to set the transpose iparm option.
+    fix_iparm!(ps, :N)
+
+    # Get the correct matrix to be sent into the pardiso function.
+    # :N for normal matrix, :T for transpose, :C for conjugate
+    A_pardiso = get_matrix(ps, A, :N)
+
+    # Analyze the matrix and compute a symbolic factorization.
+    set_phase!(ps, Pardiso.ANALYSIS)
+    set_perm!(ps, randperm(n))
+    pardiso(ps, A_pardiso, B)
+    @printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))
+
+    # Compute the numeric factorization.
+    set_phase!(ps, Pardiso.NUM_FACT)
+    pardiso(ps, A_pardiso, B)
+    @printf("The matrix has %d positive and %d negative eigenvalues.\n",
+            get_iparm(ps, 22), get_iparm(ps, 23))
+
+    # Compute the solutions X using the symbolic factorization.
+    set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
+    X = similar(B) # Solution is stored in X
+    pardiso(ps, X, A_pardiso, B)
+    @printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))
+
+    # Compute the residuals.
+    R = maximum(abs.(A*X - B))
+    @printf("The maximum residual for the solution X is %0.3g.\n", R)
+    @test R < 1e-10
+
+    # Free the PARDISO data structures.
+    set_phase!(ps, Pardiso.RELEASE_ALL)
+    pardiso(ps)
 end
-
-# If we want, we could just solve the system right now.
-# Pardiso.jl will automatically detect the correct matrix type,
-# solve the system and free the data
-X1 = solve(ps, A, B)
-
-# We also show how to do this in incremental steps.
-
-# ps = PardisoSolver()
-ps = MKLPardisoSolver()
-
-# First set the matrix type to handle general real symmetric matrices
-set_matrixtype!(ps, Pardiso.REAL_SYM_INDEF)
-
-# Initialize the default settings with the current matrix type
-pardisoinit(ps)
-
-# Remember that we pass in a CSC matrix to Pardiso, so need
-# to set the transpose iparm option.
-fix_iparm!(ps, :N)
-
-# Get the correct matrix to be sent into the pardiso function.
-# :N for normal matrix, :T for transpose, :C for conjugate
-A_pardiso = get_matrix(ps, A, :N)
-
-# Analyze the matrix and compute a symbolic factorization.
-set_phase!(ps, Pardiso.ANALYSIS)
-set_perm!(ps, randperm(n))
-pardiso(ps, A_pardiso, B)
-@printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))
-
-# Compute the numeric factorization.
-set_phase!(ps, Pardiso.NUM_FACT)
-pardiso(ps, A_pardiso, B)
-@printf("The matrix has %d positive and %d negative eigenvalues.\n",
-        get_iparm(ps, 22), get_iparm(ps, 23))
-
-# Compute the solutions X using the symbolic factorization.
-set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
-X = similar(B) # Solution is stored in X
-pardiso(ps, X, A_pardiso, B)
-@printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))
-
-# Compute the residuals.
-R = maximum(abs.(A*X - B))
-@printf("The maximum residual for the solution X is %0.3g.\n", R)
-@test R < 1e-10
-
-# Free the PARDISO data structures.
-set_phase!(ps, Pardiso.RELEASE_ALL)
-pardiso(ps)