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

MODULE Solver_LAPACK

USE ABQINTERFACE
USE type_meshparameters
USE type_macro_GP_DATA
USE iso_c_binding

IMPLICIT NONE

PUBLIC

TYPE solve_soe_LAPACK
    
        REAL(KIND=AbqRK), DIMENSION(:,:),POINTER:: A
        INTEGER, DIMENSION(:),POINTER:: ipiv
    
    CONTAINS
    
        PROCEDURE :: initialize => initialize_procedure
        PROCEDURE :: factor => factor_procedure
        PROCEDURE :: solve => solve_procedure
        PROCEDURE :: finish => finish_procedure

END TYPE solve_soe_LAPACK

PRIVATE:: initialize_procedure,factor_procedure,solve_procedure,finish_procedure

CONTAINS

    SUBROUTINE initialize_procedure(solver,para,analysis,GP_DATA)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_LAPACK), INTENT(INOUT):: solver
        TYPE(meshparameters),INTENT(IN)::para
        TYPE (analysisparameters),INTENT(IN)::analysis
        TYPE(macro_GP_DATA),INTENT(INOUT)::GP_DATA

        IF (analysis%monolithic) THEN !system of equations is being saved
            IF (GP_DATA%step_indicator(1)==-1) THEN !beginning of analysis, allocate memory
                CALL GP_DATA%allocate_solver_data(para%n_ROM_modes**2+para%n_ROM_modes,0)
            END IF
            !Assign Pointers to existing storage
            CALL C_F_POINTER(C_LOC(GP_DATA%solver_real(1)),solver%A,[para%n_ROM_modes,para%n_ROM_modes])
            CALL C_F_POINTER(C_LOC(GP_DATA%solver_real(para%n_ROM_modes**2+1)),solver%ipiv,[para%n_ROM_modes])
        ELSE
            ALLOCATE(solver%A(para%n_ROM_modes,para%n_ROM_modes))
            ALLOCATE(solver%ipiv(para%n_ROM_modes))
        END IF
                
    END SUBROUTINE initialize_procedure
        
    SUBROUTINE factor_procedure(solver,para,analysis,k_matrix)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_LAPACK), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: error !returns the error code
        !this type bundles all convergence/analysis parameters
        TYPE (analysisparameters),INTENT(IN)::analysis
        TYPE(meshparameters),INTENT(IN)::para!matrix structure
        REAL(KIND=AbqRK), DIMENSION(:,:),INTENT(IN):: k_matrix
        INTEGER, DIMENSION(64*para%n_ROM_modes):: work
        INTEGER:: lwork
        
        !----------------------factorize the ROM matrix-------------------------
        solver%A=k_matrix
        solver%ipiv=0
        
        lwork=size(work) !size of work as suggested by manual
        IF (analysis%symmetric_matrix) THEN
            CALL dsytrf('U',para%n_ROM_modes,solver%A,para%n_ROM_modes,solver%ipiv,work,lwork,error)
        ELSE
            CALL dgetrf(para%n_ROM_modes,para%n_ROM_modes,solver%A,para%n_ROM_modes,solver%ipiv,error)
        END IF

        IF (error/=0) THEN !output error
            CALL STDB_ABQERR(-3_AbqIK, 'The LAPACK solver exited with error %I'//&
                             ' at the factorization stage',error, 0.0_AbqRK, ' ')
        END IF
                        
    END SUBROUTINE factor_procedure
    
    SUBROUTINE solve_procedure(solver,para,analysis,rhs,solution)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_LAPACK), INTENT(INOUT):: solver !solve_soe object
        INTEGER:: error !returns the error code
        !this type bundles all convergence/analysis parameters
        TYPE (analysisparameters),INTENT(IN)::analysis
        !matrix structure
        TYPE(meshparameters),INTENT(IN)::para
        REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(OUT):: solution !solution array
        REAL(KIND=AbqRK),DIMENSION(:,:),INTENT(IN):: rhs !righthandside
        INTEGER:: nrhs !number of righthandsides
        
        solution=rhs !solution is overwritten by the solver
        nrhs=size(rhs,2)
        IF (analysis%symmetric_matrix) THEN
            CALL dsytrs('U',para%n_ROM_modes,nrhs,solver%A,para%n_ROM_modes,solver%ipiv,solution,para%n_ROM_modes,error)
        ELSE
            CALL dgetrs('N',para%n_ROM_modes,nrhs,solver%A,para%n_ROM_modes,solver%ipiv,solution,para%n_ROM_modes,error)
        END IF
        
        IF (error/=0) THEN !output error
            CALL STDB_ABQERR(-3_AbqIK, 'The LAPACK solver exited with error %I'//&
                             'at the solve stage',[error], 0.0_AbqRK, ' ')
        END IF

    END SUBROUTINE solve_procedure
    
    
    SUBROUTINE finish_procedure(solver,analysis)
    
        IMPLICIT NONE
        
        CLASS(solve_soe_LAPACK), INTENT(INOUT):: solver !solve_soe object
        !this type bundles all convergence/analysis parameters
        TYPE (analysisparameters),INTENT(IN)::analysis

        IF (.NOT. analysis%monolithic) THEN
            DEALLOCATE(solver%A)
            DEALLOCATE(solver%ipiv)
        END IF

    END SUBROUTINE finish_procedure
    
END MODULE Solver_LAPACK
