!Example. pardiso_sym.f for Symmetric Linear Systems
!------------------------------------------------------------------
! Example program to show the use of the "PARDISO" routine
! for symmetric linear systems
!---------------------------------------------------------------------
! This program can be downloaded from the following site:
! http://www.pardiso-project.org
! (C) Olaf Schenk, Department of Computer Science,
! University of Basel, Switzerland.
! Email: olaf.schenk@unibas.ch
!
!---------------------------------------------------------------------
PROGRAM pardiso_sym
!INCLUDE 'mkl_pardiso.f90'
!!DEC$ OBJCOMMENT LIB:"gidpost_MT64.lib"
!DEC$ OBJCOMMENT LIB:"mkl_intel_lp64.lib"
!DEC$ OBJCOMMENT LIB:"mkl_intel_thread.lib"
!DEC$ OBJCOMMENT LIB:"mkl_core.lib"
!!#include "mkl_solver.h"

IMPLICIT NONE
!.. Internal solver memory pointer for 64-bit architectures
!.. INTEGER*8 pt(64)
!.. Internal solver memory pointer for 32-bit architectures
!.. INTEGER*4 pt(64)
!.. This is OK in both cases
INTEGER*8 pt(64)
!.. All other variables
INTEGER maxfct, mnum, mtype, phase, n, nrhs, error, msglvl
INTEGER iparm(64)
INTEGER ia(9)
INTEGER ja(18)
REAL*8 a(18)
REAL*8 b(8)
REAL*8 x(8)
INTEGER i, idum
REAL*8 waltime1, waltime2, ddum
!.. Fill all arrays containing matrix data.
DATA n /8/, nrhs /1/, maxfct /1/, mnum /1/
DATA ia /1,5,8,10,12,15,17,18,19/
DATA ja     /1,  3,    6,7,  &
    2,3,  5,         &
      3,        8,   &
        4,    7,     &
          5,6,7,     &
            6,  8,   &
              7,     &
                8/
DATA a  /7.d0,      1.d0,          2.d0,7.d0,            &
       -4.d0,8.d0,     2.d0,                      &
             1.d0,                    5.d0,       &
                  7.d0,     9.d0,                 &
                       5.d0,1.d0,5.d0,            &
                           -1.d0,     5.d0,       &
                                11.d0,            &
                                      5.d0/
integer omp_get_max_threads
external omp_get_max_threads
integer :: perm(8) !lmz0911

!..
!.. Set up PARDISO control parameter
!..
do i = 1, 64
   iparm(i) = 0
end do
iparm(1) = 1 ! no solver default
iparm(2) = 2 ! fill-in reordering from METIS
iparm(3) =6 !mkl_get_max_threads() !numbers of processors, value of MKL_NUM_THREADS
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) = 16 ! default logical fortran unit number for output
iparm(8) = 9 ! numbers of iterative refinement steps
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) = 0 ! not in use
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
iparm(23)=1   !lmz 0904 Inertia: number of negative eigenvalues.
iparm(30)=2
!iparm(31)=1  !lmz0910
iparm(36)=0  !lmz0910
error = 0 ! initialize error flag
msglvl = 0 ! don't print statistical information
mtype = -2 ! unsymmetric matrix symmetric, indefinite, no pivoting
!.. Initiliaze the internal solver memory pointer. This is only
! necessary for the FIRST call of the PARDISO solver.
  do i = 1, 64
     pt(i) = 0
  end do



!.. Reordering and Symbolic Factorization, This step also allocates
! all memory that is necessary for the factorization
phase = 11 ! only reordering and symbolic factorization
iparm(31)=3 !lmz0914

do i=1,8   !lmz0914
  perm(i)=0
end do
!perm(1)=1
!perm(7)=0 !lmz0914
perm(8)=1 !lmz0914

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,    &
    perm(8), nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Reordering completed ... '
IF (error .NE. 0) THEN
   WRITE(*,*) 'The following ERROR was detected: ', error
   STOP
END IF
WRITE(*,*) 'Number of nonzeros in factors = ',iparm(18)
WRITE(*,*) 'Number of factorization MFLOPS = ',iparm(19)
!.. Factorization.
phase = 22 ! only factorization
iparm(31)=3 !lmz0914

do i=1,8  !lmz0914
  perm(i)=0
end do
!perm(1)=1
!perm(7)=0 !lmz0914
perm(8)=1 !lmz0914
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,  &
     perm(8), nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Factorization completed ... '
IF (error .NE. 0) THEN
   WRITE(*,*) 'The following ERROR was detected: ', error
   STOP
ENDIF
!.. Back substitution and iterative refinement
iparm(8) = 0 ! max numbers of iterative refinement steps
phase = 33 ! only solve
do i = 1, n
   b(i) = 1.d0
end do

iparm(31)=3 !lmz0914

do i=1,8   !lmz0914
  perm(i)=0
end do
!perm(1)=1
!perm(7)=0 !lmz0914
perm(8)=1 !lmz0914

CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,  &
     perm(8), nrhs, iparm, msglvl, b, x, error)
WRITE(*,*) 'Solve completed ... '
WRITE(*,*) 'The solution of the system is '
DO i = 1, n
   WRITE(*,*) ' x(',i,') = ', x(i)
END DO
!.. Termination and release of memory
phase = -1 ! release internal memory
CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,  &
     perm, nrhs, iparm, msglvl, b, x, error)
END