! Tests functionality of the PARDISO Solver (in MKL)
! On multicore machines set the OMP_NUM_THREADS variable to 2 or more
! 
! stali@purdue.edu

program main
  implicit none

  integer(8) :: pt(64)
  integer :: n, nnz, nrhs, maxfct, mnum, mtype, phase, perm, error, msglvl
  integer :: iparm(64)
  real(8), allocatable :: a(:), b(:), x(:)
  integer, allocatable :: ja(:), ia(:)

  n=5 ; nnz=13; nrhs=1 ; maxfct=1 ; mnum=1
  mtype = 1  ! matrix is real and symmetric

  allocate(a(nnz),ja(nnz),ia(n+1),b(n),x(n))

  a  = (/2.0, -1.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, -1.0, 2.0, -1.0, &
       -1.0, 2.0/)
  ja = (/1, 2, 1, 2, 3, 2, 3, 4, 3, 4, 5, 4, 5/)
  ia = (/1, 3, 6, 9, 12, 14/)
  b  = (/0.0, 0.0, 0.0, 0.0, 100.0/)

  iparm    = 0   ! explicitly set all solver parameter values to zero
  iparm(1) = 0   ! use solver defaults for iparm(2) and iparm(4) through iparm(64)
  msglvl   = 0   ! dont print any statistical information
  pt       = 0   ! explicitly initialize all pt values to zero

  phase    = 13  ! analyze, factorize and solve
  call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
       & perm, nrhs, iparm, msglvl, b, x, error)
  print*, x

  phase    = -1  ! release internal memory
  call pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, &
       & perm, nrhs, iparm, msglvl, b, x, error)

end program main