module minput

USE MKL95_PRECISION, ONLY: WP => DP
USE MKL95_LAPACK, ONLY: GBSV


implicit none
save

    integer,parameter :: dp = KIND(1.0D0)

    integer,parameter :: N = 8                                                        !4 - 8 - 144
    integer,parameter :: SIZE = 128
    
    integer,dimension(N+1) :: IA
    integer,dimension(18) :: JA                                                       !12 - 18 - 488
    real(dp),dimension(18) :: A, BILU0                                                !12 - 18 - 488      
    real(dp),dimension(N) :: TRVEC
    
    integer,dimension(SIZE) :: IPAR
    integer :: IERR
    real(dp),dimension(SIZE) :: DPAR
    real(dp),dimension(N*(2*min(150,N)+1)+(min(150,N)*(min(150,N)+9))/2+1) :: TMP
    real(dp),dimension(N) :: EXPECTED_SOLUTION, RHS, B, COMPUTED_SOLUTION, RESIDUAL
    
    integer,parameter :: MATSIZE = 12
    integer,parameter :: INCX = 1
    integer,parameter :: REF_NIT = 2
    real(dp),parameter :: REF_NORM2 = 7.772387d0
    real(dp) :: NRM2
    
    integer :: ITERCOUNT, RCI_REQUEST, I
    real(dp) :: DVAR
    
    real(dp),external :: DNRM2
    
      

    contains
    
    subroutine subrinput
    implicit none
    
!        IA = (/ 1,4,7,10,13 /)
!        IA = (/ 1, 4, 7, 10, 12, 14, 16, 18, 19 /)
        IA = (/ 1, 5, 8, 10, 12, 15, 17, 18, 19 /)  !from Pardiso example
        
!        open(104,file='innzele.txt')                      
!            read(104,*) IA
!        close(104)
        
!        JA = (/ 1,2,3,1,2,4,1,3,4,2,3,4 /)
!        JA = (/ 1, 2, 5, 2, 3, 6, 3, 4, 7, 4, 8, 5, 6, 6, 7, 7, 8, 8 /)

        JA = (/ 1,    3,       6, 7,    &
                   2, 3,    5,          &
                      3,             8, &
                         4,       7,    &
                            5, 6, 7,    &
                               6,    8, &
                                  7,    &
                                     8 /)

!        open(101,file='colnzele.txt')                      
!            read(101,*) JA
!        close(101)
        
!        A = (/ 4.,-1.,-1.,-1.,4.,-1.,-1.,4.,-1.,-1.,-1.,4. /)
!        A = (/ 5.48529E-2_dp, -3.664746E-8_dp, -1.30303E-11_dp, 0.1097058_dp, -3.664746E-8_dp, -2.60606E-11_dp, 0.1097058_dp, -3.664746E-8_dp, -2.60606E-11_dp, 5.48529E-2_dp, -1.30303E-11_dp, 0.740804885_dp, -4.949348E-7_dp, 1.4816097709_dp, -4.949348E-7_dp, 1.4816097709_dp, -4.949348E-7_dp, 0.740804885_dp /)
    
        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 /)    
    
!        open(102,file='nzele.txt')                      
!            read(102,*) A
!        close(102)
    
!        EXPECTED_SOLUTION = (/ 1.0,1.0,1.0,1.0 /)
!        EXPECTED_SOLUTION = (/ 1890986.34329, 3781972.68659, 3781972.68659, 1890986.34329, 25538332.3725, 51076664.7450, 51076664.7450, 25538332.3725 /)
    
        !---------------------------------------------------------------------------
        ! Initialize variables and the right hand side through matrix-vector product
        !---------------------------------------------------------------------------
!        CALL MKL_DCSRGEMV('N', N, A, IA, JA, EXPECTED_SOLUTION, RHS)

!        RHS = (/ 1890986.34329, 3781972.68659, 3781972.68659, 1890986.34329, 25538332.3725, 51076664.7450, 51076664.7450, 25538332.3725 /)

        do i = 1, N
           RHS(i) = 1.d0
        end do

