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(4), dimension(:,:),intent(in) :: a,b
    real(4), 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(8), dimension(:,:),intent(in) :: a,b
    real(8), 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(4), dimension(:,:),intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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(4), dimension(:,:), intent(in) :: a,b
    real(4), 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(8), dimension(:,:), intent(in) :: a,b
    real(8), 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(4), dimension(:,:), intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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(4), dimension(:,:), intent(in) :: a,b
    real(4), 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(8), dimension(:,:), intent(in) :: a,b
    real(8), 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(4), dimension(:,:), intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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(4),dimension(:,:),intent(in) :: a
    real(4),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(8),dimension(:,:),intent(in) :: a
    real(8),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(4),dimension(:,:),intent(in) :: a
    complex(4),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(8),dimension(:,:),intent(in) :: a
    complex(8),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(4),dimension(:,:),intent(in) :: a
    real(4),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(8),dimension(:,:),intent(in) :: a
    real(8),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(4),dimension(:,:),intent(in) :: a
    complex(4),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(8),dimension(:,:),intent(in) :: a
    complex(8),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(4), dimension(:,:), intent(in) :: a,b
    real(4), 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(8), dimension(:,:), intent(in) :: a,b
    real(8), 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(4), dimension(:,:), intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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(4), dimension(:,:), intent(in) :: a,b
    real(4), 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(8), dimension(:,:), intent(in) :: a,b
    real(8), 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(4), dimension(:,:), intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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(4),dimension(:,:),intent(in) :: a
    complex(4),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(8),dimension(:,:),intent(in) :: a
    complex(8),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(4), dimension(:,:), intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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(4), dimension(:,:), intent(in) :: a,b
    complex(4), 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(8), dimension(:,:), intent(in) :: a,b
    complex(8), 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=4) :: n
      real(kind=8), dimension(n,n) :: fvn_d_ident

      real(kind=8),dimension(n*n) :: vect
      integer(kind=4) :: i

      vect=0._8
      vect(1:n*n:n+1) = 1._8
      fvn_d_ident=reshape(vect, shape = (/ n,n /))
end function

function fvn_s_ident(n)
      implicit none
      integer(kind=4) :: n
      real(kind=4), dimension(n,n) :: fvn_s_ident

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

     vect=0._4
      vect(1:n*n:n+1) = 1._4
      fvn_s_ident=reshape(vect, shape = (/ n,n /))
end function

function fvn_c_ident(n)
      implicit none
      integer(kind=4) :: n
      complex(kind=4), dimension(n,n) :: fvn_c_ident

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

      vect=(0._4,0._4)
      vect(1:n*n:n+1) = (1._4,0._4)
      fvn_c_ident=reshape(vect, shape = (/ n,n /))
end function

function fvn_z_ident(n)
      implicit none
      integer(kind=4) :: n
      complex(kind=8), dimension(n,n) :: fvn_z_ident

      complex(kind=8),dimension(n*n) :: vect
      integer(kind=4) :: i

      vect=(0._8,0._8)
      vect(1:n*n:n+1) = (1._8,0._8)
      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 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, intent(in) :: d
    real, intent(in) :: a(d,d)
    real, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    real, allocatable :: work(:)
    real twork(1)
    integer :: info
    integer :: 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 double precision 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, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double precision, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    double precision, allocatable :: work(:)
    double precision :: twork(1)
    integer :: info
    integer :: 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 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, intent(in) :: d
    complex, intent(in) :: a(d,d)
    complex, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    complex, allocatable :: work(:)
    complex :: twork(1)
    integer :: info
    integer :: 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 double complex 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, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double complex, intent(out) :: inva(d,d)
    integer, intent(out),optional :: status

    integer, allocatable :: ipiv(:)
    double complex, allocatable :: work(:)
    double complex :: twork(1)
    integer :: info
    integer :: 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, intent(in) :: d
    real, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    real :: fvn_s_det
    
    real, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: 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, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    double precision :: fvn_d_det
    
    double precision, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: 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, intent(in) :: d
    complex, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    complex :: fvn_c_det
    
    complex, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: 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, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    integer, intent(out), optional :: status
    double complex :: fvn_z_det
    
    double complex, allocatable :: wc_a(:,:)
    integer, allocatable :: ipiv(:)
    integer :: 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, intent(in) :: d
    real, intent(in) :: a(d,d)
    real, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    real, allocatable :: work(:)
    integer, allocatable :: iwork(:)
    real :: anorm
    real, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    real, 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, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double precision, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    double precision, allocatable :: work(:)
    integer, allocatable :: iwork(:)
    double precision :: anorm
    double precision, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    double precision, 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, intent(in) :: d
    complex, intent(in) :: a(d,d)
    real, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    real, allocatable :: rwork(:)
    complex, allocatable :: work(:)
    integer, allocatable :: iwork(:)
    real :: anorm
    complex, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    real, 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, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double precision, intent(out) :: rcond
    integer, intent(out), optional :: status
    
    double complex, allocatable :: work(:)
    double precision, allocatable :: rwork(:)
    double precision :: anorm
    double complex, allocatable :: wc_a(:,:) ! working copy of a
    integer :: info
    integer, allocatable :: ipiv(:)
    
    double precision, 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 :: d,i,j
