!DEC$ FREEFORM
!=============================================================================
! Monolithic FE^2
! Nils Lange, Geralf Huetter, Bjoern Kiefer
!   Nils.Lange@imfd.tu-freiberg.de, Geralf.Huetter@imfd.tu-freiberg.de, 
!   Bjoern.Kiefer@imfd.tu-freiberg.de
! distributed under CC BY-NC-SA 4.0 license
! (https://creativecommons.org/licenses/by-nc-sa/4.0/)
! Reference: 
!   N. Lange, G. Huetter, B. Kiefer: "An efficient monolithic solution scheme for FE2 problems",
!   https://arxiv.org/abs/2101.01802
!
! Further information on the implementation, structure of the source code,
! examples and tutorials can be found in the file doc/documentation.pdf
! 
!=============================================================================
INCLUDE 'mkl_pardiso.f90'          !includes the PARDISO solver interface
!INCLUDE 'mkl_service.f90'          !includes some mkl service functions

MODULE Solver_PARDISO

USE mkl_pardiso
USE mkl_service
USE ABQINTERFACE
USE type_meshparameters
USE type_macro_GP_DATA

IMPLICIT NONE

PUBLIC

!declare a public object of name "solve_soe", which contains the pointer for
!the MKL PARDISO storage; the object has the methods initialize, factor
!(to factor the matrix), solve (to compute the actual solution) and finish
!(to deallocate solver storage);the routine get_permutation_matrix returns
!the permutation matrix "perm"

TYPE solve_soe_PARDISO
    
        TYPE(MKL_PARDISO_HANDLE), DIMENSION(64):: handle
    
    CONTAINS
    
        PROCEDURE :: initialize => initialize_procedure
        PROCEDURE :: factor => factor_procedure
        PROCEDURE :: solve => solve_procedure
        PROCEDURE :: finish => finish_procedure

END TYPE solve_soe_PARDISO

PRIVATE:: initialize_procedure,factor_procedure,solve_procedure,finish_procedure

CONTAINS

    SUBROUTINE initialize_procedure(solver,para,analysis,GP_DATA)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_PARDISO), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: error !returns the error code
        TYPE (analysisparameters)::analysis
        TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
        !--------------------------matrix structure----------------------------
        TYPE(meshparameters),INTENT(IN)::para
        INTEGER,DIMENSION(para%nRows_reduced):: perm
        !-------------------PARDISO control parameters-------------------------
        INTEGER:: mtype,phase,i
        INTEGER,DIMENSION(64):: iparm
        REAL(KIND=AbqRK),DIMENSION(1):: k_matrix_values,rhs,solution !arbitrary
        
        CALL mkl_free_buffers
        
        IF ((GP_DATA%step_indicator(1)==-1) .OR. (.NOT. analysis%monolithic) &
            .OR. (analysis%monolithic .AND. (.NOT. analysis%save_soe))) THEN
                
            !make copy
            perm=para%perm

            !initialize pointer with zeros
            DO i = 1, 64
                solver%handle(i)%dummy =  0 
            END DO

            !get matrix type
            IF (analysis%symmetric_matrix) THEN
                IF (analysis%indefinite_matrix) THEN
                    mtype = -2 !symmetric, indefinite matrix
                ELSE
                    mtype = 2 !symmetric, positiv definite matrix
                END IF
            ELSE
                mtype = 1 !real and structurally symmetric
            END IF

            phase = 11 !->define structure

            !--------------------------iparm settings------------------------------
            iparm     = 0
            iparm(1)  = 1 !do not use default values
            iparm(2)  = 0 !fill-in reordering from METIS
            iparm(5)  = 1 !use the supplied perm
            iparm(8)  = 2 !numbers of iterative refinement steps
            iparm(10) = 8 !Zero-Pivit handeling
            iparm(28) = 0 !double precision
            !----------------------------------------------------------------------

            !define struture
            
            CALL pardiso(solver%handle,1,1,mtype,phase,para%nRows_reduced,&
                        k_matrix_values,para%rowIndex_reduced,para%columnIndex_reduced,&
                        perm,1,iparm,0,rhs,solution,error)

            IF (error/=0) THEN !display possible error
                CALL interpret_error(error)
            END IF

            !beginning of analysis, allocate memory for the MKL handle
            IF (analysis%monolithic .AND. analysis%save_soe) THEN
                CALL GP_DATA%allocate_solver_data(0,64)
                forall(i=1:64) GP_DATA%solver_int(i)=solver%handle(i)%dummy
            END IF

        ELSE !stored system of equations, only assign pointer to data
        
            forall(i=1:64) solver%handle(i)%dummy=GP_DATA%solver_int(i)
           
        END IF
                
    END SUBROUTINE initialize_procedure
    
    
    SUBROUTINE factor_procedure(solver,para,analysis,k_matrix_values)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_PARDISO), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: error !returns the error code
        !this type bundles all convergence/analysis parameters
        TYPE (analysisparameters)::analysis
        !--------------------------matrix structure----------------------------
        TYPE(meshparameters),INTENT(IN)::para
        INTEGER,DIMENSION(para%nRows_reduced):: perm
        REAL(KIND=AbqRK), DIMENSION(:),INTENT(IN):: k_matrix_values
        !-------------------PARDISO control parameters-------------------------
        INTEGER:: mtype,phase
        INTEGER,DIMENSION(64):: iparm
        REAL(KIND=AbqRK),DIMENSION(1):: rhs,solution !arbitrary in this step
        
        !make copy
        perm=para%perm

        !get matrix type
        IF (analysis%symmetric_matrix) THEN
            IF (analysis%indefinite_matrix) THEN
                mtype = -2 !symmetric, indefinite matrix
            ELSE
                mtype = 2 !symmetric, positiv definite matrix
            END IF
        ELSE
            mtype = 1 !real and structurally symmetric
        END IF
        
        phase = 22 !->factorization
        
        !--------------------------iparm settings------------------------------
        iparm = 0
        iparm(1) = 1 !do not use default values
        iparm(2) = 0 ! fill-in reordering from METIS
        iparm(5) = 1 !use the supplied perm
        iparm(8) = 2 ! numbers of iterative refinement steps
        iparm(10) = 8 !Zero-Pivit handeling
        iparm(28) = 0 !double precision
        !----------------------------------------------------------------------

        !get factorization
        CALL pardiso(solver%handle,1,1,mtype,phase,para%nRows_reduced,&
                     k_matrix_values,para%rowIndex_reduced,para%columnIndex_reduced,&
                     perm,1,iparm,0,rhs,solution,error)
                     
        IF (error/=0) THEN !display possible error
            CALL interpret_error(error)
        END IF
                
    END SUBROUTINE factor_procedure
    
    
    SUBROUTINE solve_procedure(solver,para,analysis,k_matrix_values,rhs,solution)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_PARDISO), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: error !returns the error code
        !this type bundles all convergence/analysis parameters
        TYPE (analysisparameters)::analysis
        !--------------------------matrix structure----------------------------
        TYPE(meshparameters),INTENT(IN)::para
        INTEGER,DIMENSION(para%nRows_reduced):: perm
        REAL(KIND=AbqRK), DIMENSION(:),INTENT(IN):: k_matrix_values
        INTEGER:: nrhs !number of right hand sides
        REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: solution !sol. array
        REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN):: rhs !righthandside
        REAL(KIND=AbqRK),DIMENSION(:,:),ALLOCATABLE:: rhs_in
        INTEGER,DIMENSION(2):: rhs_shape
        !-------------------PARDISO control parameters-------------------------
        INTEGER,DIMENSION(64):: iparm
        INTEGER::phase,mtype
                
        !make copys, to avoid errors while running in parallel mode
        perm=para%perm
        rhs_shape=SHAPE(rhs)
        nrhs=rhs_shape(2)
        ALLOCATE(rhs_in(para%nRows_reduced,nrhs))
        rhs_in=rhs
        
        !get matrix type
        IF (analysis%symmetric_matrix) THEN
            IF (analysis%indefinite_matrix) THEN
                mtype = -2 !symmetric, indefinite matrix
            ELSE
                mtype = 2 !symmetric, positiv definite matrix
            END IF
        ELSE
            mtype = 1 !real and structurally symmetric
        END IF

        phase = 33 !->forward and backward substitution
        
        !--------------------------iparm settings------------------------------
        iparm = 0
        iparm(1) = 1 !do not use default values
        iparm(5) = 1 !use the supplied perm
        iparm(8) = 2 ! numbers of iterative refinement steps
        iparm(10) = 8 !Zero-Pivit handeling
        iparm(28) = 0 !double precision
        !----------------------------------------------------------------------

        solution=0.0 !initialize solution array
    
        ! Get the solution vector
        CALL pardiso(solver%handle,1,1,mtype,phase,para%nRows_reduced,&
                     k_matrix_values,para%rowIndex_reduced,para%columnIndex_reduced,&
                     perm,nrhs,iparm,0,rhs_in,solution,error)
                     
        DEALLOCATE(rhs_in)

        IF (error/=0) THEN !display possible error
            CALL interpret_error(error)
        END IF
        
    END SUBROUTINE solve_procedure
    
    
    SUBROUTINE finish_procedure(solver,para,analysis,GP_DATA)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_PARDISO), INTENT(INOUT):: solver !solve_soe object
        TYPE (analysisparameters),INTENT(IN)::analysis
        TYPE(meshparameters),INTENT(IN)::para
        TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA
        !-------------------PARDISO control parameters-------------------------
        INTEGER:: phase
        INTEGER,DIMENSION(64):: iparm
        INTEGER,DIMENSION(1)::perm !arbitrary
        INTEGER::nRows_reduced !arbitrary
        REAL(KIND=AbqRK),DIMENSION(1)::x,rhs !arbitrary
        INTEGER:: error
        
        phase = -1 !-> free memory
    
        ! Deallocate solver storage
        CALL pardiso(solver%handle,1,1,-2,phase,1,[1.0d0],[1],[1],perm,1,iparm,&
                     0,rhs,x,error)
        
        IF (error/=0) THEN !display possible error
            CALL interpret_error(error)
        END IF
    
    END SUBROUTINE finish_procedure
    
    
    SUBROUTINE get_permutation_matrix(para,analysis)
    
        IMPLICIT NONE
        
        TYPE(solve_soe_PARDISO):: solver !solve_soe object
        INTEGER:: error !returns the error code
        !this type bundles all convergence/analysis parameters
        TYPE (analysisparameters)::analysis
        !--------------------------matrix structure----------------------------
        TYPE(meshparameters),INTENT(INOUT)::para
        INTEGER,DIMENSION(para%nRows_reduced):: perm !permutation matrix
        !-------------------PARDISO control parameters-------------------------
        INTEGER:: mtype,phase,i
        INTEGER,DIMENSION(64):: iparm
        REAL(KIND=AbqRK),DIMENSION(1):: values !arbitrary
        REAL(KIND=AbqRK),DIMENSION(1):: rhs,solution !arbitrary 
        
        !initialize pointer with zeros
        FORALL(i=1:64) solver%handle(i)%dummy=0
        
        !get matrix type
        IF (analysis%symmetric_matrix) THEN
            IF (analysis%indefinite_matrix) THEN
                mtype = -2 !symmetric, indefinite matrix
            ELSE
                mtype = 2 !symmetric, positiv definite matrix
            END IF
        ELSE
            mtype = 1 !real and structurally symmetric
        END IF
        
        phase = 11 !->get permutation matrix
        
        !--------------------------iparm settings------------------------------
        iparm     =  0 !initialize iparm
        iparm(1)  =  1 !do not use default values
        iparm(2)  =  0 ! fill-in reordering from METIS
        iparm(5)  =  2 !compute perm
        iparm(19) = -1 !return number of necessary floating point operations
        iparm(28) =  0 !double precision
        !----------------------------------------------------------------------

        !get parmutation matrix
        CALL pardiso(solver%handle,1,1,mtype,phase,para%nRows_reduced,values,&
                     para%rowIndex_reduced,para%columnIndex_reduced,perm,1,&
                     iparm,0,rhs,solution,error)
        para%perm=perm

        IF (error/=0) THEN !display possible error
            CALL interpret_error(error)
        END IF
        
        !output statistical information to the message file:
        CALL STDB_ABQERR(1_AbqIK,  'Linear equation solver type: direct sparse ',0_AbqIK, 0.0_AbqRK, ' ')
        CALL STDB_ABQERR(1_AbqIK,  'Reordering of equations to minimize the wavefront.',0_AbqIK, 0.0_AbqRK, ' ')
        CALL STDB_ABQERR(1_AbqIK,  'Total number of degrees of freedom: %I ',[para%nRows_reduced], 0.0_AbqRK, ' ')
        CALL STDB_ABQERR(1_AbqIK,  'Numer of floating point operations for factorization: %I 10^6',[iparm(19)], 0.0_AbqRK, ' ')
        
        !Deallocate solver storage
        phase = -1 !-> free memory
        CALL pardiso(solver%handle,1,1,mtype,phase,para%nRows_reduced,values,&
                     para%rowIndex_reduced,para%columnIndex_reduced,perm,1,&
                     iparm,0,rhs,solution,error)
        
        RETURN
        
    END SUBROUTINE get_permutation_matrix
    
    
    SUBROUTINE interpret_error(error)
    
    IMPLICIT NONE
    
    INTEGER:: error
    CHARACTER(LEN=100),DIMENSION(12):: message_str
    
    message_str=['input inconsistent','not enough memory','reordering problem',&
                 'zero pivot, numerical factorization or iterative refinement problem',&
                 'unclassified (internal) error','reordering failed',&
                 'diagonal matrix is singular','32-bit integer overflow problem',&
                 'not enough memory for OOC','error opening OOC files',&
                 'read/write error with OOC files', 'pardiso_64 called from 32-bit library']
    
    IF (error/=0) THEN
        !write the error code in message file and end the analysis
        CALL STDB_ABQERR(-3_AbqIK, 'The PARDISO solver exited with error: '//&
                                   message_str(abs(error)),0, 0.0_AbqRK, ' ')
    END IF
    
    END SUBROUTINE
    
    
END MODULE Solver_PARDISO
