Program test
   Implicit None
   Integer, Parameter :: N=13068, non_zeros=106245
   Integer :: error, iparm(64)
   Integer(8) :: pt(64)
   ! a, ja, ia: Matrix in crs format, a_full: Matrix in dense format
   Real(8), Allocatable :: a(:), a_full(:), rhs(:), sol1(:), sol2(:), r1(:), r2(:)
   Integer, Allocatable :: ja(:), ia(:), perm(:)
   Allocate(a(non_zeros), ja(non_zeros), ia(N+1), a_full(N**2), rhs(N), sol1(N), sol2(N), r1(N), r2(N), perm(N))
   ! Read matrix in crs format and right-hand side from file
   Write (*,*) "Reading matrix in crs format and right-hand side from file..."
   Open (1111, file="matrix.txt", form="unformatted", Status="old")
   Read (1111) a, ja, ia, rhs
   Close (1111)
   ! Convert matrix to dense format
   Write (*,*) "Converting matrix to dense format..."
   Call mkl_ddnscsr([1,1,1,2,0,0,0,0], N, N, a_full, N, a, ja, ia, error)
   If (error /= 0) Stop "Error!"
   ! Solve A*sol1 = rhs with DGESV
   Write (*,*) "Solving with DGESV..."
   sol1 = rhs
   Call dgesv(N,1,a_full,N,perm,sol1,N,error)
   If (error /= 0) Stop "Error!"
   ! Solve A*sol2 = rhs with pardiso: 
   Write (*,*) "Solving with Pardiso..."
   Call pardisoinit(pt, 11, iparm)
   iparm(27) = 1   ! Use matrix checker
   Call pardiso (pt, 1, 1, 11, 13, N, a, ia, ja, perm, 1, iparm, 0, rhs, sol2, error)
   If (error /= 0) Stop "Error!"
   ! Calculate r = A*sol-rhs
   Write (*,*) "Calculating residuum..."
   Call mkl_dcsrgemv('N', N, a, ia, ja, sol1, r1)
   r1 = r1-rhs
   Call mkl_dcsrgemv('N', N, a, ia, ja, sol2, r2)
   r2 = r2-rhs
   ! Output norm of residual
   Write (*,*) "Norm of residuum for DGESV:", Sum(Abs(r1))
   Write (*,*) "Norm of residuum for Pardiso:", Sum(Abs(r2))
End Program
