*   Modified version of MKL pardiso_unsym.f for problem given in MKL Forum
*       https://software.intel.com/en-us/forums/topic/588531
********************************************************************************
*   Copyright(C) 2004-2015 Intel Corporation. All Rights Reserved.
*
*   The source code, information  and  material ("Material") contained herein is
*   owned  by Intel Corporation or its suppliers or licensors, and title to such
*   Material remains  with Intel Corporation  or its suppliers or licensors. The
*   Material  contains proprietary information  of  Intel or  its  suppliers and
*   licensors. The  Material is protected by worldwide copyright laws and treaty
*   provisions. No  part  of  the  Material  may  be  used,  copied, reproduced,
*   modified, published, uploaded, posted, transmitted, distributed or disclosed
*   in any way  without Intel's  prior  express written  permission. No  license
*   under  any patent, copyright  or  other intellectual property rights  in the
*   Material  is  granted  to  or  conferred  upon  you,  either  expressly,  by
*   implication, inducement,  estoppel or  otherwise.  Any  license  under  such
*   intellectual  property  rights must  be express  and  approved  by  Intel in
*   writing.
*
*   *Third Party trademarks are the property of their respective owners.
*
*   Unless otherwise  agreed  by Intel  in writing, you may not remove  or alter
*   this  notice or  any other notice embedded  in Materials by Intel or Intel's
*   suppliers or licensors in any way.
*
********************************************************************************
*   Content : MKL PARDISO Fortran-77 example
*
********************************************************************************
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 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_unsym
        use mkl_pardiso
        IMPLICIT NONE
        !include 'mkl_pardiso.f77'
        INTEGER, PARAMETER :: NMAX=72408, NZMAX=1899415
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,nnz
        INTEGER iparm(64)
        INTEGER ia(NMAX+1)
        INTEGER ja(NZMAX),IR(NZMAX),JC(NZMAX)
        REAL*8 a(NZMAX),V(NZMAX)
        REAL*8 b(NMAX)
        REAL*8 x(NMAX)
        INTEGER i, idum(1), job(6), info, ibase
        REAL*8  ddum(1)
C.. Fill all arrays containing matrix data.
        DATA nrhs /1/, maxfct /1/, mnum /1/
        DATA job/2,1,0,-1,-1,3/
        open(11,file='matrix.txt',status='old')
        read(11,*)n
        read(11,*)ibase
        read(11,*)nnz
        read(11,*)(ir(i),jc(i),v(i),i=1,nnz)
        close(11)
        job(5)=nnz
        call mkl_dcsrcoo(job,n,a,ja,ia,nnz,v,ir,jc,info)
        write(*,*)' Info from dcsrcoo = ',info
        if(info.ne.0)stop
C..
C.. Setup PARDISO control parameter
C..
        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) = 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 components of x
        iparm(7) = 0 ! not in use
        iparm(8) = 9 ! numbers of iterative refinement steps
        iparm(9) = 0 ! not in use
        iparm(10) = 13 ! perturb 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
        iparm(35) = 0 ! zero-based indexing of ia,ja,perm
        error = 0 ! initialize error flag
        msglvl = 1 ! print statistical information
        mtype = 11 ! real unsymmetric
C.. Initialize the internal solver memory pointer. This is only
C necessary for the FIRST call of the PARDISO solver.
        DO i = 1, 64
            pt(i) = 0
        END DO
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
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     &  idum, nrhs, iparm, msglvl, ddum, ddum, error)
        WRITE(*,*) 'Reordering completed ... '
        IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP 1
        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,
     &  idum, nrhs, iparm, msglvl, ddum, ddum, error)
        WRITE(*,*) 'Factorization completed ... '
        IF (error .NE. 0) THEN
            WRITE(*,*) 'The following ERROR was detected: ', error
            STOP 1
        END IF
C.. Back substitution and iterative refinement
        iparm(8) = 2 ! max numbers of iterative refinement steps
        phase = 33 ! only factorization
        DO i = 1, n
            b(i) = 1.0D0
        END DO
        CALL pardiso (pt, maxfct, mnum, mtype, phase, n, a, ia, ja,
     &  idum, nrhs, iparm, msglvl, b, x, error)
        WRITE(*,*) 'Solve completed ... '
        WRITE(*,*) 'The solution of the system is '
        DO i = 1, 10
            WRITE(*,*) ' x(',i,') = ', x(i)
        END DO
        WRITE(*,*) '......'
        DO i = n-9, 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, idum, nrhs, iparm, msglvl, ddum, ddum, error)
      END PROGRAM pardiso_unsym