!        open(103,file='RHSnn.txt')                      
!            read(103,*) RHS
!        close(103)
        
        !---------------------------------------------------------------------------
        ! Save the right-hand side in vector B for future use
        !---------------------------------------------------------------------------
	    CALL DCOPY(N, RHS, 1, B, 1)
        !---------------------------------------------------------------------------
        ! Initialize the initial guess
        !---------------------------------------------------------------------------
        DO I=1,N
            COMPUTED_SOLUTION(I)=0.d0
        END DO
!        COMPUTED_SOLUTION(1)=100.d0
        
!        DO I=1,N
!            COMPUTED_SOLUTION(I)=0.d0    !Initial Pressure as initial guess  !34473786
!        END DO
        
!        open(103,file='PardisoSol.txt')                      
!            read(103,*) COMPUTED_SOLUTION   !THIS IS THE ACTUAL TRUE SOLUTION!!!!
!        close(103)
        
!        print *, COMPUTED_SOLUTION/6894.7572931684
        
        !---------------------------------------------------------------------------
        ! Initialize the solver
        !---------------------------------------------------------------------------
        CALL DFGMRES_INIT(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)
        IF (RCI_REQUEST.NE.0) GOTO 999

        !---------------------------------------------------------------------------
        ! Calculate ILU0 preconditioner.
        !                      !ATTENTION!
        ! DCSRILU0 routine uses some IPAR, DPAR set by DFGMRES_INIT routine.
        ! Important for DCSRILU0 default entries set by DFGMRES_INIT are
        ! ipar(2) = 6 - output of error messages to the screen,
        ! ipar(6) = 1 - allow output of error messages,
        ! ipar(31)= 0 - abort DCSRILU0 calculations if routine meets zero diagonal element.
        !
        ! If ILU0 is going to be used out of MKL FGMRES context, than the values
        ! of ipar(2), ipar(6), ipar(31), dpar(31), and dpar(32) should be user
        ! provided before the DCSRILU0 routine call.
        !
        ! In this example, specific for DCSRILU0 entries are set in turn:
        ! ipar(31)= 1 - change small diagonal value to that given by dpar(32),
        ! dpar(31)= 1.D-20 instead of the default value set by DFGMRES_INIT.
        !                  It is a small value to compare a diagonal entry with it.
        ! dpar(32)= 1.D-16 instead of the default value set by DFGMRES_INIT.
        !                  It is the target value of the diagonal value if it is
        !                  small as compared to dpar(31) and the routine should change
        !                  it rather than abort DCSRILU0 calculations.
        !---------------------------------------------------------------------------

        IPAR(31)=1
        DPAR(31)=1.D-20
        DPAR(32)=1.D-16
        CALL DCSRILU0(N, A, IA, JA, BILU0, IPAR, DPAR, IERR)