complex(8),dimension(d) :: v
complex(8),dimension(d,d) :: vec
complex(8) :: cur_v
complex(8), 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 :: d,i,j
complex(4),dimension(d) :: v
complex(4),dimension(d,d) :: vec
complex(4) :: cur_v
complex(4), 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 d (in) : matrice rank
    ! real a(d,d) (in) : The Matrix
    ! complex evala(d) (out) : eigenvalues
    ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine SGEEV
    implicit none
    integer, intent(in) :: d
    real, intent(in) :: a(d,d)
    complex, intent(out) :: evala(d)
    complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status
    logical, intent(in), optional :: sortval
    
    real, allocatable :: wc_a(:,:)  ! a working copy of a
    integer :: info
    integer :: lwork
    real, allocatable :: wr(:),wi(:)
    real :: vl      ! unused but necessary for the call
    real, allocatable :: vr(:,:)
    real, allocatable :: work(:)
    real :: twork(1)
    integer i
    integer j
    
     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))
    ! query optimal work size
    call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,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)

    if (info /= 0) then
         if (present(status)) status=0
        deallocate(work)
        deallocate(vr)
        deallocate(wi)
        deallocate(wr)
        deallocate(wc_a)
        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)

    ! 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 d (in) : matrice rank
    ! double precision a(d,d) (in) : The Matrix
    ! double complex evala(d) (out) : eigenvalues
    ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine DGEEV
    implicit none
    integer, intent(in) :: d
    double precision, intent(in) :: a(d,d)
    double complex, intent(out) :: evala(d)
    double complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status
    logical, intent(in), optional :: sortval
    
    double precision, allocatable :: wc_a(:,:)  ! a working copy of a
    integer :: info
    integer :: lwork
    double precision, allocatable :: wr(:),wi(:)
    double precision :: vl      ! unused but necessary for the call
    double precision, allocatable :: vr(:,:)
    double precision, allocatable :: work(:)
    double precision :: twork(1)
    integer i
    integer j
    
    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))
    ! query optimal work size
    call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,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)

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

    ! now fill in the results
    i=1
    do while(i<=d)
        evala(i)=dcmplx(wr(i),wi(i))
        if (wi(i) == 0.) then ! eigenvalue is real
            eveca(:,i)=dcmplx(vr(:,i),0.)
        else ! eigenvalue is complex
            evala(i+1)=dcmplx(wr(i+1),wi(i+1))
            eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
            eveca(:,i+1)=dcmplx(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)

    ! 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 d (in) : matrice rank
    ! complex a(d,d) (in) : The Matrix
    ! complex evala(d) (out) : eigenvalues
    ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine CGEEV
    implicit none
    integer, intent(in) :: d
    complex, intent(in) :: a(d,d)
    complex, intent(out) :: evala(d)
    complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status
    logical, intent(in), optional :: sortval

    complex, allocatable :: wc_a(:,:) ! a working copy of a
    integer :: info
    integer :: lwork
    complex, allocatable :: work(:)
    complex :: twork(1)
    real, allocatable :: rwork(:)
    complex :: vl   ! unused but necessary for the call

     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(:,:)
    
    
    ! query optimal work size
    call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    allocate(rwork(2*d))
    call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
    
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(wc_a)

    ! 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 d (in) : matrice rank
    ! double complex a(d,d) (in) : The Matrix
    ! double complex evala(d) (out) : eigenvalues
    ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
    ! integer (out) : status =0 if something went wrong
    !
    ! interfacing Lapack routine ZGEEV
    implicit none
    integer, intent(in) :: d
    double complex, intent(in) :: a(d,d)
    double complex, intent(out) :: evala(d)
    double complex, intent(out) :: eveca(d,d)
    integer, intent(out), optional :: status
    logical, intent(in), optional :: sortval

    double complex, allocatable :: wc_a(:,:) ! a working copy of a
    integer :: info
    integer :: lwork
    double complex, allocatable :: work(:)
    double complex :: twork(1)
    double precision, allocatable :: rwork(:)
    double complex :: vl   ! unused but necessary for the call
    
     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(:,:)
    
    ! query optimal work size
    call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
    lwork=int(twork(1))
    allocate(work(lwork))
    allocate(rwork(2*d))
    call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
    
    if (info /= 0) then
         if (present(status)) status=0
    end if
    deallocate(rwork)
    deallocate(work)
    deallocate(wc_a)

    ! 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, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

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

integer :: 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=twork(1)
allocate(work(int(lwork)))
! real 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, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

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

integer :: 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=twork(1)
allocate(work(int(lwork)))
! real 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, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

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

integer :: 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=twork(1)
allocate(work(int(lwork)))
! real 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, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status

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

integer :: 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=twork(1)
allocate(work(int(lwork)))
! real 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