module fvn_linear
  use fvn_common
  implicit none
  
  
  !
  ! Interfaces for matrix operators
  
  ! .x. matrix multiplication
  interface operator(.x.)
      module procedure fvn_op_s_matmul,fvn_op_d_matmul,fvn_op_c_matmul,fvn_op_z_matmul
  end interface
  
  ! .t. matrix transposition
  interface operator(.t.)
      module procedure fvn_op_s_transpose,fvn_op_d_transpose,fvn_op_c_transpose,fvn_op_z_transpose
  end interface
  
  ! .tx. transpose first operand and multiply
  interface operator(.tx.)
      module procedure fvn_op_s_tx,fvn_op_d_tx,fvn_op_c_tx,fvn_op_z_tx
  end interface
  
  ! .xt. transpose second operand and multiply
  interface operator(.xt.)
      module procedure fvn_op_s_xt,fvn_op_d_xt,fvn_op_c_xt,fvn_op_z_xt
  end interface
  
  ! .i. inverse matrix
  interface operator(.i.)
      module procedure fvn_op_s_matinv,fvn_op_d_matinv,fvn_op_c_matinv,fvn_op_z_matinv
  end interface
  
  ! .ix. inverse first operand and multiply
  interface operator(.ix.)
      module procedure fvn_op_s_ix,fvn_op_d_ix,fvn_op_c_ix,fvn_op_z_ix
  end interface
  
  ! .xi. inverse second operand and multiply
  interface operator(.xi.)
      module procedure fvn_op_s_xi,fvn_op_d_xi,fvn_op_c_xi,fvn_op_z_xi
  end interface
  
  ! .h. transpose conjugate (adjoint)
  interface operator(.h.)
      module procedure fvn_op_s_transpose,fvn_op_d_transpose,fvn_op_c_conj_transpose,fvn_op_z_conj_transpose
  end interface
  
  ! .hx. transpose conjugate first operand and multiply
  interface operator(.hx.)
      module procedure fvn_op_s_tx,fvn_op_d_tx,fvn_op_c_hx,fvn_op_z_hx
  end interface
  
  ! .xh. transpose conjugate second operand and multiply
  interface operator(.xh.)
      module procedure fvn_op_s_xt,fvn_op_d_xt,fvn_op_c_xh,fvn_op_z_xh
  end interface
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! Generic interface Definition
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  ! Matrix inversion
  interface fvn_matinv
        module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv
  end interface fvn_matinv
  
  ! Determinant
  interface fvn_det
      module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det
  end interface fvn_det
  
  ! Condition
  interface fvn_matcon
      module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon
  end interface fvn_matcon
  
  ! Eigen
  interface fvn_matev
      module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev
  end interface fvn_matev
  
  ! Least square polynomial
  interface fvn_lspoly
      module procedure fvn_s_lspoly,fvn_d_lspoly
  end interface fvn_lspoly
  
  
  contains
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Linear operators
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  !
  ! .x.
  !
  function fvn_op_s_matmul(a,b)
      implicit none

      real(kind=sp_kind), dimension(:,:),intent(in) :: a,b
      real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_matmul

      fvn_op_s_matmul=matmul(a,b)
  end function
  function fvn_op_d_matmul(a,b)
      implicit none

      real(kind=dp_kind), dimension(:,:),intent(in) :: a,b
      real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_matmul

      fvn_op_d_matmul=matmul(a,b)
  end function
  function fvn_op_c_matmul(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:),intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_matmul

      fvn_op_c_matmul=matmul(a,b)
  end function
  function fvn_op_z_matmul(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_matmul

      fvn_op_z_matmul=matmul(a,b)
  end function
  
  !
  ! .tx.
  !
  function fvn_op_s_tx(a,b)
      implicit none

      real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_s_tx

      fvn_op_s_tx=matmul(transpose(a),b)
  end function
  function fvn_op_d_tx(a,b)
      implicit none

      real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_d_tx

      fvn_op_d_tx=matmul(transpose(a),b)
  end function
  function fvn_op_c_tx(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_c_tx

      fvn_op_c_tx=matmul(transpose(a),b)
  end function
  function fvn_op_z_tx(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_z_tx

      fvn_op_z_tx=matmul(transpose(a),b)
  end function
  
  
  !
  ! .xt.
  !
  function fvn_op_s_xt(a,b)
      implicit none

      real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_s_xt

      fvn_op_s_xt=matmul(a,transpose(b))
  end function
  function fvn_op_d_xt(a,b)
      implicit none

      real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_d_xt

      fvn_op_d_xt=matmul(a,transpose(b))
  end function
  function fvn_op_c_xt(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_c_xt

      fvn_op_c_xt=matmul(a,transpose(b))
  end function
  function fvn_op_z_xt(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_z_xt

      fvn_op_z_xt=matmul(a,transpose(b))
  end function
  
  !
  ! .t.
  !
  function fvn_op_s_transpose(a)
      implicit none

      real(kind=sp_kind),dimension(:,:),intent(in) :: a
      real(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_s_transpose

      fvn_op_s_transpose=transpose(a)
  end function
  function fvn_op_d_transpose(a)
      implicit none

      real(kind=dp_kind),dimension(:,:),intent(in) :: a
      real(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_d_transpose

      fvn_op_d_transpose=transpose(a)
  end function
  function fvn_op_c_transpose(a)
      implicit none

      complex(kind=sp_kind),dimension(:,:),intent(in) :: a
      complex(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_c_transpose

      fvn_op_c_transpose=transpose(a)
  end function
  function fvn_op_z_transpose(a)
      implicit none

      complex(kind=dp_kind),dimension(:,:),intent(in) :: a
      complex(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_z_transpose

      fvn_op_z_transpose=transpose(a)
  end function
  
  !
  ! .i.
  !

  ! It seems that there's a problem with automatic arrays with gfortran
  ! in some circumstances. To allow compilation with gfortran we use here a temporary array
  ! for the call. Without that there's a warning at compile time and a segmentation fault
  ! during execution. This is odd as we double memory use.

  function fvn_op_s_matinv(a)
      implicit none

      real(kind=sp_kind),dimension(:,:),intent(in) :: a
      real(kind=sp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_s_matinv,tmp_array

      call fvn_s_matinv(size(a,1),a,tmp_array,fvn_status)
      fvn_op_s_matinv=tmp_array

  end function
  function fvn_op_d_matinv(a)
      implicit none

      real(kind=dp_kind),dimension(:,:),intent(in) :: a
      real(kind=dp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_d_matinv,tmp_array

      call fvn_d_matinv(size(a,1),a,tmp_array,fvn_status)
      fvn_op_d_matinv=tmp_array

  end function
  function fvn_op_c_matinv(a)
      implicit none

      complex(kind=sp_kind),dimension(:,:),intent(in) :: a
      complex(kind=sp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_c_matinv,tmp_array

      call fvn_c_matinv(size(a,1),a,tmp_array,fvn_status)
      fvn_op_c_matinv=tmp_array

  end function
  function fvn_op_z_matinv(a)
      implicit none

      complex(kind=dp_kind),dimension(:,:),intent(in) :: a
      complex(kind=dp_kind),dimension(size(a,1),size(a,1)) :: fvn_op_z_matinv,tmp_array

      call fvn_z_matinv(size(a,1),a,tmp_array,fvn_status)
      fvn_op_z_matinv=tmp_array

  end function
  
  !
  ! .ix.
  !
  function fvn_op_s_ix(a,b)
      implicit none

      real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_ix

      fvn_op_s_ix=matmul(fvn_op_s_matinv(a),b)
  end function
  function fvn_op_d_ix(a,b)
      implicit none

      real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_ix

      fvn_op_d_ix=matmul(fvn_op_d_matinv(a),b)
  end function
  function fvn_op_c_ix(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_ix

      fvn_op_c_ix=matmul(fvn_op_c_matinv(a),b)
  end function
  function fvn_op_z_ix(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_ix

      fvn_op_z_ix=matmul(fvn_op_z_matinv(a),b)
  end function
  
  !
  ! .xi.
  !
  function fvn_op_s_xi(a,b)
      implicit none

      real(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_s_xi

      fvn_op_s_xi=matmul(a,fvn_op_s_matinv(b))
  end function
  function fvn_op_d_xi(a,b)
      implicit none

      real(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      real(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_d_xi

      fvn_op_d_xi=matmul(a,fvn_op_d_matinv(b))
  end function
  function fvn_op_c_xi(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_c_xi

      fvn_op_c_xi=matmul(a,fvn_op_c_matinv(b))
  end function
  function fvn_op_z_xi(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,1),size(b,2)) :: fvn_op_z_xi

      fvn_op_z_xi=matmul(a,fvn_op_z_matinv(b))
  end function
  
  !
  ! .h.
  !
  function fvn_op_c_conj_transpose(a)
      implicit none

      complex(kind=sp_kind),dimension(:,:),intent(in) :: a
      complex(kind=sp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_c_conj_transpose

      fvn_op_c_conj_transpose=transpose(conjg(a))
  end function
  function fvn_op_z_conj_transpose(a)
      implicit none

      complex(kind=dp_kind),dimension(:,:),intent(in) :: a
      complex(kind=dp_kind),dimension(size(a,2),size(a,1)) :: fvn_op_z_conj_transpose

      fvn_op_z_conj_transpose=transpose(conjg(a))
  end function
  
  
  !
  ! .hx.
  !
  function fvn_op_c_hx(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_c_hx

      fvn_op_c_hx=matmul(transpose(conjg(a)),b)
  end function
  function fvn_op_z_hx(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,2),size(b,2)) :: fvn_op_z_hx

      fvn_op_z_hx=matmul(transpose(conjg(a)),b)
  end function
  
  
  !
  ! .xh.
  !
  function fvn_op_c_xh(a,b)
      implicit none

      complex(kind=sp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=sp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_c_xh

      fvn_op_c_xh=matmul(a,transpose(conjg(b)))
  end function
  function fvn_op_z_xh(a,b)
      implicit none

      complex(kind=dp_kind), dimension(:,:), intent(in) :: a,b
      complex(kind=dp_kind), dimension(size(a,1),size(b,1)) :: fvn_op_z_xh

      fvn_op_z_xh=matmul(a,transpose(conjg(b)))
  end function
  
  
  
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Identity Matrix
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function fvn_d_ident(n)
        implicit none

        integer(kind=ip_kind) :: n
        real(kind=dp_kind), dimension(n,n) :: fvn_d_ident

        real(kind=dp_kind),dimension(n*n) :: vect
        integer(kind=ip_kind) :: i

        vect=0._dp_kind
        vect(1:n*n:n+1) = 1._dp_kind

        fvn_d_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  function fvn_s_ident(n)
        implicit none

        integer(kind=ip_kind) :: n
        real(kind=sp_kind), dimension(n,n) :: fvn_s_ident

        real(kind=sp_kind),dimension(n*n) :: vect
        integer(kind=ip_kind) :: i

       vect=0._sp_kind
        vect(1:n*n:n+1) = 1._sp_kind

        fvn_s_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  function fvn_c_ident(n)
        implicit none

        integer(kind=ip_kind) :: n
        complex(kind=sp_kind), dimension(n,n) :: fvn_c_ident

        complex(kind=sp_kind),dimension(n*n) :: vect
        integer(kind=ip_kind) :: i

        vect=(0._sp_kind,0._sp_kind)
        vect(1:n*n:n+1) = (1._sp_kind,0._sp_kind)

        fvn_c_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  function fvn_z_ident(n)
        implicit none

        integer(kind=ip_kind) :: n
        complex(kind=dp_kind), dimension(n,n) :: fvn_z_ident

        complex(kind=dp_kind),dimension(n*n) :: vect
        integer(kind=ip_kind) :: i

        vect=(0._dp_kind,0._dp_kind)
        vect(1:n*n:n+1) = (1._dp_kind,0._dp_kind)

        fvn_z_ident=reshape(vect, shape = (/ n,n /))
  end function
  
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Matrix inversion subroutines
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  
  subroutine fvn_s_matinv(d,a,inva,status)
      !

      ! Matrix inversion of a real(kind=sp_kind) matrix using BLAS and LAPACK

      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=sp_kind), intent(in) :: a(d,d)
      real(kind=sp_kind), intent(out) :: inva(d,d)
      integer(kind=ip_kind), intent(out),optional :: status

      integer(kind=ip_kind), allocatable :: ipiv(:)
      real(kind=sp_kind), allocatable :: work(:)
      real(kind=sp_kind) twork(1)
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork

  
      if (present(status)) status=1
  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call scopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      ! LU factorization using LAPACK
      call sgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call sgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call sgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  subroutine fvn_d_matinv(d,a,inva,status)
      !

      ! Matrix inversion of a real(kind=dp_kind) matrix using BLAS and LAPACK

      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=dp_kind), intent(in) :: a(d,d)
      real(kind=dp_kind), intent(out) :: inva(d,d)
      integer(kind=ip_kind), intent(out),optional :: status

      integer(kind=ip_kind), allocatable :: ipiv(:)
      real(kind=dp_kind), allocatable :: work(:)
      real(kind=dp_kind) :: twork(1)
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork

  
      if (present(status)) status=1
  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call dcopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      ! LU factorization using LAPACK
      call dgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call dgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call dgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  subroutine fvn_c_matinv(d,a,inva,status)
      !

      ! Matrix inversion of a complex(kind=sp_kind) matrix using BLAS and LAPACK

      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=sp_kind), intent(in) :: a(d,d)
      complex(kind=sp_kind), intent(out) :: inva(d,d)
      integer(kind=ip_kind), intent(out),optional :: status

      integer(kind=ip_kind), allocatable :: ipiv(:)
      complex(kind=sp_kind), allocatable :: work(:)
      complex(kind=sp_kind) :: twork(1)
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork

  
       if (present(status)) status=1
  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call ccopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      
      ! LU factorization using LAPACK
      call cgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call cgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call cgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  subroutine fvn_z_matinv(d,a,inva,status)
      !

      ! Matrix inversion of a complex(kind=dp_kind)(kind=sp_kind) matrix using BLAS and LAPACK

      !
      ! d (in) : matrix rank
      ! a (in) : input matrix
      ! inva (out) : inversed matrix
      ! status (ou) : =0 if something failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=dp_kind), intent(in) :: a(d,d)
      complex(kind=dp_kind), intent(out) :: inva(d,d)
      integer(kind=ip_kind), intent(out),optional :: status

      integer(kind=ip_kind), allocatable :: ipiv(:)
      complex(kind=dp_kind), allocatable :: work(:)
      complex(kind=dp_kind) :: twork(1)
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork

  
       if (present(status)) status=1
  
      allocate(ipiv(d))
      ! copy a into inva using BLAS
      !call zcopy(d*d,a,1,inva,1)
      inva(:,:)=a(:,:)
      
      ! LU factorization using LAPACK
      call zgetrf(d,d,inva,d,ipiv,info)
      ! if info is not equal to 0, something went wrong we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
          deallocate(ipiv)
          return
      end if
      ! we use the query fonction of xxxtri to obtain the optimal workspace size
      call zgetri(d,inva,d,ipiv,twork,-1,info)
      lwork=int(twork(1))
      allocate(work(lwork))
      ! Matrix inversion using LAPACK
      call zgetri(d,inva,d,ipiv,work,lwork,info)
      ! again if info is not equal to 0, we exit setting status to 0
      if (info /= 0) then
          if (present(status)) status=0
      end if
      deallocate(work)
      deallocate(ipiv)
  end subroutine
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Determinants
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  function fvn_s_det(d,a,status)
      !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! status (out) : =0 if LU factorization failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=sp_kind), intent(in) :: a(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      real(kind=sp_kind) :: fvn_s_det

      real(kind=sp_kind), allocatable :: wc_a(:,:)
      integer(kind=ip_kind), allocatable :: ipiv(:)
      integer(kind=ip_kind) :: info,i

      
       if (present(status)) status=1
      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call sgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then
           if (present(status)) status=0
          fvn_s_det=0.e0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_s_det=1.e0
      do i=1,d
          if (ipiv(i)==i) then
              fvn_s_det=fvn_s_det*wc_a(i,i)
          else
              fvn_s_det=-fvn_s_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  function fvn_d_det(d,a,status)
      !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! status (out) : =0 if LU factorization failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=dp_kind), intent(in) :: a(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      real(kind=dp_kind) :: fvn_d_det

      real(kind=dp_kind), allocatable :: wc_a(:,:)
      integer(kind=ip_kind), allocatable :: ipiv(:)
      integer(kind=ip_kind) :: info,i

      
       if (present(status)) status=1
      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call dgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then
           if (present(status)) status=0
          fvn_d_det=0.d0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_d_det=1.d0
      do i=1,d
          if (ipiv(i)==i) then
              fvn_d_det=fvn_d_det*wc_a(i,i)
          else
              fvn_d_det=-fvn_d_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  function fvn_c_det(d,a,status)    !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! status (out) : =0 if LU factorization failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=sp_kind), intent(in) :: a(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      complex(kind=sp_kind) :: fvn_c_det

      complex(kind=sp_kind), allocatable :: wc_a(:,:)
      integer(kind=ip_kind), allocatable :: ipiv(:)
      integer(kind=ip_kind) :: info,i

      
       if (present(status)) status=1
      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call cgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then
           if (present(status)) status=0
          fvn_c_det=(0.e0,0.e0)
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_c_det=(1.e0,0.e0)
      do i=1,d
          if (ipiv(i)==i) then
              fvn_c_det=fvn_c_det*wc_a(i,i)
          else
              fvn_c_det=-fvn_c_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  function fvn_z_det(d,a,status)
      !
      ! Evaluate the determinant of a square matrix using lapack LU factorization
      !
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! det (out) : determinant
      ! status (out) : =0 if LU factorization failed
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=dp_kind), intent(in) :: a(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      complex(kind=dp_kind) :: fvn_z_det

      complex(kind=dp_kind), allocatable :: wc_a(:,:)
      integer(kind=ip_kind), allocatable :: ipiv(:)
      integer(kind=ip_kind) :: info,i

      
       if (present(status)) status=1
      allocate(wc_a(d,d))
      allocate(ipiv(d))
      wc_a(:,:)=a(:,:)
      call zgetrf(d,d,wc_a,d,ipiv,info)
      if (info/= 0) then
           if (present(status)) status=0
          fvn_z_det=(0.d0,0.d0)
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      fvn_z_det=(1.d0,0.d0)
      do i=1,d
          if (ipiv(i)==i) then
              fvn_z_det=fvn_z_det*wc_a(i,i)
          else
              fvn_z_det=-fvn_z_det*wc_a(i,i)
          end if
      end do
      deallocate(ipiv)
      deallocate(wc_a)
  
  end function
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Condition test
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 1-norm 
  ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
  ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
  !
  subroutine fvn_s_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=sp_kind), intent(in) :: a(d,d)
      real(kind=sp_kind), intent(out) :: rcond
      integer(kind=ip_kind), intent(out), optional :: status

      real(kind=sp_kind), allocatable :: work(:)
      integer(kind=ip_kind), allocatable :: iwork(:)
      real(kind=sp_kind) :: anorm
      real(kind=sp_kind), allocatable :: wc_a(:,:) ! working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind), allocatable :: ipiv(:)

      real(kind=sp_kind), external :: slange

      
      
       if (present(status)) status=1
      
      anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call scopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call sgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then
           if (present(status)) status=0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      allocate(work(4*d))
      allocate(iwork(d))
      call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
      if (info /= 0) then
           if (present(status)) status=0
      end if
      deallocate(iwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_d_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=dp_kind), intent(in) :: a(d,d)
      real(kind=dp_kind), intent(out) :: rcond
      integer(kind=ip_kind), intent(out), optional :: status

      real(kind=dp_kind), allocatable :: work(:)
      integer(kind=ip_kind), allocatable :: iwork(:)
      real(kind=dp_kind) :: anorm
      real(kind=dp_kind), allocatable :: wc_a(:,:) ! working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind), allocatable :: ipiv(:)

      real(kind=dp_kind), external :: dlange

      
      
       if (present(status)) status=1
      
      anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call dcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call dgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then
           if (present(status)) status=0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
  
      allocate(work(4*d))
      allocate(iwork(d))
      call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
      if (info /= 0) then
           if (present(status)) status=0
      end if
      deallocate(iwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  
  end subroutine
  
  subroutine fvn_c_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=sp_kind), intent(in) :: a(d,d)
      real(kind=sp_kind), intent(out) :: rcond
      integer(kind=ip_kind), intent(out), optional :: status

      real(kind=sp_kind), allocatable :: rwork(:)
      complex(kind=sp_kind), allocatable :: work(:)
      integer(kind=ip_kind), allocatable :: iwork(:)
      real(kind=sp_kind) :: anorm
      complex(kind=sp_kind), allocatable :: wc_a(:,:) ! working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind), allocatable :: ipiv(:)

      real(kind=sp_kind), external :: clange

      
      
       if (present(status)) status=1
      
      anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call ccopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call cgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then
           if (present(status)) status=0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
      allocate(work(2*d))
      allocate(rwork(2*d))
      call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
      if (info /= 0) then
           if (present(status)) status=0
      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  end subroutine
  
  subroutine fvn_z_matcon(d,a,rcond,status)
      ! Matrix condition (reciprocal of condition number)
      ! 
      ! d (in) : matrix rank
      ! a (in) : The Matrix
      ! rcond (out) : guess what
      ! status (out) : =0 if something went wrong
      !
      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=dp_kind), intent(in) :: a(d,d)
      real(kind=dp_kind), intent(out) :: rcond
      integer(kind=ip_kind), intent(out), optional :: status

      complex(kind=dp_kind), allocatable :: work(:)
      real(kind=dp_kind), allocatable :: rwork(:)
      real(kind=dp_kind) :: anorm
      complex(kind=dp_kind), allocatable :: wc_a(:,:) ! working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind), allocatable :: ipiv(:)

      real(kind=dp_kind), external :: zlange

      
      
       if (present(status)) status=1
      
      anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
      
      allocate(wc_a(d,d))
      !call zcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(ipiv(d))
      call zgetrf(d,d,wc_a,d,ipiv,info)
      if (info /= 0) then
           if (present(status)) status=0
          deallocate(ipiv)
          deallocate(wc_a)
          return
      end if
  
      allocate(work(2*d))
      allocate(rwork(2*d))
      call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
      if (info /= 0) then
           if (present(status)) status=0
      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(ipiv)
      deallocate(wc_a)
  end subroutine
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  ! Valeurs propres/ Vecteurs propre
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  ! August 2009
  ! William Daniau
  ! Adding sorting of eigenvalues and vectors
  subroutine fvn_z_sort_eigen(d,v,vec)
  ! this routine takes in input :
  ! v : a vector containing unsorted eigenvalues
  ! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
  ! 
  ! At the end of subroutine the eigenvalues are sorted by decreasing order of 
  ! modulus so as the vectors.
  implicit none

  integer(kind=ip_kind) :: d,i,j
  complex(kind=dp_kind),dimension(d) :: v
  complex(kind=dp_kind),dimension(d,d) :: vec
  complex(kind=dp_kind) :: cur_v
  complex(kind=dp_kind), dimension(d) :: cur_vec 

  
  do i=2,d
    cur_v=v(i)
    cur_vec=vec(:,i)
    j=i-1
    do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
      v(j+1)=v(j)
      vec(:,j+1)=vec(:,j)
      j=j-1
    end do
    v(j+1)=cur_v
    vec(:,j+1)=cur_vec
  end do
  
  end subroutine
  
  subroutine fvn_c_sort_eigen(d,v,vec)
  ! this routine takes in input :
  ! v : a vector containing unsorted eigenvalues
  ! vec : a matrix where vec(:,j) is the eigenvector corresponding to v(j)
  ! 
  ! At the end of subroutine the eigenvalues are sorted by decreasing order of 
  ! modulus so as the vectors.
  implicit none

  integer(kind=ip_kind) :: d,i,j
  complex(kind=sp_kind),dimension(d) :: v
  complex(kind=sp_kind),dimension(d,d) :: vec
  complex(kind=sp_kind) :: cur_v
  complex(kind=sp_kind), dimension(d) :: cur_vec 

  
  do i=2,d
    cur_v=v(i)
    cur_vec=vec(:,i)
    j=i-1
    do while( (j>=1) .and. (abs(v(j)) < abs(cur_v)) )
      v(j+1)=v(j)
      vec(:,j+1)=vec(:,j)
      j=j-1
    end do
    v(j+1)=cur_v
    vec(:,j+1)=cur_vec
  end do
  
  end subroutine

  subroutine fvn_s_matev(d,a,evala,eveca,status,sortval)

      ! 

      ! integer(kind=ip_kind) d (in) : matrice rank
      ! real(kind=sp_kind) a(d,d) (in) : The Matrix
      ! complex(kind=sp_kind) evala(d) (out) : eigenvalues
      ! complex(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer(kind=ip_kind) (out) : status =0 if something went wrong

      !
      ! interfacing Lapack routine SGEEV

      !   2020 11 03 : SGEEVX with BALANC='N'

      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=sp_kind), intent(in) :: a(d,d)
      complex(kind=sp_kind), intent(out) :: evala(d)
      complex(kind=sp_kind), intent(out) :: eveca(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      logical(kind=ip_kind), intent(in), optional :: sortval

      real(kind=sp_kind), allocatable :: wc_a(:,:)  ! a working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork
      real(kind=sp_kind), allocatable :: wr(:),wi(:)
      real(kind=sp_kind) :: vl      ! unused but necessary for the call
      real(kind=sp_kind), allocatable :: vr(:,:)
      real(kind=sp_kind), allocatable :: work(:)
      real(kind=sp_kind) :: twork(1)
      integer(kind=ip_kind) i
      integer(kind=ip_kind) j

      integer(kind=ip_kind) :: ilo,ihi, iwork
      real(kind=sp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
      real(kind=sp_kind) :: abnrm
      

       if (present(status)) status=1
  
      ! making a working copy of a
      allocate(wc_a(d,d))
      !call scopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(wr(d))
      allocate(wi(d))
      allocate(vr(d,d))

      allocate(scal(d),rconde(d),rcondv(d))    

      ! query optimal work size

      ! call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
      call sgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,IWORK,info)

      lwork=int(twork(1))
      allocate(work(lwork))

      ! call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
      call sgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,IWORK,info)

      if (info /= 0) then
           if (present(status)) status=0
          deallocate(work)
          deallocate(vr)
          deallocate(wi)
          deallocate(wr)
          deallocate(wc_a)

          deallocate(scal,rconde,rcondv)

          return
      end if
  
      ! now fill in the results
      i=1
      do while(i<=d)
          evala(i)=cmplx(wr(i),wi(i))
          if (wi(i) == 0.) then ! eigenvalue is real
              eveca(:,i)=cmplx(vr(:,i),0.)
          else ! eigenvalue is complex
              evala(i+1)=cmplx(wr(i+1),wi(i+1))
              eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
              eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
              i=i+1
          end if
          i=i+1
      enddo
      deallocate(work)
      deallocate(vr)
      deallocate(wi)
      deallocate(wr)
      deallocate(wc_a)

      deallocate(scal,rconde,rcondv)

      ! sorting

      if (present(sortval) .and. sortval) then
        call fvn_c_sort_eigen(d,evala,eveca)
      end if

  end subroutine

  subroutine fvn_d_matev(d,a,evala,eveca,status,sortval)

      ! 

      ! integer(kind=ip_kind) d (in) : matrice rank
      ! real(kind=dp_kind) a(d,d) (in) : The Matrix
      ! complex(kind=dp_kind)(kind=sp_kind) evala(d) (out) : eigenvalues
      ! complex(kind=dp_kind)(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer(kind=ip_kind) (out) : status =0 if something went wrong

      !
      ! interfacing Lapack routine DGEEV

      !   2020 11 03 : DGEEVX with BALANC='N'

      implicit none

      integer(kind=ip_kind), intent(in) :: d
      real(kind=dp_kind), intent(in) :: a(d,d)
      complex(kind=dp_kind), intent(out) :: evala(d)
      complex(kind=dp_kind), intent(out) :: eveca(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      logical(kind=ip_kind), intent(in), optional :: sortval

      real(kind=dp_kind), allocatable :: wc_a(:,:)  ! a working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork
      real(kind=dp_kind), allocatable :: wr(:),wi(:)
      real(kind=dp_kind) :: vl      ! unused but necessary for the call
      real(kind=dp_kind), allocatable :: vr(:,:)
      real(kind=dp_kind), allocatable :: work(:)
      real(kind=dp_kind) :: twork(1)
      integer(kind=ip_kind) i
      integer(kind=ip_kind) j

      integer(kind=ip_kind) :: ilo,ihi, iwork
      real(kind=dp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
      real(kind=dp_kind) :: abnrm
      

      if (present(status)) status=1
  
      ! making a working copy of a
      allocate(wc_a(d,d))
      !call dcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      
      allocate(wr(d))
      allocate(wi(d))
      allocate(vr(d,d))

      allocate(scal(d),rconde(d),rcondv(d))

      ! query optimal work size

      ! call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
      call dgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,IWORK,info)

      lwork=int(twork(1))
      allocate(work(lwork))

      ! call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
      call dgeevx('N','N','V','N',d,wc_a,d,wr,wi,vl,1,vr,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,IWORK,info)

  
      if (info /= 0) then
           if (present(status)) status=0
          deallocate(work)
          deallocate(vr)
          deallocate(wi)
          deallocate(wr)
          deallocate(wc_a)

          deallocate(scal,rconde,rcondv)

          return
      end if
  
      ! now fill in the results
      i=1
      do while(i<=d)

          evala(i)=cmplx(wr(i),wi(i),dp_kind)

          if (wi(i) == 0.) then ! eigenvalue is real

              eveca(:,i)=cmplx(vr(:,i),0._dp_kind,dp_kind)

          else ! eigenvalue is complex

              evala(i+1)=cmplx(wr(i+1),wi(i+1),dp_kind)
              eveca(:,i)=cmplx(vr(:,i),vr(:,i+1),dp_kind)
              eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1),dp_kind)

              i=i+1
          end if
          i=i+1
      enddo
  
      deallocate(work)
      deallocate(vr)
      deallocate(wi)
      deallocate(wr)
      deallocate(wc_a)

      deallocate(scal,rconde,rcondv)

      ! sorting

      if (present(sortval) .and. sortval) then 
        call fvn_z_sort_eigen(d,evala,eveca)
      end if

  end subroutine

  subroutine fvn_c_matev(d,a,evala,eveca,status,sortval)

      ! 

      ! integer(kind=ip_kind) d (in) : matrice rank
      ! complex(kind=sp_kind) a(d,d) (in) : The Matrix
      ! complex(kind=sp_kind) evala(d) (out) : eigenvalues
      ! complex(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer(kind=ip_kind) (out) : status =0 if something went wrong

      !
      ! interfacing Lapack routine CGEEV

      !   2020 11 03 : CGEEVX with BALANC='N'

      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=sp_kind), intent(in) :: a(d,d)
      complex(kind=sp_kind), intent(out) :: evala(d)
      complex(kind=sp_kind), intent(out) :: eveca(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      logical(kind=ip_kind), intent(in), optional :: sortval
  
      complex(kind=sp_kind), allocatable :: wc_a(:,:) ! a working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork
      complex(kind=sp_kind), allocatable :: work(:)
      complex(kind=sp_kind) :: twork(1)
      real(kind=sp_kind), allocatable :: rwork(:)
      complex(kind=sp_kind) :: vl   ! unused but necessary for the call

      integer(kind=ip_kind) :: ilo,ihi
      real(kind=sp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
      real(kind=sp_kind) :: abnrm
      

       if (present(status)) status=1
  
      ! making a working copy of a
      allocate(wc_a(d,d))
      !call ccopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      

      ! rwork must be allocated before query
      allocate(rwork(2*d))

      allocate(scal(d),rconde(d),rcondv(d))

      ! query optimal work size

      ! call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
      call cgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,rwork,info)

      lwork=int(twork(1))

      allocate(work(lwork))    

      ! call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
      call cgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,rwork,info)

      
      if (info /= 0) then
           if (present(status)) status=0
      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(wc_a)

      deallocate(scal,rconde,rcondv)

      ! sorting

      if (present(sortval) .and. sortval) then
        call fvn_c_sort_eigen(d,evala,eveca)
      end if

  end subroutine

  subroutine fvn_z_matev(d,a,evala,eveca,status,sortval)

      ! 

      ! integer(kind=ip_kind) d (in) : matrice rank
      ! complex(kind=dp_kind)(kind=sp_kind) a(d,d) (in) : The Matrix
      ! complex(kind=dp_kind)(kind=sp_kind) evala(d) (out) : eigenvalues
      ! complex(kind=dp_kind)(kind=sp_kind) eveca(d,d) (out) : eveca(:,j) = jth eigenvector
      ! integer(kind=ip_kind) (out) : status =0 if something went wrong

      !
      ! interfacing Lapack routine ZGEEV

      !   2020 11 03 : ZGEEVX with BALANC='N'

      implicit none

      integer(kind=ip_kind), intent(in) :: d
      complex(kind=dp_kind), intent(in) :: a(d,d)
      complex(kind=dp_kind), intent(out) :: evala(d)
      complex(kind=dp_kind), intent(out) :: eveca(d,d)
      integer(kind=ip_kind), intent(out), optional :: status
      logical(kind=ip_kind), intent(in), optional :: sortval
  
      complex(kind=dp_kind), allocatable :: wc_a(:,:) ! a working copy of a
      integer(kind=ip_kind) :: info
      integer(kind=ip_kind) :: lwork
      complex(kind=dp_kind), allocatable :: work(:)
      complex(kind=dp_kind) :: twork(1)
      real(kind=dp_kind), allocatable :: rwork(:)
      complex(kind=dp_kind) :: vl   ! unused but necessary for the call

      integer(kind=ip_kind) :: ilo,ihi
      real(kind=dp_kind), allocatable, dimension(:) :: scal, rconde, rcondv
      real(kind=dp_kind) :: abnrm
  
      

       if (present(status)) status=1
  
      ! making a working copy of a
      allocate(wc_a(d,d))
      !call zcopy(d*d,a,1,wc_a,1)
      wc_a(:,:)=a(:,:)
      

      ! rwork must be allocated before query
      allocate(rwork(2*d))

      allocate(scal(d),rconde(d),rcondv(d))

      ! query optimal work size

      ! call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
      call zgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,twork,-1,rwork,info)

      lwork=int(twork(1))

      allocate(work(lwork))    

      ! call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
      call zgeevx('N','N','V','N',d,wc_a,d,evala,vl,1,eveca,d,ILO,IHI,SCAL,ABNRM,RCONDE,RCONDV,work,lwork,rwork,info)

      
      if (info /= 0) then
           if (present(status)) status=0
      end if
      deallocate(rwork)
      deallocate(work)
      deallocate(wc_a)

      deallocate(scal,rconde,rcondv)

      ! sorting

      if (present(sortval) .and. sortval) then
        call fvn_z_sort_eigen(d,evala,eveca)
      end if

  end subroutine
  
  
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  ! 
  ! Least square problem
  !
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  !
  !
  
  
  
  
  subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
  !
  !   Least square polynomial fitting
  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none

  integer(kind=ip_kind), intent(in) :: np,deg
  real(kind=dp_kind), intent(in), dimension(np) :: x,y
  real(kind=dp_kind), intent(out), dimension(deg+1) :: coeff
  integer(kind=ip_kind), intent(out), optional :: status

  real(kind=dp_kind), allocatable, dimension(:,:) :: mat,bmat
  real(kind=dp_kind),dimension(:),allocatable :: work
  real(kind=dp_kind),dimension(1) :: twork
  integer(kind=ip_kind) :: lwork,info

  integer(kind=ip_kind) :: i,j

  
   if (present(status)) status=1
  allocate(mat(np,deg+1),bmat(np,1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

  lwork=int(twork(1))
  allocate(work(lwork))

  ! real(kind=sp_kind) call

  call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
  
  if (info /= 0) then
       if (present(status)) status=0
  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat)
  end subroutine
  
  subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)
  !
  !   Least square polynomial fitting
  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none

  integer(kind=ip_kind), intent(in) :: np,deg
  real(kind=sp_kind), intent(in), dimension(np) :: x,y
  real(kind=sp_kind), intent(out), dimension(deg+1) :: coeff
  integer(kind=ip_kind), intent(out), optional :: status

  real(kind=sp_kind), allocatable, dimension(:,:) :: mat,bmat
  real(kind=sp_kind),dimension(:),allocatable :: work
  real(kind=sp_kind),dimension(1) :: twork
  integer(kind=ip_kind) :: lwork,info

  integer(kind=ip_kind) :: i,j

  
   if (present(status)) status=1
  allocate(mat(np,deg+1),bmat(np,1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)

  lwork=int(twork(1))
  allocate(work(lwork))

  ! real(kind=sp_kind) call

  call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
  
  if (info /= 0) then
       if (present(status)) status=0
  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat)
  end subroutine
  
  
  
  
  
  
  
  
  subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)
  !
  !   Least square polynomial fitting using singular value decomposition
  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none

  integer(kind=ip_kind), intent(in) :: np,deg
  real(kind=dp_kind), intent(in), dimension(np) :: x,y
  real(kind=dp_kind), intent(out), dimension(deg+1) :: coeff
  integer(kind=ip_kind), intent(out), optional :: status

  real(kind=dp_kind), allocatable, dimension(:,:) :: mat,bmat
  real(kind=dp_kind),dimension(:),allocatable :: work,singval
  real(kind=dp_kind),dimension(1) :: twork
  integer(kind=ip_kind) :: lwork,info,rank

  integer(kind=ip_kind) :: i,j

  
   if (present(status)) status=1
  allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

  lwork=int(twork(1))
  allocate(work(lwork))

  ! real(kind=sp_kind) call

  call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
  
  if (info /= 0) then
       if (present(status)) status=0
  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat,singval)
  end subroutine
  
  subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)
  !
  !   Least square polynomial fitting using singular value decomposition
  !
  !   Find the coefficients of the least square polynomial of a given degree
  !   for a set of coordinates.
  !
  !   The degree must be lower than the number of points
  !
  !   np (in)             : number of points
  !   x(np) (in)          : x data
  !   y(np) (in)          : y data
  !   deg (in)            : polynomial's degree
  !   coeff(deg+1) (out)  : polynomial coefficients
  !   status (out)        : =0 if a problem occurs
  !
  implicit none

  integer(kind=ip_kind), intent(in) :: np,deg
  real(kind=sp_kind), intent(in), dimension(np) :: x,y
  real(kind=sp_kind), intent(out), dimension(deg+1) :: coeff
  integer(kind=ip_kind), intent(out), optional :: status

  real(kind=sp_kind), allocatable, dimension(:,:) :: mat,bmat
  real(kind=sp_kind),dimension(:),allocatable :: work,singval
  real(kind=sp_kind),dimension(1) :: twork
  integer(kind=ip_kind) :: lwork,info,rank

  integer(kind=ip_kind) :: i,j

  
   if (present(status)) status=1
  allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
  
  ! Design matrix valorisation
  mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
  
  ! second member valorisation
  bmat=reshape ( (/  (y(i),i=1,np)   /)  ,shape = (/ np,1 /))
  
  ! query workspace size
  call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)

  lwork=int(twork(1))
  allocate(work(lwork))

  ! real(kind=sp_kind) call

  call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
  
  if (info /= 0) then
       if (present(status)) status=0
  end if
  
   coeff = (/ (bmat(i,1),i=1,deg+1) /)
  
  deallocate(work)
  deallocate(mat,bmat,singval)
  end subroutine

  end module fvn_linear