!        NRM2=DNRM2(MATSIZE, BILU0, INCX)

        IF(IERR.ne.0) THEN
            PRINT *, 'Error after calculation of the',' preconditioner DCSRILU0',IERR
            PRINT *, 'Unfortunately, FGMRES+ILU0 Fortran example',' has FAILED',' during the preconditioning of the matrix'
        END IF

        !---------------------------------------------------------------------------
        ! Set the desired parameters:
        ! do the restart after 2 iterations
        ! LOGICAL parameters:
        ! do not do the stopping test for the maximal number of iterations
        ! do the Preconditioned iterations of FGMRES method
        ! Set parameter IPAR(11) for preconditioner call. For this example,
        ! it reduces the number of iterations.
        ! DOUBLE PRECISION parameters
        ! set the relative tolerance to 1.0D-3 instead of default value 1.0D-6
        ! NOTE. Preconditioner may increase the number of iterations for an
        ! arbitrary case of the system and initial guess and even ruin the
        ! convergence. It is user's responsibility to use a suitable preconditioner
        ! and to apply it skillfully.
        !---------------------------------------------------------------------------
        IPAR(15)=2
        IPAR(8)=0
        IPAR(11)=1
        DPAR(1)=1.0D-3
        !---------------------------------------------------------------------------
        ! Check the correctness and consistency of the newly set parameters
        !---------------------------------------------------------------------------
        CALL DFGMRES_CHECK(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)
        IF (RCI_REQUEST.NE.0) GOTO 999
        !---------------------------------------------------------------------------
        ! Print the info about the RCI FGMRES method
        !---------------------------------------------------------------------------
        
        !NOT PRINTING IT HERE - LOOK AT ORIGINAL EXAMPLE FOR MORE DETAILS
        
        !---------------------------------------------------------------------------
        ! Compute the solution by RCI (P)FGMRES solver with preconditioning
        ! Reverse Communication starts here
        !---------------------------------------------------------------------------
        1     CALL DFGMRES(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP)
        !---------------------------------------------------------------------------
        ! If RCI_REQUEST=0, then the solution was found with the required precision
        !---------------------------------------------------------------------------
        IF (RCI_REQUEST.EQ.0) GOTO 3
        !---------------------------------------------------------------------------
        ! If RCI_REQUEST=1, then compute the vector A*TMP(IPAR(22))
        ! and put the result in vector TMP(IPAR(23))
        !---------------------------------------------------------------------------
        IF (RCI_REQUEST.EQ.1) THEN
      	    CALL MKL_DCSRGEMV('N',N, A, IA, JA, TMP(IPAR(22)), TMP(IPAR(23)))
      	    GOTO 1
        END IF
        !---------------------------------------------------------------------------
        ! If RCI_request=2, then do the user-defined stopping test
        ! The residual stopping test for the computed solution is performed here
        !---------------------------------------------------------------------------
        ! NOTE: from this point vector B(N) is no longer containing the right-hand
        ! side of the problem! It contains the current FGMRES approximation to the
        ! solution. If you need to keep the right-hand side, save it in some other
        ! vector before the call to DFGMRES routine. Here we saved it in vector
        ! RHS(N). The vector B is used instead of RHS to preserve the original
        ! right-hand side of the problem and guarantee the proper restart of FGMRES
        ! method. Vector B will be altered when computing the residual stopping
        ! criterion!
        !---------------------------------------------------------------------------
        IF (RCI_REQUEST.EQ.2) THEN
        ! Request to the DFGMRES_GET routine to put the solution into B(N) via IPAR(13)
      	    IPAR(13)=1
        ! Get the current FGMRES solution in the vector B(N)
      	     CALL DFGMRES_GET(N, COMPUTED_SOLUTION, B, RCI_REQUEST, IPAR, DPAR, TMP, ITERCOUNT)
        ! Compute the current true residual via MKL (Sparse) BLAS routines
      	     CALL MKL_DCSRGEMV('N', N, A, IA, JA, B, RESIDUAL)
      	     CALL DAXPY(N, -1.0D0, RHS, 1, RESIDUAL, 1)
      	     DVAR=DNRM2(N, RESIDUAL, 1)
      	     IF (DVAR.LT.1.0E-3) THEN
      	        GOTO 3
      	     ELSE
      	        GOTO 1
      	     END IF
        END IF
        !---------------------------------------------------------------------------
        ! If RCI_REQUEST=3, then apply the preconditioner on the vector
        ! TMP(IPAR(22)) and put the result in vector TMP(IPAR(23))
        ! Here is the recommended usage of the result produced by ILU0 routine
        ! via standard MKL Sparse Blas solver routine mkl_dcsrtrsv.
        !---------------------------------------------------------------------------
        IF (RCI_REQUEST.EQ.3) THEN
            CALL MKL_DCSRTRSV('L','N','U',N,BILU0,IA,JA,TMP(IPAR(22)),TRVEC)
            CALL MKL_DCSRTRSV('U','N','N',N,BILU0,IA,JA,TRVEC,TMP(IPAR(23)))
            GOTO 1
        END IF
        !---------------------------------------------------------------------------
        ! If RCI_REQUEST=4, then check if the norm of the next generated vector is
        ! not zero up to rounding and computational errors. The norm is contained
        ! in DPAR(7) parameter
        !---------------------------------------------------------------------------
        IF (RCI_REQUEST.EQ.4) THEN
      	    IF (DPAR(7).LT.1.0D-12) THEN
      	        GOTO 3
      	    ELSE
      	        GOTO 1
      	    ENDIF
        !---------------------------------------------------------------------------
        ! If RCI_REQUEST=anything else, then DFGMRES subroutine failed
        ! to compute the solution vector: COMPUTED_SOLUTION(N)
        !---------------------------------------------------------------------------
        ELSE
      	    GOTO 999
        END IF
        !---------------------------------------------------------------------------
        ! Reverse Communication ends here
        ! Get the current iteration number and the FGMRES solution. (DO NOT FORGET to
        ! call DFGMRES_GET routine as computed_solution is still containing
        ! the initial guess!). Request to DFGMRES_GET to put the solution into
        ! vector COMPUTED_SOLUTION(N) via IPAR(13)
        !---------------------------------------------------------------------------
        3   IPAR(13)=0
        CALL DFGMRES_GET(N, COMPUTED_SOLUTION, RHS, RCI_REQUEST, IPAR, DPAR, TMP, ITERCOUNT)
        !---------------------------------------------------------------------------
        ! Print solution vector: COMPUTED_SOLUTION(N) and
        ! the number of iterations: ITERCOUNT
        !---------------------------------------------------------------------------
        PRINT *, 'The system has been solved'
        PRINT *, 'The following solution has been obtained:'
        DO I=1,N
            PRINT *, 'COMPUTED_SOLUTION(',I,')=',COMPUTED_SOLUTION(I)*0.00014503773773021      !format '(A18,I1,A2,E10.3)'
        END DO
