C************************************************************************************************
C Copyright (C) 2022 until present, by Panua Technologies Sagl, Switzerland. All Rights Reserved.                      
C 
C Example program to show the use of the "PARDISO" routine for symmetric linear systems
C
C This program can be downloaded from the following site: http://www.panua.ch 
C The above copyright and permission notice must appear in all partial and complete
C copies of this file.               
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  dparm(64) 
        REAL*8  a(18) 
        REAL*8  b(8)
        REAL*8  x(8)

        INTEGER i, idum, solver
        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  .. set right hand side
      do i = 1, n
        b(i) = 1.d0
      end do

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.                                     
     
      mtype     = -2  ! unsymmetric matrix symmetric, indefinite
      solver    =  1  ! use sparse direct method

C  .. pardiso_chk_matrix(...)
C     Checks the consistency of the given matrix.
C     Use this functionality only for debugging purposes
      CALL pardiso_chkmatrix  (mtype, n, a, ia, ja, error);
      IF (error .NE. 0) THEN
         WRITE(*,*) 'The following ERROR was detected: ', error
        STOP
      ENDIF

C ..  pardiso_chkvec(...)
C     Checks the given vectors for infinite and NaN values
C     Input parameters (see PARDISO user manual for a description):
C     Use this functionality only for debugging purposes
      CALL pardiso_chkvec (n, nrhs, b, error);
      IF (error .NE. 0) THEN
         WRITE(*,*) 'The following ERROR was detected: ', error
        STOP
      ENDIF
 
C ..  pardiso_printstats(...) 
C     prints information on the matrix to STDOUT.
C     Use this functionality only for debugging purposes

      CALL pardiso_printstats (mtype, n, a, ia, ja, nrhs, b, error);
      IF (error .NE. 0) THEN
         WRITE(*,*) 'The following ERROR was detected: ', error
        STOP
       ENDIF
      
C  .. Pardiso license check and initialize solver
      call pardisoinit(pt, mtype, solver, iparm, dparm, error)
      IF (error .NE. 0) THEN
        IF (error.EQ.-10 ) WRITE(*,*) 'No license file found'
        IF (error.EQ.-11 ) WRITE(*,*) 'License is expired'
        IF (error.EQ.-12 ) WRITE(*,*) 'Wrong username or hostname'
        STOP
      ELSE
        WRITE(*,*) '[PARDISO]: License check was successful ... '
      END IF
     
C  .. Numbers of Processors ( value of OMP_NUM_THREADS )
      iparm(3) = 1

      phase     = 12  ! set up of the preconditioner
      msglvl    = 1

      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, 
     1              idum, nrhs, iparm, msglvl, ddum, ddum, error, dparm) 

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

C  .. Solve with SQMR
      phase     = 33  ! only solve

C  ..Set initial solution to zero
      DO i = 1, n
        x(i) = 0
      END DO

      CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja, 
     1              idum, nrhs, iparm, msglvl, b, x, error, dparm) 

      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, dparm) 
      END
