program test

include 'mkl_spblas.fi'

character(LEN=19):: symm
integer:: n=4,nz=7
integer:: jrn(7), jcn(7)
doubleprecision:: matA_coo(7)

doubleprecision::vB(4),vX(4)

! PARDISO parameters and settings
integer nrhs, maxfct,m_type, mnum, phase, error, msglvl, idum
integer iparm(64)
integer*8 pt(64)

! Row index and colum pointer - CSR format
integer,allocatable::row_ptr(:),col_ind(:)
! Sparse matrix data - CSR format
double precision,allocatable::matA_csr(:)

! For conversion utility
integer info, job(8)

! Timing
double precision tm1,tm2,tm3,tm_sec

! Sample symmetric matrix - data in upper half
symm='symmetric'
jrn = (/1,2,3,4,1,2,3/)
jcn = (/1,2,3,4,2,3,4/)
matA_coo = (/1.0d0,2.0d0,3.0d0,4.0d0,5.0d0,6.0d0,7.0d0/)

! From Intel example - use this as default settings
! Note: use this as a consistent set of parameters for comparison to 
! GPU solver since the latter does not require this kind of fine-tuning.
iparm = 0
iparm(1) = 1 ! no solver default
iparm(2) = 3 ! fill-in reordering from METIS (parallel version)
iparm(3) = 1 ! numbers of processors
iparm(4) = 0 ! no iterative-direct algorithm
iparm(5) = 0 ! no user fill-in reducing permutation
iparm(6) = 0 ! =0 solution on the first n compoments of x
iparm(7) = 0 ! not in use
iparm(8) = 10 ! numbers of iterative refinement steps: Modified - same as for GPU
iparm(9) = 0 ! not in use
iparm(10) = 13 ! perturbe the pivot elements with 1E-13
iparm(11) = 1 ! use nonsymmetric permutation and scaling MPS
iparm(12) = 0 ! not in use
iparm(13) = 1 ! maximum weighted matching algorithm is switched-on (default for non-symmetric)
iparm(14) = 0 ! Output: number of perturbed pivots
iparm(15) = 0 ! not in use
iparm(16) = 0 ! not in use
iparm(17) = 0 ! not in use
iparm(18) = -1 ! Output: number of nonzeros in the factor LU
iparm(19) = -1 ! Output: Mflops for LU factorization
iparm(20) = 0 ! Output: Numbers of CG Iterations
error = 0 ! initialize error flag
msglvl = 0 ! print statistical information
nrhs=1 ! # of right-hand sides to solve
maxfct=1 ! Maximal number of factors in memory
mnum=1 ! The number of matrix (from 1 to maxfct) to solve

if (symm.eq.'symmetric') then
  m_type = 1 ! real symmetric
else
  m_type = 11 ! real unsymmetric
endif  

! Local storage of sparse matrix
allocate(col_ind(nz),row_ptr(n+1),matA_csr(nz))

! Convert triplet/COO data to CSR format
! Note: this requires Intel MKL 10.3 or above as 
! sort option was not available in previous versions.
job(1)=2 ! convert COO to CSR and sort
job(2)=1 ! 1-based indexing in CSR
job(3)=1 ! 1-based indexing in COO
job(5)=nz ! # of non-zeros in matrix A
job(6)=0 ! Fill all CSR arrays
call mkl_dcsrcoo(job, n, matA_csr, col_ind, row_ptr, nz, matA_coo, jrn, jcn, info)

write(6,*) n,nz
write(6,*) col_ind
write(6,*) row_ptr
write(6,*) matA_csr

! Init internal solver memory 
pt = 0

phase=11
tm_sec=SECNDS(0.0)
call pardiso(pt, maxfct, mnum, m_type,phase, n, matA_csr, row_ptr, col_ind, idum, nrhs, iparm, msglvl, vB, vX, error)
tm1=SECNDS(tm_sec)
if (error .ne. 0) then
  write(6,*) 'The following PARDISO error was detected: ', error
  stop
endif
write(6,*) '--->'
write(6,*) 'Time for symbolic analysis (sec)'
write(6,*) tm1

phase=22
tm_sec=SECNDS(0.0)
call pardiso(pt, maxfct, mnum, m_type,phase, n, matA_csr, row_ptr, col_ind, idum, nrhs, iparm, msglvl, vB, vX, error)
tm2=SECNDS(tm_sec)
if (error .ne. 0) then
  write(6,*) 'The following PARDISO error was detected: ', error
  stop
endif
write(6,*) '--->'
write(6,*) 'Time for numeric factorization (sec)'
write(6,*) tm2

phase=33
tm_sec=SECNDS(0.0)
call pardiso(pt, maxfct, mnum, m_type,phase, n, matA_csr, row_ptr, col_ind, idum, nrhs, iparm, msglvl, vB, vX, error)
tm3=SECNDS(tm_sec)
if (error .ne. 0) then
  write(6,*) 'The following PARDISO error was detected: ', error
  stop
endif
write(6,*) '--->'
write(6,*) 'Time for backsolve+iterative refinement (sec)'
write(6,*) tm3

! Clean-up
phase=-1
call pardiso(pt, maxfct, mnum, m_type,phase, n, matA_csr, row_ptr, col_ind, idum, nrhs, iparm, msglvl, vB, vX, error)
deallocate(col_ind,row_ptr,matA_csr)

endprogram