!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 class_Solver

USE mkl_pardiso
USE mkl_service
USE ABQINTERFACE
USE type_meshparameters

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

CONTAINS

    FUNCTION initialize_procedure(solver,para,analysis)
    
        IMPLICIT NONE
        
        CLASS(solve_soe), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: initialize_procedure !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
        !-------------------PARDISO control parameters-------------------------
        INTEGER:: mtype,phase,i
        INTEGER,DIMENSION(64):: iparm
        REAL(KIND=AbqRK),DIMENSION(1):: k_matrix_values,rhs,solution !arbitrary
        
        !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)  = 2 !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,initialize_procedure)

        IF (initialize_procedure/=0) THEN !display possible error
            CALL interpret_error(initialize_procedure)
        END IF
                
    END FUNCTION initialize_procedure
    
    
    FUNCTION factor_procedure(solver,para,analysis,k_matrix_values)
    
        IMPLICIT NONE
        
        CLASS(solve_soe), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: factor_procedure !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(para%nNonZero_reduced),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) = 2 ! 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,factor_procedure)

        IF (factor_procedure/=0) THEN !display possible error
            CALL interpret_error(factor_procedure)
        END IF
                
    END FUNCTION factor_procedure
    
    
    FUNCTION solve_procedure(solver,para,analysis,k_matrix_values,nrhs,rhs,solution)
    
        IMPLICIT NONE
        
        CLASS(solve_soe), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: solve_procedure !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(para%nNonZero_reduced),INTENT(IN):: k_matrix_values
        INTEGER,INTENT(IN):: nrhs !number of right hand sides
        REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced,nrhs),INTENT(OUT):: solution !sol. array
        REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced,nrhs),INTENT(IN):: rhs !righthandside
        REAL(KIND=AbqRK),DIMENSION(para%nRows_reduced,nrhs):: rhs_in
        !-------------------PARDISO control parameters-------------------------
        INTEGER,DIMENSION(64):: iparm
        INTEGER::phase,mtype
        
        !make copys, to avoid errors while running in parallel mode
        perm=para%perm
        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,solve_procedure)

        IF (solve_procedure/=0) THEN !display possible error
            CALL interpret_error(solve_procedure)
        END IF
        
    END FUNCTION solve_procedure
    
    
    FUNCTION finish_procedure(solver)
    
        IMPLICIT NONE
        
        CLASS(solve_soe), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: finish_procedure !returns the error code
        !-------------------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
                
        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,finish_procedure)      

        IF (finish_procedure/=0) THEN !display possible error
            CALL interpret_error(finish_procedure)
        END IF
    
    END FUNCTION finish_procedure
    
    
    FUNCTION get_permutation_matrix(para,analysis)
    
        IMPLICIT NONE
        
        TYPE(solve_soe):: solver !solve_soe object
        INTEGER:: get_permutation_matrix !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
        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 !->get permutation matrix
        
        !--------------------------iparm settings------------------------------
        iparm     =  0 !initialize iparm
        iparm(1)  =  1 !do not use default values
        iparm(2)  =  2 ! 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,get_permutation_matrix)
        para%perm=perm
        
        !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,get_permutation_matrix)
                     
        IF (get_permutation_matrix/=0) THEN !display possible error
            CALL interpret_error(get_permutation_matrix)
        END IF
        
        RETURN
        
    END FUNCTION 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 sparse solver exited with error: '//&
                                   message_str(abs(error)),0, 0.0_AbqRK, ' ')
    END IF
    
    END SUBROUTINE
    
    
END MODULE class_Solver
