C----------------------------------------------------------------------
C       Example program to show the use of the "PARDISO" routine
C       for symmetric linear systems
C---------------------------------------------------------------------
C       This program can be downloaded from the following site:
C       http://www.pardiso-project.org 
C
C   (C) Olaf Schenk, Department of Computer Science,
C       University of Basel, Switzerland.
C       Email: olaf.schenk@unibas.ch
C
C---------------------------------------------------------------------
        PROGRAM pardiso_sym
        IMPLICIT NONE

C..     Internal solver memory pointer for 64-bit architectures
C..     INTEGER*8 pt(64)
C..     Internal solver memory pointer for 32-bit architectures
C..     INTEGER*4 pt(64)
C..     This is OK in both cases
        INTEGER*8 pt(64)

C..     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

C.. 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        /1,          3,                 6,    7,
     2               2,    3,          5,                   
     3                     3,                             8,            
     4                          4,                  7,     
     5                                 5,     6,    7,
     6                                        6,          8,         
     7                                              7,      
     8                                                    8/
      
        DATA a
     1     /7.d0,       1.d0,              2.d0, 7.d0,
     2           -4.d0, 8.d0,        2.d0,                   
     3                  1.d0,                          5.d0,            
     4                        7.d0,              9.d0,     
     5                               5.d0, 1.d0, 5.d0,
     6                                     0.d0,       5.d0,         
     7                                           11.d0,      
     8                                                 5.d0/
       

C
C  .. Setup Pardiso control parameters und initialize the solvers     
C     internal adress pointers. This is only necessary for the FIRST   
C     call of the PARDISO solver.                                     
C     
      mtype     = -2  ! unsymmetric matrix symmetric, indefinite
      call pardisoinit(pt, mtype, iparm)
     
C  .. Numbers of Processors ( value of OMP_NUM_THREADS )
      iparm(3) = 1

C..   Reordering and Symbolic Factorization, This step also allocates
C     all memory that is necessary for the factorization
      
      phase     = 11     ! only reordering and symbolic factorization
      msglvl    = 1       ! with statistical information
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     1              idum, nrhs, iparm, msglvl, ddum, ddum, 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)

C.. Factorization.
      phase     = 22  ! only factorization
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, 
     1              idum, nrhs, iparm, msglvl, ddum, ddum, error) 

      WRITE(*,*) 'Factorization completed ... '
      IF (error .NE. 0) THEN
         WRITE(*,*) 'The following ERROR was detected: ', error
        STOP
      ENDIF 

C.. Back substitution and iterative refinement
      iparm(8)  = 1   ! max numbers of iterative refinement steps
      phase     = 33  ! only factorization
      do i = 1, n
        b(i) = 1.d0
      end do

      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, 
     1              idum, 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

C.. Termination and release of memory
      phase     = -1           ! release internal memory
      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, ddum, idum, idum, 
     1              idum, nrhs, iparm, msglvl, ddum, ddum, error) 
      END