!        WRITE( *,'(A)') ' '
!        WRITE( *,'(A)') 'The expected solution is:'
!        DO I=1,N
!            WRITE(*,'(A18,I1,A2,E10.3)') 'EXPECTED_SOLUTION(',I,')=',EXPECTED_SOLUTION(I)
!        END DO
        PRINT *, 'Number of iterations: ',ITERCOUNT
!        WRITE( *,'(A)') ' '

        !---------------------------------------------------------------------------
        ! Release internal MKL memory that might be used for computations
        ! NOTE: It is important to call the routine below to avoid memory leaks
        ! unless you disable MKL Memory Manager
        !---------------------------------------------------------------------------
        CALL MKL_FREEBUFFERS

!        IF(ITERCOUNT.EQ.REF_NIT.AND.DABS(REF_NORM2-NRM2).LT.1.D-6) THEN
!            WRITE( *,'(A)') ' '
!            WRITE( *,'(A)') '---------------------------------------------'
!            WRITE( *,'(A,A)') 'Fortran example of FGMRES with ILU0',' preconditioner '
!            WRITE( *,'(A,A)') 'has successfully PASSED all stages of',' computations'
!            WRITE( *,'(A)') '---------------------------------------------'
!            WRITE( *,'(A)') ' '
!            STOP
!        ELSE
!            WRITE( *,'(A,A)') 'Probably, the preconditioner was computed',' incorrectly:'
!            WRITE( *,'(A,F9.6,A,F9.6)') 'Either the preconditioner norm',NRM2,' differs from the expected norm',REF_NORM2
!            WRITE( *,'(A,I2,A,I2)'),'and/or the number of iterations ', ITERCOUNT,' differs from the expected number ', REF_NIT
!            WRITE( *,'(A)') ' '
!            WRITE( *,'(A,A)') '------------------------------------------','----------------------'
!            WRITE( *,'(A,A)') 'Unfortunately, FGMRES+ILU0 Fortran example',' has FAILED'
!            WRITE( *,'(A,A)') '------------------------------------------','----------------------'
!            WRITE( *,'(A)') ' '
!            STOP
!        END IF
        999   PRINT *, 'The solver has returned the ERROR code ',RCI_REQUEST
        !---------------------------------------------------------------------------
        ! Release internal MKL memory that might be used for computations
        ! NOTE: It is important to call the routine below to avoid memory leaks
        ! unless you disable MKL Memory Manager
        !---------------------------------------------------------------------------
!        998   WRITE( *,'(A)') ' '
!        WRITE( *,'(A,A)') '---------------------------------------------','----------------------'
!        WRITE( *,'(A,A)') 'Unfortunately, FGMRES+ILU0 Fortran example',' has FAILED'
!        WRITE( *,'(A,A)') '---------------------------------------------','----------------------'
!        WRITE( *,'(A)') ' '
        
        CALL MKL_FREEBUFFERS
        STOP 

   
    end subroutine subrinput    


end module minput