MODULE Expansion_Coefficients
  
USE data_structures_module, ONLY: source, MoM, post, SizeMatrix


CONTAINS


    SUBROUTINE Expansion_Coefficients_Inversion_routine(L_rnl,pivot,Gamma_qr,A)

    USE LAPACK95, ONLY: getrs
    
    IMPLICIT NONE
    
    INTEGER,INTENT(IN)     :: pivot(SizeMatrix)
    INTEGER                :: n_source, ll, n, index
    COMPLEX, INTENT(IN)    :: L_rnl(SizeMatrix,SizeMatrix)
    COMPLEX, INTENT(INOUT) :: Gamma_qr(SizeMatrix,source%N)
    COMPLEX, INTENT(OUT)   :: A(source%N,2*MoM%Max_n+1,post%N)

    CALL getrs( L_rnl, pivot, Gamma_qr )
    ! Save results
    do n_source=1,source%N
        index=1
        do ll=1,post%N
        do n=-MoM%Max_n,MoM%Max_n
            A(n_source,n+MoM%Max_n+1,ll)=Gamma_qr(index,n_source)
            index=index+1
        enddo !n
        enddo !ll
    enddo ! n_source
        
    END SUBROUTINE Expansion_Coefficients_Inversion_routine
    
    
    

    SUBROUTINE Expansion_Coefficients_Interative_routine(num_elem,L,jL,iL,Gamma_qr,A)

    !USE linear_operators
    USE MKL_PARDISO, ONLY: pardisoinit,pardiso
    USE MKL_PARDISO_PRIVATE
    
    IMPLICIT NONE
    
    INTEGER, INTENT(IN)      :: num_elem
    INTEGER                  :: ll, n, n_source, index
    INTEGER                  :: iparm(64),msglvl,error,perm(SizeMatrix)
    TYPE(MKL_PARDISO_HANDLE) :: PT(64)
    COMPLEX, INTENT(IN)      :: L(num_elem)
    INTEGER, INTENT(IN)      :: jL(num_elem),iL(SizeMatrix+1)
    COMPLEX, INTENT(INOUT)   :: Gamma_qr(SizeMatrix,source%N)
    COMPLEX                  :: solution(SizeMatrix,source%N)
    COMPLEX, INTENT(OUT)     :: A(source%N,2*MoM%Max_n+1,post%N)
    COMPLEX, ALLOCATABLE     :: vec(:,:)
    
        
        perm = 0; msglvl = 1
        CALL pardisoinit(pt, 13, iparm)
        CALL pardiso_sc(pt,1,1,13,13,SizeMatrix,L,iL,jL,perm, source%N ,iparm,msglvl, Gamma_qr ,solution,error)
        !Save results
        do n_source=1,source%N
          index=1
          do ll=1,post%N
            do n=-MoM%Max_n,MoM%Max_n
              A(n_source,n+MoM%Max_n+1,ll)=solution(index,n_source)
              index=index+1
            enddo !n
          enddo !ll
        enddo ! n_source
 
    END SUBROUTINE Expansion_Coefficients_Interative_routine



END MODULE Expansion_Coefficients