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