# This is an example script demonstrating how PARDISO works on a
# medium-sized Hermitian positive definite matrix.
using Pardiso
# 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,1)
e2 = ones(n-1,1)
A = spdiagm((im*e2, lambda*e, -im*e2), (-1,0,1))
b = rand(n,1) + im * zeros(n,1)
# Initialize the PARDISO internal data structures.
ps = PardisoSolver()
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)
# 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)
# 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))
# Free the PARDISO data structures.
set_phase!(ps, Pardiso.RELEASE_ALL)
pardiso(ps)