lxsd / fvn

Download as
Browse Code ยป

Commit c90ee9d32e829cd43f6e37c15eb52c112a4e401f

Authored by daniau
1 parent 81b5d24e17
Exists in master and in 3 other branches geevx, geevx_2020, multi

git-svn-id: https://lxsd.femto-st.fr/svn/fvn@32 b657c933-2333-4658-acf2-d3c7c2708721

Showing 1 changed file with 2 additions and 2 deletions Inline Diff

1 1
module fvn 2 2 module fvn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 4 4 !
! fvn : a f95 module replacement for some imsl routines 5 5 ! fvn : a f95 module replacement for some imsl routines
! it uses lapack for linear algebra 6 6 ! it uses lapack for linear algebra
! it uses modified quadpack for integration 7 7 ! it uses modified quadpack for integration
! 8 8 !
! William Daniau 2007->today 9 9 ! William Daniau 2007->today
! william.daniau@femto-st.fr 10 10 ! william.daniau@femto-st.fr
! 11 11 !
! Routines naming scheme : 12 12 ! Routines naming scheme :
! 13 13 !
! fvn_x_name 14 14 ! fvn_x_name
! where x can be s : real 15 15 ! where x can be s : real
! d : real double precision 16 16 ! d : real double precision
! c : complex 17 17 ! c : complex
! z : double complex 18 18 ! z : double complex
! 19 19 !
! 20 20 !
! This piece of code is totally free! Do whatever you want with it. However 21 21 ! This piece of code is totally free! Do whatever you want with it. However
! if you find it usefull it would be kind to give credits ;-) 22 22 ! if you find it usefull it would be kind to give credits ;-)
! 23 23 !
! svn version 24 24 ! svn version
! February 2008 : added fnlib to repository so there's no use to have 25 25 ! February 2008 : added fnlib to repository so there's no use to have
! special functions and trigonometry here. Some functions 26 26 ! special functions and trigonometry here. Some functions
! are then removed 27 27 ! are then removed
! January 2008 : added quadratic interpolation, gamma/factorial function, 28 28 ! January 2008 : added quadratic interpolation, gamma/factorial function,
! a function which return identity matrix, 29 29 ! a function which return identity matrix,
! evaluation of nterm chebyshev series 30 30 ! evaluation of nterm chebyshev series
! September 2007 : added sparse system solving by interfacing umfpack 31 31 ! September 2007 : added sparse system solving by interfacing umfpack
! June 2007 : added some complex trigonometric functions 32 32 ! June 2007 : added some complex trigonometric functions
! 33 33 !
! TO DO LIST : 34 34 ! TO DO LIST :
! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm 35 35 ! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm
! eigenvalues are given with no particular order. 36 36 ! eigenvalues are given with no particular order.
! + Generic interface for fvn_x_name family -> fvn_name 37 37 ! + Generic interface for fvn_x_name family -> fvn_name
! + Make some parameters optional, status for example 38 38 ! + Make some parameters optional, status for example
! + use f95 kinds "double complex" -> complex(kind=8) 39 39 ! + use f95 kinds "double complex" -> complex(kind=8)
! + unify quadpack routines 40 40 ! + unify quadpack routines
! + ... 41 41 ! + ...
! 42 42 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 43 43 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
44 44
implicit none 45 45 implicit none
! We define pi and i for the module 46 46 ! We define pi and i for the module
real(kind=8),parameter :: fvn_pi = 3.141592653589793_8 47 47 real(kind=8),parameter :: fvn_pi = 3.141592653589793_8
complex(kind=8),parameter :: fvn_i = (0._8,1._8) 48 48 complex(kind=8),parameter :: fvn_i = (0._8,1._8)
49 49
50 50
! All quadpack routines are private to the module 51 51 ! All quadpack routines are private to the module
private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, & 52 52 private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, &
dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, & 53 53 dqage_2d_outer,dqk15,dqk15_2d_inner,dqk15_2d_outer,dqk21,dqk21_2d_inner,dqk21_2d_outer, &
dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, & 54 54 dqk31,dqk31_2d_inner,dqk31_2d_outer,dqk41,dqk41_2d_inner,dqk41_2d_outer, &
dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt 55 55 dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt
56 56
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 57 57 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 58 58 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Generic interface Definition 59 59 ! Generic interface Definition
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 60 60 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 61 61 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
62 62
63 63
! Matrix inversion 64 64 ! Matrix inversion
interface fvn_matinv 65 65 interface fvn_matinv
module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv 66 66 module procedure fvn_s_matinv,fvn_d_matinv,fvn_c_matinv,fvn_z_matinv
end interface fvn_matinv 67 67 end interface fvn_matinv
68 68
! Determinant 69 69 ! Determinant
interface fvn_det 70 70 interface fvn_det
module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det 71 71 module procedure fvn_s_det,fvn_d_det,fvn_c_det,fvn_z_det
end interface fvn_det 72 72 end interface fvn_det
73 73
! Condition 74 74 ! Condition
interface fvn_matcon 75 75 interface fvn_matcon
module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon 76 76 module procedure fvn_s_matcon,fvn_d_matcon,fvn_c_matcon,fvn_z_matcon
end interface fvn_matcon 77 77 end interface fvn_matcon
78 78
! Eigen 79 79 ! Eigen
interface fvn_matev 80 80 interface fvn_matev
module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev 81 81 module procedure fvn_s_matev,fvn_d_matev,fvn_c_matev,fvn_z_matev
end interface fvn_matev 82 82 end interface fvn_matev
83 83
! Utility procedure find interval 84 84 ! Utility procedure find interval
interface fvn_find_interval 85 85 interface fvn_find_interval
module procedure fvn_s_find_interval,fvn_d_find_interval 86 86 module procedure fvn_s_find_interval,fvn_d_find_interval
end interface fvn_find_interval 87 87 end interface fvn_find_interval
88 88
! Quadratic 1D interpolation 89 89 ! Quadratic 1D interpolation
interface fvn_quad_interpol 90 90 interface fvn_quad_interpol
module procedure fvn_s_quad_interpol,fvn_d_quad_interpol 91 91 module procedure fvn_s_quad_interpol,fvn_d_quad_interpol
end interface fvn_quad_interpol 92 92 end interface fvn_quad_interpol
93 93
! Quadratic 2D interpolation 94 94 ! Quadratic 2D interpolation
interface fvn_quad_2d_interpol 95 95 interface fvn_quad_2d_interpol
module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol 96 96 module procedure fvn_s_quad_2d_interpol,fvn_d_quad_2d_interpol
end interface fvn_quad_2d_interpol 97 97 end interface fvn_quad_2d_interpol
98 98
! Quadratic 3D interpolation 99 99 ! Quadratic 3D interpolation
interface fvn_quad_3d_interpol 100 100 interface fvn_quad_3d_interpol
module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol 101 101 module procedure fvn_s_quad_3d_interpol,fvn_d_quad_3d_interpol
end interface fvn_quad_3d_interpol 102 102 end interface fvn_quad_3d_interpol
103 103
! Akima interpolation 104 104 ! Akima interpolation
interface fvn_akima 105 105 interface fvn_akima
module procedure fvn_s_akima,fvn_d_akima 106 106 module procedure fvn_s_akima,fvn_d_akima
end interface fvn_akima 107 107 end interface fvn_akima
108 108
! Akima evaluation 109 109 ! Akima evaluation
interface fvn_spline_eval 110 110 interface fvn_spline_eval
module procedure fvn_s_spline_eval,fvn_d_spline_eval 111 111 module procedure fvn_s_spline_eval,fvn_d_spline_eval
end interface fvn_spline_eval 112 112 end interface fvn_spline_eval
113 113
! Least square polynomial 114 114 ! Least square polynomial
interface fvn_lspoly 115 115 interface fvn_lspoly
module procedure fvn_s_lspoly,fvn_d_lspoly 116 116 module procedure fvn_s_lspoly,fvn_d_lspoly
end interface fvn_lspoly 117 117 end interface fvn_lspoly
118 118
! Muller 119 119 ! Muller
interface fvn_muller 120 120 interface fvn_muller
module procedure fvn_z_muller 121 121 module procedure fvn_z_muller
end interface fvn_muller 122 122 end interface fvn_muller
123 123
! Gauss legendre 124 124 ! Gauss legendre
interface fvn_gauss_legendre 125 125 interface fvn_gauss_legendre
module procedure fvn_d_gauss_legendre 126 126 module procedure fvn_d_gauss_legendre
end interface fvn_gauss_legendre 127 127 end interface fvn_gauss_legendre
128 128
! Simple Gauss Legendre integration 129 129 ! Simple Gauss Legendre integration
interface fvn_gl_integ 130 130 interface fvn_gl_integ
module procedure fvn_d_gl_integ 131 131 module procedure fvn_d_gl_integ
end interface fvn_gl_integ 132 132 end interface fvn_gl_integ
133 133
! Adaptative Gauss Kronrod integration f(x) 134 134 ! Adaptative Gauss Kronrod integration f(x)
interface fvn_integ_1_gk 135 135 interface fvn_integ_1_gk
module procedure fvn_d_integ_1_gk 136 136 module procedure fvn_d_integ_1_gk
end interface fvn_integ_1_gk 137 137 end interface fvn_integ_1_gk
138 138
! Adaptative Gauss Kronrod integration f(x,y) 139 139 ! Adaptative Gauss Kronrod integration f(x,y)
interface fvn_integ_2_gk 140 140 interface fvn_integ_2_gk
module procedure fvn_d_integ_2_gk 141 141 module procedure fvn_d_integ_2_gk
end interface fvn_integ_2_gk 142 142 end interface fvn_integ_2_gk
143 143
! Sparse solving 144 144 ! Sparse solving
interface fvn_sparse_solve 145 145 interface fvn_sparse_solve
module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve 146 146 module procedure fvn_zl_sparse_solve,fvn_zi_sparse_solve,fvn_dl_sparse_solve,fvn_di_sparse_solve
end interface fvn_sparse_solve 147 147 end interface fvn_sparse_solve
148 148
149 149
contains 150 150 contains
151 151
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 152 152 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 153 153 !
! Identity Matrix 154 154 ! Identity Matrix
! 155 155 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 156 156 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_d_ident(n) 157 157 function fvn_d_ident(n)
implicit none 158 158 implicit none
integer(kind=4) :: n 159 159 integer(kind=4) :: n
real(kind=8), dimension(n,n) :: fvn_d_ident 160 160 real(kind=8), dimension(n,n) :: fvn_d_ident
161 161
real(kind=8),dimension(n*n) :: vect 162 162 real(kind=8),dimension(n*n) :: vect
integer(kind=4) :: i 163 163 integer(kind=4) :: i
164 164
vect=0._8 165 165 vect=0._8
vect(1:n*n:n+1) = 1._8 166 166 vect(1:n*n:n+1) = 1._8
fvn_d_ident=reshape(vect, shape = (/ n,n /)) 167 167 fvn_d_ident=reshape(vect, shape = (/ n,n /))
end function 168 168 end function
169 169
function fvn_s_ident(n) 170 170 function fvn_s_ident(n)
implicit none 171 171 implicit none
integer(kind=4) :: n 172 172 integer(kind=4) :: n
real(kind=4), dimension(n,n) :: fvn_s_ident 173 173 real(kind=4), dimension(n,n) :: fvn_s_ident
174 174
real(kind=4),dimension(n*n) :: vect 175 175 real(kind=4),dimension(n*n) :: vect
integer(kind=4) :: i 176 176 integer(kind=4) :: i
177 177
vect=0._4 178 178 vect=0._4
vect(1:n*n:n+1) = 1._4 179 179 vect(1:n*n:n+1) = 1._4
fvn_s_ident=reshape(vect, shape = (/ n,n /)) 180 180 fvn_s_ident=reshape(vect, shape = (/ n,n /))
end function 181 181 end function
182 182
function fvn_c_ident(n) 183 183 function fvn_c_ident(n)
implicit none 184 184 implicit none
integer(kind=4) :: n 185 185 integer(kind=4) :: n
complex(kind=4), dimension(n,n) :: fvn_c_ident 186 186 complex(kind=4), dimension(n,n) :: fvn_c_ident
187 187
complex(kind=4),dimension(n*n) :: vect 188 188 complex(kind=4),dimension(n*n) :: vect
integer(kind=4) :: i 189 189 integer(kind=4) :: i
190 190
vect=(0._4,0._4) 191 191 vect=(0._4,0._4)
vect(1:n*n:n+1) = (1._4,0._4) 192 192 vect(1:n*n:n+1) = (1._4,0._4)
fvn_c_ident=reshape(vect, shape = (/ n,n /)) 193 193 fvn_c_ident=reshape(vect, shape = (/ n,n /))
end function 194 194 end function
195 195
function fvn_z_ident(n) 196 196 function fvn_z_ident(n)
implicit none 197 197 implicit none
integer(kind=4) :: n 198 198 integer(kind=4) :: n
complex(kind=8), dimension(n,n) :: fvn_z_ident 199 199 complex(kind=8), dimension(n,n) :: fvn_z_ident
200 200
complex(kind=8),dimension(n*n) :: vect 201 201 complex(kind=8),dimension(n*n) :: vect
integer(kind=4) :: i 202 202 integer(kind=4) :: i
203 203
vect=(0._8,0._8) 204 204 vect=(0._8,0._8)
vect(1:n*n:n+1) = (1._8,0._8) 205 205 vect(1:n*n:n+1) = (1._8,0._8)
fvn_z_ident=reshape(vect, shape = (/ n,n /)) 206 206 fvn_z_ident=reshape(vect, shape = (/ n,n /))
end function 207 207 end function
208 208
209 209
210 210
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 211 211 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 212 212 !
! Matrix inversion subroutines 213 213 ! Matrix inversion subroutines
! 214 214 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 215 215 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
216 216
subroutine fvn_s_matinv(d,a,inva,status) 217 217 subroutine fvn_s_matinv(d,a,inva,status)
! 218 218 !
! Matrix inversion of a real matrix using BLAS and LAPACK 219 219 ! Matrix inversion of a real matrix using BLAS and LAPACK
! 220 220 !
! d (in) : matrix rank 221 221 ! d (in) : matrix rank
! a (in) : input matrix 222 222 ! a (in) : input matrix
! inva (out) : inversed matrix 223 223 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 224 224 ! status (ou) : =0 if something failed
! 225 225 !
implicit none 226 226 implicit none
integer, intent(in) :: d 227 227 integer, intent(in) :: d
real, intent(in) :: a(d,d) 228 228 real, intent(in) :: a(d,d)
real, intent(out) :: inva(d,d) 229 229 real, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 230 230 integer, intent(out),optional :: status
231 231
integer, allocatable :: ipiv(:) 232 232 integer, allocatable :: ipiv(:)
real, allocatable :: work(:) 233 233 real, allocatable :: work(:)
real twork(1) 234 234 real twork(1)
integer :: info 235 235 integer :: info
integer :: lwork 236 236 integer :: lwork
237 237
if (present(status)) status=1 238 238 if (present(status)) status=1
239 239
allocate(ipiv(d)) 240 240 allocate(ipiv(d))
! copy a into inva using BLAS 241 241 ! copy a into inva using BLAS
!call scopy(d*d,a,1,inva,1) 242 242 !call scopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 243 243 inva(:,:)=a(:,:)
! LU factorization using LAPACK 244 244 ! LU factorization using LAPACK
call sgetrf(d,d,inva,d,ipiv,info) 245 245 call sgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 246 246 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 247 247 if (info /= 0) then
if (present(status)) status=0 248 248 if (present(status)) status=0
deallocate(ipiv) 249 249 deallocate(ipiv)
return 250 250 return
end if 251 251 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 252 252 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call sgetri(d,inva,d,ipiv,twork,-1,info) 253 253 call sgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 254 254 lwork=int(twork(1))
allocate(work(lwork)) 255 255 allocate(work(lwork))
! Matrix inversion using LAPACK 256 256 ! Matrix inversion using LAPACK
call sgetri(d,inva,d,ipiv,work,lwork,info) 257 257 call sgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 258 258 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 259 259 if (info /= 0) then
if (present(status)) status=0 260 260 if (present(status)) status=0
end if 261 261 end if
deallocate(work) 262 262 deallocate(work)
deallocate(ipiv) 263 263 deallocate(ipiv)
end subroutine 264 264 end subroutine
265 265
subroutine fvn_d_matinv(d,a,inva,status) 266 266 subroutine fvn_d_matinv(d,a,inva,status)
! 267 267 !
! Matrix inversion of a double precision matrix using BLAS and LAPACK 268 268 ! Matrix inversion of a double precision matrix using BLAS and LAPACK
! 269 269 !
! d (in) : matrix rank 270 270 ! d (in) : matrix rank
! a (in) : input matrix 271 271 ! a (in) : input matrix
! inva (out) : inversed matrix 272 272 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 273 273 ! status (ou) : =0 if something failed
! 274 274 !
implicit none 275 275 implicit none
integer, intent(in) :: d 276 276 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 277 277 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: inva(d,d) 278 278 double precision, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 279 279 integer, intent(out),optional :: status
280 280
integer, allocatable :: ipiv(:) 281 281 integer, allocatable :: ipiv(:)
double precision, allocatable :: work(:) 282 282 double precision, allocatable :: work(:)
double precision :: twork(1) 283 283 double precision :: twork(1)
integer :: info 284 284 integer :: info
integer :: lwork 285 285 integer :: lwork
286 286
if (present(status)) status=1 287 287 if (present(status)) status=1
288 288
allocate(ipiv(d)) 289 289 allocate(ipiv(d))
! copy a into inva using BLAS 290 290 ! copy a into inva using BLAS
!call dcopy(d*d,a,1,inva,1) 291 291 !call dcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 292 292 inva(:,:)=a(:,:)
! LU factorization using LAPACK 293 293 ! LU factorization using LAPACK
call dgetrf(d,d,inva,d,ipiv,info) 294 294 call dgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 295 295 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 296 296 if (info /= 0) then
if (present(status)) status=0 297 297 if (present(status)) status=0
deallocate(ipiv) 298 298 deallocate(ipiv)
return 299 299 return
end if 300 300 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 301 301 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call dgetri(d,inva,d,ipiv,twork,-1,info) 302 302 call dgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 303 303 lwork=int(twork(1))
allocate(work(lwork)) 304 304 allocate(work(lwork))
! Matrix inversion using LAPACK 305 305 ! Matrix inversion using LAPACK
call dgetri(d,inva,d,ipiv,work,lwork,info) 306 306 call dgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 307 307 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 308 308 if (info /= 0) then
if (present(status)) status=0 309 309 if (present(status)) status=0
end if 310 310 end if
deallocate(work) 311 311 deallocate(work)
deallocate(ipiv) 312 312 deallocate(ipiv)
end subroutine 313 313 end subroutine
314 314
subroutine fvn_c_matinv(d,a,inva,status) 315 315 subroutine fvn_c_matinv(d,a,inva,status)
! 316 316 !
! Matrix inversion of a complex matrix using BLAS and LAPACK 317 317 ! Matrix inversion of a complex matrix using BLAS and LAPACK
! 318 318 !
! d (in) : matrix rank 319 319 ! d (in) : matrix rank
! a (in) : input matrix 320 320 ! a (in) : input matrix
! inva (out) : inversed matrix 321 321 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 322 322 ! status (ou) : =0 if something failed
! 323 323 !
implicit none 324 324 implicit none
integer, intent(in) :: d 325 325 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 326 326 complex, intent(in) :: a(d,d)
complex, intent(out) :: inva(d,d) 327 327 complex, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 328 328 integer, intent(out),optional :: status
329 329
integer, allocatable :: ipiv(:) 330 330 integer, allocatable :: ipiv(:)
complex, allocatable :: work(:) 331 331 complex, allocatable :: work(:)
complex :: twork(1) 332 332 complex :: twork(1)
integer :: info 333 333 integer :: info
integer :: lwork 334 334 integer :: lwork
335 335
if (present(status)) status=1 336 336 if (present(status)) status=1
337 337
allocate(ipiv(d)) 338 338 allocate(ipiv(d))
! copy a into inva using BLAS 339 339 ! copy a into inva using BLAS
!call ccopy(d*d,a,1,inva,1) 340 340 !call ccopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 341 341 inva(:,:)=a(:,:)
342 342
! LU factorization using LAPACK 343 343 ! LU factorization using LAPACK
call cgetrf(d,d,inva,d,ipiv,info) 344 344 call cgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 345 345 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 346 346 if (info /= 0) then
if (present(status)) status=0 347 347 if (present(status)) status=0
deallocate(ipiv) 348 348 deallocate(ipiv)
return 349 349 return
end if 350 350 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 351 351 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call cgetri(d,inva,d,ipiv,twork,-1,info) 352 352 call cgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 353 353 lwork=int(twork(1))
allocate(work(lwork)) 354 354 allocate(work(lwork))
! Matrix inversion using LAPACK 355 355 ! Matrix inversion using LAPACK
call cgetri(d,inva,d,ipiv,work,lwork,info) 356 356 call cgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 357 357 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 358 358 if (info /= 0) then
if (present(status)) status=0 359 359 if (present(status)) status=0
end if 360 360 end if
deallocate(work) 361 361 deallocate(work)
deallocate(ipiv) 362 362 deallocate(ipiv)
end subroutine 363 363 end subroutine
364 364
subroutine fvn_z_matinv(d,a,inva,status) 365 365 subroutine fvn_z_matinv(d,a,inva,status)
! 366 366 !
! Matrix inversion of a double complex matrix using BLAS and LAPACK 367 367 ! Matrix inversion of a double complex matrix using BLAS and LAPACK
! 368 368 !
! d (in) : matrix rank 369 369 ! d (in) : matrix rank
! a (in) : input matrix 370 370 ! a (in) : input matrix
! inva (out) : inversed matrix 371 371 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 372 372 ! status (ou) : =0 if something failed
! 373 373 !
implicit none 374 374 implicit none
integer, intent(in) :: d 375 375 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 376 376 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: inva(d,d) 377 377 double complex, intent(out) :: inva(d,d)
integer, intent(out),optional :: status 378 378 integer, intent(out),optional :: status
379 379
integer, allocatable :: ipiv(:) 380 380 integer, allocatable :: ipiv(:)
double complex, allocatable :: work(:) 381 381 double complex, allocatable :: work(:)
double complex :: twork(1) 382 382 double complex :: twork(1)
integer :: info 383 383 integer :: info
integer :: lwork 384 384 integer :: lwork
385 385
if (present(status)) status=1 386 386 if (present(status)) status=1
387 387
allocate(ipiv(d)) 388 388 allocate(ipiv(d))
! copy a into inva using BLAS 389 389 ! copy a into inva using BLAS
!call zcopy(d*d,a,1,inva,1) 390 390 !call zcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 391 391 inva(:,:)=a(:,:)
392 392
! LU factorization using LAPACK 393 393 ! LU factorization using LAPACK
call zgetrf(d,d,inva,d,ipiv,info) 394 394 call zgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 395 395 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 396 396 if (info /= 0) then
if (present(status)) status=0 397 397 if (present(status)) status=0
deallocate(ipiv) 398 398 deallocate(ipiv)
return 399 399 return
end if 400 400 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 401 401 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call zgetri(d,inva,d,ipiv,twork,-1,info) 402 402 call zgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 403 403 lwork=int(twork(1))
allocate(work(lwork)) 404 404 allocate(work(lwork))
! Matrix inversion using LAPACK 405 405 ! Matrix inversion using LAPACK
call zgetri(d,inva,d,ipiv,work,lwork,info) 406 406 call zgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 407 407 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 408 408 if (info /= 0) then
if (present(status)) status=0 409 409 if (present(status)) status=0
end if 410 410 end if
deallocate(work) 411 411 deallocate(work)
deallocate(ipiv) 412 412 deallocate(ipiv)
end subroutine 413 413 end subroutine
414 414
415 415
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 416 416 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 417 417 !
! Determinants 418 418 ! Determinants
! 419 419 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 420 420 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_s_det(d,a,status) 421 421 function fvn_s_det(d,a,status)
! 422 422 !
! Evaluate the determinant of a square matrix using lapack LU factorization 423 423 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 424 424 !
! d (in) : matrix rank 425 425 ! d (in) : matrix rank
! a (in) : The Matrix 426 426 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 427 427 ! status (out) : =0 if LU factorization failed
! 428 428 !
implicit none 429 429 implicit none
integer, intent(in) :: d 430 430 integer, intent(in) :: d
real, intent(in) :: a(d,d) 431 431 real, intent(in) :: a(d,d)
integer, intent(out), optional :: status 432 432 integer, intent(out), optional :: status
real :: fvn_s_det 433 433 real :: fvn_s_det
434 434
real, allocatable :: wc_a(:,:) 435 435 real, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 436 436 integer, allocatable :: ipiv(:)
integer :: info,i 437 437 integer :: info,i
438 438
if (present(status)) status=1 439 439 if (present(status)) status=1
allocate(wc_a(d,d)) 440 440 allocate(wc_a(d,d))
allocate(ipiv(d)) 441 441 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 442 442 wc_a(:,:)=a(:,:)
call sgetrf(d,d,wc_a,d,ipiv,info) 443 443 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 444 444 if (info/= 0) then
if (present(status)) status=0 445 445 if (present(status)) status=0
fvn_s_det=0.e0 446 446 fvn_s_det=0.e0
deallocate(ipiv) 447 447 deallocate(ipiv)
deallocate(wc_a) 448 448 deallocate(wc_a)
return 449 449 return
end if 450 450 end if
fvn_s_det=1.e0 451 451 fvn_s_det=1.e0
do i=1,d 452 452 do i=1,d
if (ipiv(i)==i) then 453 453 if (ipiv(i)==i) then
fvn_s_det=fvn_s_det*wc_a(i,i) 454 454 fvn_s_det=fvn_s_det*wc_a(i,i)
else 455 455 else
fvn_s_det=-fvn_s_det*wc_a(i,i) 456 456 fvn_s_det=-fvn_s_det*wc_a(i,i)
end if 457 457 end if
end do 458 458 end do
deallocate(ipiv) 459 459 deallocate(ipiv)
deallocate(wc_a) 460 460 deallocate(wc_a)
461 461
end function 462 462 end function
463 463
function fvn_d_det(d,a,status) 464 464 function fvn_d_det(d,a,status)
! 465 465 !
! Evaluate the determinant of a square matrix using lapack LU factorization 466 466 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 467 467 !
! d (in) : matrix rank 468 468 ! d (in) : matrix rank
! a (in) : The Matrix 469 469 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 470 470 ! status (out) : =0 if LU factorization failed
! 471 471 !
implicit none 472 472 implicit none
integer, intent(in) :: d 473 473 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 474 474 double precision, intent(in) :: a(d,d)
integer, intent(out), optional :: status 475 475 integer, intent(out), optional :: status
double precision :: fvn_d_det 476 476 double precision :: fvn_d_det
477 477
double precision, allocatable :: wc_a(:,:) 478 478 double precision, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 479 479 integer, allocatable :: ipiv(:)
integer :: info,i 480 480 integer :: info,i
481 481
if (present(status)) status=1 482 482 if (present(status)) status=1
allocate(wc_a(d,d)) 483 483 allocate(wc_a(d,d))
allocate(ipiv(d)) 484 484 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 485 485 wc_a(:,:)=a(:,:)
call dgetrf(d,d,wc_a,d,ipiv,info) 486 486 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 487 487 if (info/= 0) then
if (present(status)) status=0 488 488 if (present(status)) status=0
fvn_d_det=0.d0 489 489 fvn_d_det=0.d0
deallocate(ipiv) 490 490 deallocate(ipiv)
deallocate(wc_a) 491 491 deallocate(wc_a)
return 492 492 return
end if 493 493 end if
fvn_d_det=1.d0 494 494 fvn_d_det=1.d0
do i=1,d 495 495 do i=1,d
if (ipiv(i)==i) then 496 496 if (ipiv(i)==i) then
fvn_d_det=fvn_d_det*wc_a(i,i) 497 497 fvn_d_det=fvn_d_det*wc_a(i,i)
else 498 498 else
fvn_d_det=-fvn_d_det*wc_a(i,i) 499 499 fvn_d_det=-fvn_d_det*wc_a(i,i)
end if 500 500 end if
end do 501 501 end do
deallocate(ipiv) 502 502 deallocate(ipiv)
deallocate(wc_a) 503 503 deallocate(wc_a)
504 504
end function 505 505 end function
506 506
function fvn_c_det(d,a,status) ! 507 507 function fvn_c_det(d,a,status) !
! Evaluate the determinant of a square matrix using lapack LU factorization 508 508 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 509 509 !
! d (in) : matrix rank 510 510 ! d (in) : matrix rank
! a (in) : The Matrix 511 511 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 512 512 ! status (out) : =0 if LU factorization failed
! 513 513 !
implicit none 514 514 implicit none
integer, intent(in) :: d 515 515 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 516 516 complex, intent(in) :: a(d,d)
integer, intent(out), optional :: status 517 517 integer, intent(out), optional :: status
complex :: fvn_c_det 518 518 complex :: fvn_c_det
519 519
complex, allocatable :: wc_a(:,:) 520 520 complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 521 521 integer, allocatable :: ipiv(:)
integer :: info,i 522 522 integer :: info,i
523 523
if (present(status)) status=1 524 524 if (present(status)) status=1
allocate(wc_a(d,d)) 525 525 allocate(wc_a(d,d))
allocate(ipiv(d)) 526 526 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 527 527 wc_a(:,:)=a(:,:)
call cgetrf(d,d,wc_a,d,ipiv,info) 528 528 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 529 529 if (info/= 0) then
if (present(status)) status=0 530 530 if (present(status)) status=0
fvn_c_det=(0.e0,0.e0) 531 531 fvn_c_det=(0.e0,0.e0)
deallocate(ipiv) 532 532 deallocate(ipiv)
deallocate(wc_a) 533 533 deallocate(wc_a)
return 534 534 return
end if 535 535 end if
fvn_c_det=(1.e0,0.e0) 536 536 fvn_c_det=(1.e0,0.e0)
do i=1,d 537 537 do i=1,d
if (ipiv(i)==i) then 538 538 if (ipiv(i)==i) then
fvn_c_det=fvn_c_det*wc_a(i,i) 539 539 fvn_c_det=fvn_c_det*wc_a(i,i)
else 540 540 else
fvn_c_det=-fvn_c_det*wc_a(i,i) 541 541 fvn_c_det=-fvn_c_det*wc_a(i,i)
end if 542 542 end if
end do 543 543 end do
deallocate(ipiv) 544 544 deallocate(ipiv)
deallocate(wc_a) 545 545 deallocate(wc_a)
546 546
end function 547 547 end function
548 548
function fvn_z_det(d,a,status) 549 549 function fvn_z_det(d,a,status)
! 550 550 !
! Evaluate the determinant of a square matrix using lapack LU factorization 551 551 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 552 552 !
! d (in) : matrix rank 553 553 ! d (in) : matrix rank
! a (in) : The Matrix 554 554 ! a (in) : The Matrix
! det (out) : determinant 555 555 ! det (out) : determinant
! status (out) : =0 if LU factorization failed 556 556 ! status (out) : =0 if LU factorization failed
! 557 557 !
implicit none 558 558 implicit none
integer, intent(in) :: d 559 559 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 560 560 double complex, intent(in) :: a(d,d)
integer, intent(out), optional :: status 561 561 integer, intent(out), optional :: status
double complex :: fvn_z_det 562 562 double complex :: fvn_z_det
563 563
double complex, allocatable :: wc_a(:,:) 564 564 double complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 565 565 integer, allocatable :: ipiv(:)
integer :: info,i 566 566 integer :: info,i
567 567
if (present(status)) status=1 568 568 if (present(status)) status=1
allocate(wc_a(d,d)) 569 569 allocate(wc_a(d,d))
allocate(ipiv(d)) 570 570 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 571 571 wc_a(:,:)=a(:,:)
call zgetrf(d,d,wc_a,d,ipiv,info) 572 572 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 573 573 if (info/= 0) then
if (present(status)) status=0 574 574 if (present(status)) status=0
fvn_z_det=(0.d0,0.d0) 575 575 fvn_z_det=(0.d0,0.d0)
deallocate(ipiv) 576 576 deallocate(ipiv)
deallocate(wc_a) 577 577 deallocate(wc_a)
return 578 578 return
end if 579 579 end if
fvn_z_det=(1.d0,0.d0) 580 580 fvn_z_det=(1.d0,0.d0)
do i=1,d 581 581 do i=1,d
if (ipiv(i)==i) then 582 582 if (ipiv(i)==i) then
fvn_z_det=fvn_z_det*wc_a(i,i) 583 583 fvn_z_det=fvn_z_det*wc_a(i,i)
else 584 584 else
fvn_z_det=-fvn_z_det*wc_a(i,i) 585 585 fvn_z_det=-fvn_z_det*wc_a(i,i)
end if 586 586 end if
end do 587 587 end do
deallocate(ipiv) 588 588 deallocate(ipiv)
deallocate(wc_a) 589 589 deallocate(wc_a)
590 590
end function 591 591 end function
592 592
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 593 593 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 594 594 !
! Condition test 595 595 ! Condition test
! 596 596 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 597 597 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1-norm 598 598 ! 1-norm
! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm 599 599 ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond 600 600 ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
! 601 601 !
subroutine fvn_s_matcon(d,a,rcond,status) 602 602 subroutine fvn_s_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 603 603 ! Matrix condition (reciprocal of condition number)
! 604 604 !
! d (in) : matrix rank 605 605 ! d (in) : matrix rank
! a (in) : The Matrix 606 606 ! a (in) : The Matrix
! rcond (out) : guess what 607 607 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 608 608 ! status (out) : =0 if something went wrong
! 609 609 !
implicit none 610 610 implicit none
integer, intent(in) :: d 611 611 integer, intent(in) :: d
real, intent(in) :: a(d,d) 612 612 real, intent(in) :: a(d,d)
real, intent(out) :: rcond 613 613 real, intent(out) :: rcond
integer, intent(out), optional :: status 614 614 integer, intent(out), optional :: status
615 615
real, allocatable :: work(:) 616 616 real, allocatable :: work(:)
integer, allocatable :: iwork(:) 617 617 integer, allocatable :: iwork(:)
real :: anorm 618 618 real :: anorm
real, allocatable :: wc_a(:,:) ! working copy of a 619 619 real, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 620 620 integer :: info
integer, allocatable :: ipiv(:) 621 621 integer, allocatable :: ipiv(:)
622 622
real, external :: slange 623 623 real, external :: slange
624 624
625 625
if (present(status)) status=1 626 626 if (present(status)) status=1
627 627
anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 628 628 anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
629 629
allocate(wc_a(d,d)) 630 630 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 631 631 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 632 632 wc_a(:,:)=a(:,:)
633 633
allocate(ipiv(d)) 634 634 allocate(ipiv(d))
call sgetrf(d,d,wc_a,d,ipiv,info) 635 635 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 636 636 if (info /= 0) then
if (present(status)) status=0 637 637 if (present(status)) status=0
deallocate(ipiv) 638 638 deallocate(ipiv)
deallocate(wc_a) 639 639 deallocate(wc_a)
return 640 640 return
end if 641 641 end if
allocate(work(4*d)) 642 642 allocate(work(4*d))
allocate(iwork(d)) 643 643 allocate(iwork(d))
call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 644 644 call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 645 645 if (info /= 0) then
if (present(status)) status=0 646 646 if (present(status)) status=0
end if 647 647 end if
deallocate(iwork) 648 648 deallocate(iwork)
deallocate(work) 649 649 deallocate(work)
deallocate(ipiv) 650 650 deallocate(ipiv)
deallocate(wc_a) 651 651 deallocate(wc_a)
652 652
end subroutine 653 653 end subroutine
654 654
subroutine fvn_d_matcon(d,a,rcond,status) 655 655 subroutine fvn_d_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 656 656 ! Matrix condition (reciprocal of condition number)
! 657 657 !
! d (in) : matrix rank 658 658 ! d (in) : matrix rank
! a (in) : The Matrix 659 659 ! a (in) : The Matrix
! rcond (out) : guess what 660 660 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 661 661 ! status (out) : =0 if something went wrong
! 662 662 !
implicit none 663 663 implicit none
integer, intent(in) :: d 664 664 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 665 665 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 666 666 double precision, intent(out) :: rcond
integer, intent(out), optional :: status 667 667 integer, intent(out), optional :: status
668 668
double precision, allocatable :: work(:) 669 669 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 670 670 integer, allocatable :: iwork(:)
double precision :: anorm 671 671 double precision :: anorm
double precision, allocatable :: wc_a(:,:) ! working copy of a 672 672 double precision, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 673 673 integer :: info
integer, allocatable :: ipiv(:) 674 674 integer, allocatable :: ipiv(:)
675 675
double precision, external :: dlange 676 676 double precision, external :: dlange
677 677
678 678
if (present(status)) status=1 679 679 if (present(status)) status=1
680 680
anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 681 681 anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
682 682
allocate(wc_a(d,d)) 683 683 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 684 684 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 685 685 wc_a(:,:)=a(:,:)
686 686
allocate(ipiv(d)) 687 687 allocate(ipiv(d))
call dgetrf(d,d,wc_a,d,ipiv,info) 688 688 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 689 689 if (info /= 0) then
if (present(status)) status=0 690 690 if (present(status)) status=0
deallocate(ipiv) 691 691 deallocate(ipiv)
deallocate(wc_a) 692 692 deallocate(wc_a)
return 693 693 return
end if 694 694 end if
695 695
allocate(work(4*d)) 696 696 allocate(work(4*d))
allocate(iwork(d)) 697 697 allocate(iwork(d))
call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 698 698 call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 699 699 if (info /= 0) then
if (present(status)) status=0 700 700 if (present(status)) status=0
end if 701 701 end if
deallocate(iwork) 702 702 deallocate(iwork)
deallocate(work) 703 703 deallocate(work)
deallocate(ipiv) 704 704 deallocate(ipiv)
deallocate(wc_a) 705 705 deallocate(wc_a)
706 706
end subroutine 707 707 end subroutine
708 708
subroutine fvn_c_matcon(d,a,rcond,status) 709 709 subroutine fvn_c_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 710 710 ! Matrix condition (reciprocal of condition number)
! 711 711 !
! d (in) : matrix rank 712 712 ! d (in) : matrix rank
! a (in) : The Matrix 713 713 ! a (in) : The Matrix
! rcond (out) : guess what 714 714 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 715 715 ! status (out) : =0 if something went wrong
! 716 716 !
implicit none 717 717 implicit none
integer, intent(in) :: d 718 718 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 719 719 complex, intent(in) :: a(d,d)
real, intent(out) :: rcond 720 720 real, intent(out) :: rcond
integer, intent(out), optional :: status 721 721 integer, intent(out), optional :: status
722 722
real, allocatable :: rwork(:) 723 723 real, allocatable :: rwork(:)
complex, allocatable :: work(:) 724 724 complex, allocatable :: work(:)
integer, allocatable :: iwork(:) 725 725 integer, allocatable :: iwork(:)
real :: anorm 726 726 real :: anorm
complex, allocatable :: wc_a(:,:) ! working copy of a 727 727 complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 728 728 integer :: info
integer, allocatable :: ipiv(:) 729 729 integer, allocatable :: ipiv(:)
730 730
real, external :: clange 731 731 real, external :: clange
732 732
733 733
if (present(status)) status=1 734 734 if (present(status)) status=1
735 735
anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 736 736 anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
737 737
allocate(wc_a(d,d)) 738 738 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 739 739 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 740 740 wc_a(:,:)=a(:,:)
741 741
allocate(ipiv(d)) 742 742 allocate(ipiv(d))
call cgetrf(d,d,wc_a,d,ipiv,info) 743 743 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 744 744 if (info /= 0) then
if (present(status)) status=0 745 745 if (present(status)) status=0
deallocate(ipiv) 746 746 deallocate(ipiv)
deallocate(wc_a) 747 747 deallocate(wc_a)
return 748 748 return
end if 749 749 end if
allocate(work(2*d)) 750 750 allocate(work(2*d))
allocate(rwork(2*d)) 751 751 allocate(rwork(2*d))
call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 752 752 call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 753 753 if (info /= 0) then
if (present(status)) status=0 754 754 if (present(status)) status=0
end if 755 755 end if
deallocate(rwork) 756 756 deallocate(rwork)
deallocate(work) 757 757 deallocate(work)
deallocate(ipiv) 758 758 deallocate(ipiv)
deallocate(wc_a) 759 759 deallocate(wc_a)
end subroutine 760 760 end subroutine
761 761
subroutine fvn_z_matcon(d,a,rcond,status) 762 762 subroutine fvn_z_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 763 763 ! Matrix condition (reciprocal of condition number)
! 764 764 !
! d (in) : matrix rank 765 765 ! d (in) : matrix rank
! a (in) : The Matrix 766 766 ! a (in) : The Matrix
! rcond (out) : guess what 767 767 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 768 768 ! status (out) : =0 if something went wrong
! 769 769 !
implicit none 770 770 implicit none
integer, intent(in) :: d 771 771 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 772 772 double complex, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 773 773 double precision, intent(out) :: rcond
integer, intent(out), optional :: status 774 774 integer, intent(out), optional :: status
775 775
double complex, allocatable :: work(:) 776 776 double complex, allocatable :: work(:)
double precision, allocatable :: rwork(:) 777 777 double precision, allocatable :: rwork(:)
double precision :: anorm 778 778 double precision :: anorm
double complex, allocatable :: wc_a(:,:) ! working copy of a 779 779 double complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 780 780 integer :: info
integer, allocatable :: ipiv(:) 781 781 integer, allocatable :: ipiv(:)
782 782
double precision, external :: zlange 783 783 double precision, external :: zlange
784 784
785 785
if (present(status)) status=1 786 786 if (present(status)) status=1
787 787
anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 788 788 anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
789 789
allocate(wc_a(d,d)) 790 790 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 791 791 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 792 792 wc_a(:,:)=a(:,:)
793 793
allocate(ipiv(d)) 794 794 allocate(ipiv(d))
call zgetrf(d,d,wc_a,d,ipiv,info) 795 795 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 796 796 if (info /= 0) then
if (present(status)) status=0 797 797 if (present(status)) status=0
deallocate(ipiv) 798 798 deallocate(ipiv)
deallocate(wc_a) 799 799 deallocate(wc_a)
return 800 800 return
end if 801 801 end if
802 802
allocate(work(2*d)) 803 803 allocate(work(2*d))
allocate(rwork(2*d)) 804 804 allocate(rwork(2*d))
call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 805 805 call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 806 806 if (info /= 0) then
if (present(status)) status=0 807 807 if (present(status)) status=0
end if 808 808 end if
deallocate(rwork) 809 809 deallocate(rwork)
deallocate(work) 810 810 deallocate(work)
deallocate(ipiv) 811 811 deallocate(ipiv)
deallocate(wc_a) 812 812 deallocate(wc_a)
end subroutine 813 813 end subroutine
814 814
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 815 815 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 816 816 !
! Valeurs propres/ Vecteurs propre 817 817 ! Valeurs propres/ Vecteurs propre
! 818 818 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 819 819 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
820 820
subroutine fvn_s_matev(d,a,evala,eveca,status) 821 821 subroutine fvn_s_matev(d,a,evala,eveca,status)
! 822 822 !
! integer d (in) : matrice rank 823 823 ! integer d (in) : matrice rank
! real a(d,d) (in) : The Matrix 824 824 ! real a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 825 825 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 826 826 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 827 827 ! integer (out) : status =0 if something went wrong
! 828 828 !
! interfacing Lapack routine SGEEV 829 829 ! interfacing Lapack routine SGEEV
implicit none 830 830 implicit none
integer, intent(in) :: d 831 831 integer, intent(in) :: d
real, intent(in) :: a(d,d) 832 832 real, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 833 833 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 834 834 complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 835 835 integer, intent(out), optional :: status
836 836
real, allocatable :: wc_a(:,:) ! a working copy of a 837 837 real, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 838 838 integer :: info
integer :: lwork 839 839 integer :: lwork
real, allocatable :: wr(:),wi(:) 840 840 real, allocatable :: wr(:),wi(:)
real :: vl ! unused but necessary for the call 841 841 real :: vl ! unused but necessary for the call
real, allocatable :: vr(:,:) 842 842 real, allocatable :: vr(:,:)
real, allocatable :: work(:) 843 843 real, allocatable :: work(:)
real :: twork(1) 844 844 real :: twork(1)
integer i 845 845 integer i
integer j 846 846 integer j
847 847
if (present(status)) status=1 848 848 if (present(status)) status=1
849 849
! making a working copy of a 850 850 ! making a working copy of a
allocate(wc_a(d,d)) 851 851 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 852 852 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 853 853 wc_a(:,:)=a(:,:)
854 854
allocate(wr(d)) 855 855 allocate(wr(d))
allocate(wi(d)) 856 856 allocate(wi(d))
allocate(vr(d,d)) 857 857 allocate(vr(d,d))
! query optimal work size 858 858 ! query optimal work size
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 859 859 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 860 860 lwork=int(twork(1))
allocate(work(lwork)) 861 861 allocate(work(lwork))
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 862 862 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
863 863
if (info /= 0) then 864 864 if (info /= 0) then
if (present(status)) status=0 865 865 if (present(status)) status=0
deallocate(work) 866 866 deallocate(work)
deallocate(vr) 867 867 deallocate(vr)
deallocate(wi) 868 868 deallocate(wi)
deallocate(wr) 869 869 deallocate(wr)
deallocate(wc_a) 870 870 deallocate(wc_a)
return 871 871 return
end if 872 872 end if
873 873
! now fill in the results 874 874 ! now fill in the results
i=1 875 875 i=1
do while(i<=d) 876 876 do while(i<=d)
evala(i)=cmplx(wr(i),wi(i)) 877 877 evala(i)=cmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 878 878 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0.) 879 879 eveca(:,i)=cmplx(vr(:,i),0.)
else ! eigenvalue is complex 880 880 else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1)) 881 881 evala(i+1)=cmplx(wr(i+1),wi(i+1))
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1)) 882 882 eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1)) 883 883 eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
i=i+1 884 884 i=i+1
end if 885 885 end if
i=i+1 886 886 i=i+1
enddo 887 887 enddo
deallocate(work) 888 888 deallocate(work)
deallocate(vr) 889 889 deallocate(vr)
deallocate(wi) 890 890 deallocate(wi)
deallocate(wr) 891 891 deallocate(wr)
deallocate(wc_a) 892 892 deallocate(wc_a)
893 893
end subroutine 894 894 end subroutine
895 895
subroutine fvn_d_matev(d,a,evala,eveca,status) 896 896 subroutine fvn_d_matev(d,a,evala,eveca,status)
! 897 897 !
! integer d (in) : matrice rank 898 898 ! integer d (in) : matrice rank
! double precision a(d,d) (in) : The Matrix 899 899 ! double precision a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 900 900 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 901 901 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 902 902 ! integer (out) : status =0 if something went wrong
! 903 903 !
! interfacing Lapack routine DGEEV 904 904 ! interfacing Lapack routine DGEEV
implicit none 905 905 implicit none
integer, intent(in) :: d 906 906 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 907 907 double precision, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 908 908 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 909 909 double complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 910 910 integer, intent(out), optional :: status
911 911
double precision, allocatable :: wc_a(:,:) ! a working copy of a 912 912 double precision, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 913 913 integer :: info
integer :: lwork 914 914 integer :: lwork
double precision, allocatable :: wr(:),wi(:) 915 915 double precision, allocatable :: wr(:),wi(:)
double precision :: vl ! unused but necessary for the call 916 916 double precision :: vl ! unused but necessary for the call
double precision, allocatable :: vr(:,:) 917 917 double precision, allocatable :: vr(:,:)
double precision, allocatable :: work(:) 918 918 double precision, allocatable :: work(:)
double precision :: twork(1) 919 919 double precision :: twork(1)
integer i 920 920 integer i
integer j 921 921 integer j
922 922
if (present(status)) status=1 923 923 if (present(status)) status=1
924 924
! making a working copy of a 925 925 ! making a working copy of a
allocate(wc_a(d,d)) 926 926 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 927 927 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 928 928 wc_a(:,:)=a(:,:)
929 929
allocate(wr(d)) 930 930 allocate(wr(d))
allocate(wi(d)) 931 931 allocate(wi(d))
allocate(vr(d,d)) 932 932 allocate(vr(d,d))
! query optimal work size 933 933 ! query optimal work size
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 934 934 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 935 935 lwork=int(twork(1))
allocate(work(lwork)) 936 936 allocate(work(lwork))
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 937 937 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
938 938
if (info /= 0) then 939 939 if (info /= 0) then
if (present(status)) status=0 940 940 if (present(status)) status=0
deallocate(work) 941 941 deallocate(work)
deallocate(vr) 942 942 deallocate(vr)
deallocate(wi) 943 943 deallocate(wi)
deallocate(wr) 944 944 deallocate(wr)
deallocate(wc_a) 945 945 deallocate(wc_a)
return 946 946 return
end if 947 947 end if
948 948
! now fill in the results 949 949 ! now fill in the results
i=1 950 950 i=1
do while(i<=d) 951 951 do while(i<=d)
evala(i)=dcmplx(wr(i),wi(i)) 952 952 evala(i)=dcmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 953 953 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=dcmplx(vr(:,i),0.) 954 954 eveca(:,i)=dcmplx(vr(:,i),0.)
else ! eigenvalue is complex 955 955 else ! eigenvalue is complex
evala(i+1)=dcmplx(wr(i+1),wi(i+1)) 956 956 evala(i+1)=dcmplx(wr(i+1),wi(i+1))
eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1)) 957 957 eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1)) 958 958 eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
i=i+1 959 959 i=i+1
end if 960 960 end if
i=i+1 961 961 i=i+1
enddo 962 962 enddo
963 963
deallocate(work) 964 964 deallocate(work)
deallocate(vr) 965 965 deallocate(vr)
deallocate(wi) 966 966 deallocate(wi)
deallocate(wr) 967 967 deallocate(wr)
deallocate(wc_a) 968 968 deallocate(wc_a)
969 969
end subroutine 970 970 end subroutine
971 971
subroutine fvn_c_matev(d,a,evala,eveca,status) 972 972 subroutine fvn_c_matev(d,a,evala,eveca,status)
! 973 973 !
! integer d (in) : matrice rank 974 974 ! integer d (in) : matrice rank
! complex a(d,d) (in) : The Matrix 975 975 ! complex a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 976 976 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 977 977 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 978 978 ! integer (out) : status =0 if something went wrong
! 979 979 !
! interfacing Lapack routine CGEEV 980 980 ! interfacing Lapack routine CGEEV
implicit none 981 981 implicit none
integer, intent(in) :: d 982 982 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 983 983 complex, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 984 984 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 985 985 complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 986 986 integer, intent(out), optional :: status
987 987
complex, allocatable :: wc_a(:,:) ! a working copy of a 988 988 complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 989 989 integer :: info
integer :: lwork 990 990 integer :: lwork
complex, allocatable :: work(:) 991 991 complex, allocatable :: work(:)
complex :: twork(1) 992 992 complex :: twork(1)
real, allocatable :: rwork(:) 993 993 real, allocatable :: rwork(:)
complex :: vl ! unused but necessary for the call 994 994 complex :: vl ! unused but necessary for the call
995 995
if (present(status)) status=1 996 996 if (present(status)) status=1
997 997
! making a working copy of a 998 998 ! making a working copy of a
allocate(wc_a(d,d)) 999 999 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 1000 1000 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1001 1001 wc_a(:,:)=a(:,:)
1002 1002
1003 1003
! query optimal work size 1004 1004 ! query optimal work size
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1005 1005 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 1006 1006 lwork=int(twork(1))
allocate(work(lwork)) 1007 1007 allocate(work(lwork))
allocate(rwork(2*d)) 1008 1008 allocate(rwork(2*d))
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1009 1009 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1010 1010
if (info /= 0) then 1011 1011 if (info /= 0) then
if (present(status)) status=0 1012 1012 if (present(status)) status=0
end if 1013 1013 end if
deallocate(rwork) 1014 1014 deallocate(rwork)
deallocate(work) 1015 1015 deallocate(work)
deallocate(wc_a) 1016 1016 deallocate(wc_a)
1017 1017
end subroutine 1018 1018 end subroutine
1019 1019
subroutine fvn_z_matev(d,a,evala,eveca,status) 1020 1020 subroutine fvn_z_matev(d,a,evala,eveca,status)
! 1021 1021 !
! integer d (in) : matrice rank 1022 1022 ! integer d (in) : matrice rank
! double complex a(d,d) (in) : The Matrix 1023 1023 ! double complex a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 1024 1024 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 1025 1025 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 1026 1026 ! integer (out) : status =0 if something went wrong
! 1027 1027 !
! interfacing Lapack routine ZGEEV 1028 1028 ! interfacing Lapack routine ZGEEV
implicit none 1029 1029 implicit none
integer, intent(in) :: d 1030 1030 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 1031 1031 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 1032 1032 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 1033 1033 double complex, intent(out) :: eveca(d,d)
integer, intent(out), optional :: status 1034 1034 integer, intent(out), optional :: status
1035 1035
double complex, allocatable :: wc_a(:,:) ! a working copy of a 1036 1036 double complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 1037 1037 integer :: info
integer :: lwork 1038 1038 integer :: lwork
double complex, allocatable :: work(:) 1039 1039 double complex, allocatable :: work(:)
double complex :: twork(1) 1040 1040 double complex :: twork(1)
double precision, allocatable :: rwork(:) 1041 1041 double precision, allocatable :: rwork(:)
double complex :: vl ! unused but necessary for the call 1042 1042 double complex :: vl ! unused but necessary for the call
1043 1043
if (present(status)) status=1 1044 1044 if (present(status)) status=1
1045 1045
! making a working copy of a 1046 1046 ! making a working copy of a
allocate(wc_a(d,d)) 1047 1047 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 1048 1048 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 1049 1049 wc_a(:,:)=a(:,:)
1050 1050
! query optimal work size 1051 1051 ! query optimal work size
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 1052 1052 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 1053 1053 lwork=int(twork(1))
allocate(work(lwork)) 1054 1054 allocate(work(lwork))
allocate(rwork(2*d)) 1055 1055 allocate(rwork(2*d))
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 1056 1056 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
1057 1057
if (info /= 0) then 1058 1058 if (info /= 0) then
if (present(status)) status=0 1059 1059 if (present(status)) status=0
end if 1060 1060 end if
deallocate(rwork) 1061 1061 deallocate(rwork)
deallocate(work) 1062 1062 deallocate(work)
deallocate(wc_a) 1063 1063 deallocate(wc_a)
1064 1064
end subroutine 1065 1065 end subroutine
1066 1066
1067 1067
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1068 1068 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1069 1069 !
! Quadratic interpolation of tabulated function of 1,2 or 3 variables 1070 1070 ! Quadratic interpolation of tabulated function of 1,2 or 3 variables
! 1071 1071 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1072 1072 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1073 1073
subroutine fvn_s_find_interval(x,i,xdata,n) 1074 1074 subroutine fvn_s_find_interval(x,i,xdata,n)
implicit none 1075 1075 implicit none
! This routine find the indice i where xdata(i) <= x < xdata(i+1) 1076 1076 ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
! xdata(n) must contains a set of increasingly ordered values 1077 1077 ! xdata(n) must contains a set of increasingly ordered values
! if x < xdata(1) i=0 is returned 1078 1078 ! if x < xdata(1) i=0 is returned
! if x > xdata(n) i=n is returned 1079 1079 ! if x > xdata(n) i=n is returned
! special case is where x=xdata(n) then n-1 is returned so 1080 1080 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 1081 1081 ! we will not exclude the upper limit
! a simple dichotomy method is used 1082 1082 ! a simple dichotomy method is used
1083 1083
real(kind=4), intent(in) :: x 1084 1084 real(kind=4), intent(in) :: x
real(kind=4), intent(in), dimension(n) :: xdata 1085 1085 real(kind=4), intent(in), dimension(n) :: xdata
integer(kind=4), intent(in) :: n 1086 1086 integer(kind=4), intent(in) :: n
integer(kind=4), intent(out) :: i 1087 1087 integer(kind=4), intent(out) :: i
1088 1088
integer(kind=4) :: imin,imax,imoyen 1089 1089 integer(kind=4) :: imin,imax,imoyen
1090 1090
! special case is where x=xdata(n) then n-1 is returned so 1091 1091 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 1092 1092 ! we will not exclude the upper limit
if (x == xdata(n)) then 1093 1093 if (x == xdata(n)) then
i=n-1 1094 1094 i=n-1
return 1095 1095 return
end if 1096 1096 end if
1097 1097
! if x < xdata(1) i=0 is returned 1098 1098 ! if x < xdata(1) i=0 is returned
if (x < xdata(1)) then 1099 1099 if (x < xdata(1)) then
i=0 1100 1100 i=0
return 1101 1101 return
end if 1102 1102 end if
1103 1103
! if x > xdata(n) i=n is returned 1104 1104 ! if x > xdata(n) i=n is returned
if (x > xdata(n)) then 1105 1105 if (x > xdata(n)) then
i=n 1106 1106 i=n
return 1107 1107 return
end if 1108 1108 end if
1109 1109
! here xdata(1) <= x <= xdata(n) 1110 1110 ! here xdata(1) <= x <= xdata(n)
imin=0 1111 1111 imin=0
imax=n+1 1112 1112 imax=n+1
1113 1113
do while((imax-imin) > 1) 1114 1114 do while((imax-imin) > 1)
imoyen=(imax+imin)/2 1115 1115 imoyen=(imax+imin)/2
if (x >= xdata(imoyen)) then 1116 1116 if (x >= xdata(imoyen)) then
imin=imoyen 1117 1117 imin=imoyen
else 1118 1118 else
imax=imoyen 1119 1119 imax=imoyen
end if 1120 1120 end if
end do 1121 1121 end do
1122 1122
i=imin 1123 1123 i=imin
1124 1124
end subroutine 1125 1125 end subroutine
1126 1126
1127 1127
subroutine fvn_d_find_interval(x,i,xdata,n) 1128 1128 subroutine fvn_d_find_interval(x,i,xdata,n)
implicit none 1129 1129 implicit none
! This routine find the indice i where xdata(i) <= x < xdata(i+1) 1130 1130 ! This routine find the indice i where xdata(i) <= x < xdata(i+1)
! xdata(n) must contains a set of increasingly ordered values 1131 1131 ! xdata(n) must contains a set of increasingly ordered values
! if x < xdata(1) i=0 is returned 1132 1132 ! if x < xdata(1) i=0 is returned
! if x > xdata(n) i=n is returned 1133 1133 ! if x > xdata(n) i=n is returned
! special case is where x=xdata(n) then n-1 is returned so 1134 1134 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 1135 1135 ! we will not exclude the upper limit
! a simple dichotomy method is used 1136 1136 ! a simple dichotomy method is used
1137 1137
real(kind=8), intent(in) :: x 1138 1138 real(kind=8), intent(in) :: x
real(kind=8), intent(in), dimension(n) :: xdata 1139 1139 real(kind=8), intent(in), dimension(n) :: xdata
integer(kind=4), intent(in) :: n 1140 1140 integer(kind=4), intent(in) :: n
integer(kind=4), intent(out) :: i 1141 1141 integer(kind=4), intent(out) :: i
1142 1142
integer(kind=4) :: imin,imax,imoyen 1143 1143 integer(kind=4) :: imin,imax,imoyen
1144 1144
! special case is where x=xdata(n) then n-1 is returned so 1145 1145 ! special case is where x=xdata(n) then n-1 is returned so
! we will not exclude the upper limit 1146 1146 ! we will not exclude the upper limit
if (x == xdata(n)) then 1147 1147 if (x == xdata(n)) then
i=n-1 1148 1148 i=n-1
return 1149 1149 return
end if 1150 1150 end if
1151 1151
! if x < xdata(1) i=0 is returned 1152 1152 ! if x < xdata(1) i=0 is returned
if (x < xdata(1)) then 1153 1153 if (x < xdata(1)) then
i=0 1154 1154 i=0
return 1155 1155 return
end if 1156 1156 end if
1157 1157
! if x > xdata(n) i=n is returned 1158 1158 ! if x > xdata(n) i=n is returned
if (x > xdata(n)) then 1159 1159 if (x > xdata(n)) then
i=n 1160 1160 i=n
return 1161 1161 return
end if 1162 1162 end if
1163 1163
! here xdata(1) <= x <= xdata(n) 1164 1164 ! here xdata(1) <= x <= xdata(n)
imin=0 1165 1165 imin=0
imax=n+1 1166 1166 imax=n+1
1167 1167
do while((imax-imin) > 1) 1168 1168 do while((imax-imin) > 1)
imoyen=(imax+imin)/2 1169 1169 imoyen=(imax+imin)/2
if (x >= xdata(imoyen)) then 1170 1170 if (x >= xdata(imoyen)) then
imin=imoyen 1171 1171 imin=imoyen
else 1172 1172 else
imax=imoyen 1173 1173 imax=imoyen
end if 1174 1174 end if
end do 1175 1175 end do
1176 1176
i=imin 1177 1177 i=imin
1178 1178
end subroutine 1179 1179 end subroutine
1180 1180
1181 1181
function fvn_s_quad_interpol(x,n,xdata,ydata) 1182 1182 function fvn_s_quad_interpol(x,n,xdata,ydata)
implicit none 1183 1183 implicit none
! This function evaluate the value of a function defined by a set of points 1184 1184 ! This function evaluate the value of a function defined by a set of points
! and values, using a quadratic interpolation 1185 1185 ! and values, using a quadratic interpolation
! xdata must be increasingly ordered 1186 1186 ! xdata must be increasingly ordered
! x must be within xdata(1) and xdata(n) to actually do interpolation 1187 1187 ! x must be within xdata(1) and xdata(n) to actually do interpolation
! otherwise extrapolation is done 1188 1188 ! otherwise extrapolation is done
integer(kind=4), intent(in) :: n 1189 1189 integer(kind=4), intent(in) :: n
real(kind=4), intent(in), dimension(n) :: xdata,ydata 1190 1190 real(kind=4), intent(in), dimension(n) :: xdata,ydata
real(kind=4), intent(in) :: x 1191 1191 real(kind=4), intent(in) :: x
real(kind=4) :: fvn_s_quad_interpol 1192 1192 real(kind=4) :: fvn_s_quad_interpol
1193 1193
integer(kind=4) :: iinf,base,i,j 1194 1194 integer(kind=4) :: iinf,base,i,j
real(kind=4) :: p 1195 1195 real(kind=4) :: p
1196 1196
call fvn_s_find_interval(x,iinf,xdata,n) 1197 1197 call fvn_s_find_interval(x,iinf,xdata,n)
1198 1198
! Settings for extrapolation 1199 1199 ! Settings for extrapolation
if (iinf==0) then 1200 1200 if (iinf==0) then
! TODO -> Lower bound extrapolation warning 1201 1201 ! TODO -> Lower bound extrapolation warning
iinf=1 1202 1202 iinf=1
end if 1203 1203 end if
1204 1204
if (iinf==n) then 1205 1205 if (iinf==n) then
! TODO -> Higher bound extrapolation warning 1206 1206 ! TODO -> Higher bound extrapolation warning
iinf=n-1 1207 1207 iinf=n-1
end if 1208 1208 end if
1209 1209
! The three points we will use are iinf-1,iinf and iinf+1 with the 1210 1210 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 1211 1211 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (iinf==1) then 1212 1212 if (iinf==1) then
base=0 1213 1213 base=0
else 1214 1214 else
base=iinf-2 1215 1215 base=iinf-2
end if 1216 1216 end if
1217 1217
! The three points we will use are : 1218 1218 ! The three points we will use are :
! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3) 1219 1219 ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
1220 1220
! Straight forward Lagrange polynomial 1221 1221 ! Straight forward Lagrange polynomial
fvn_s_quad_interpol=0. 1222 1222 fvn_s_quad_interpol=0.
do i=1,3 1223 1223 do i=1,3
! polynome i 1224 1224 ! polynome i
p=ydata(base+i) 1225 1225 p=ydata(base+i)
do j=1,3 1226 1226 do j=1,3
if (j /= i) then 1227 1227 if (j /= i) then
p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j)) 1228 1228 p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
end if 1229 1229 end if
end do 1230 1230 end do
fvn_s_quad_interpol=fvn_s_quad_interpol+p 1231 1231 fvn_s_quad_interpol=fvn_s_quad_interpol+p
end do 1232 1232 end do
1233 1233
end function 1234 1234 end function
1235 1235
1236 1236
function fvn_d_quad_interpol(x,n,xdata,ydata) 1237 1237 function fvn_d_quad_interpol(x,n,xdata,ydata)
implicit none 1238 1238 implicit none
! This function evaluate the value of a function defined by a set of points 1239 1239 ! This function evaluate the value of a function defined by a set of points
! and values, using a quadratic interpolation 1240 1240 ! and values, using a quadratic interpolation
! xdata must be increasingly ordered 1241 1241 ! xdata must be increasingly ordered
! x must be within xdata(1) and xdata(n) to actually do interpolation 1242 1242 ! x must be within xdata(1) and xdata(n) to actually do interpolation
! otherwise extrapolation is done 1243 1243 ! otherwise extrapolation is done
integer(kind=4), intent(in) :: n 1244 1244 integer(kind=4), intent(in) :: n
real(kind=8), intent(in), dimension(n) :: xdata,ydata 1245 1245 real(kind=8), intent(in), dimension(n) :: xdata,ydata
real(kind=8), intent(in) :: x 1246 1246 real(kind=8), intent(in) :: x
real(kind=8) :: fvn_d_quad_interpol 1247 1247 real(kind=8) :: fvn_d_quad_interpol
1248 1248
integer(kind=4) :: iinf,base,i,j 1249 1249 integer(kind=4) :: iinf,base,i,j
real(kind=8) :: p 1250 1250 real(kind=8) :: p
1251 1251
call fvn_d_find_interval(x,iinf,xdata,n) 1252 1252 call fvn_d_find_interval(x,iinf,xdata,n)
1253 1253
! Settings for extrapolation 1254 1254 ! Settings for extrapolation
if (iinf==0) then 1255 1255 if (iinf==0) then
! TODO -> Lower bound extrapolation warning 1256 1256 ! TODO -> Lower bound extrapolation warning
iinf=1 1257 1257 iinf=1
end if 1258 1258 end if
1259 1259
if (iinf==n) then 1260 1260 if (iinf==n) then
! TODO Higher bound extrapolation warning 1261 1261 ! TODO Higher bound extrapolation warning
iinf=n-1 1262 1262 iinf=n-1
end if 1263 1263 end if
1264 1264
! The three points we will use are iinf-1,iinf and iinf+1 with the 1265 1265 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 1266 1266 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (iinf==1) then 1267 1267 if (iinf==1) then
base=0 1268 1268 base=0
else 1269 1269 else
base=iinf-2 1270 1270 base=iinf-2
end if 1271 1271 end if
1272 1272
! The three points we will use are : 1273 1273 ! The three points we will use are :
! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3) 1274 1274 ! xdata/ydata(base+1),xdata/ydata(base+2),xdata/ydata(base+3)
1275 1275
! Straight forward Lagrange polynomial 1276 1276 ! Straight forward Lagrange polynomial
fvn_d_quad_interpol=0. 1277 1277 fvn_d_quad_interpol=0.
do i=1,3 1278 1278 do i=1,3
! polynome i 1279 1279 ! polynome i
p=ydata(base+i) 1280 1280 p=ydata(base+i)
do j=1,3 1281 1281 do j=1,3
if (j /= i) then 1282 1282 if (j /= i) then
p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j)) 1283 1283 p=p*(x-xdata(base+j))/(xdata(base+i)-xdata(base+j))
end if 1284 1284 end if
end do 1285 1285 end do
fvn_d_quad_interpol=fvn_d_quad_interpol+p 1286 1286 fvn_d_quad_interpol=fvn_d_quad_interpol+p
end do 1287 1287 end do
1288 1288
end function 1289 1289 end function
1290 1290
1291 1291
function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata) 1292 1292 function fvn_s_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
implicit none 1293 1293 implicit none
! This function evaluate the value of a two variable function defined by a 1294 1294 ! This function evaluate the value of a two variable function defined by a
! set of points and values, using a quadratic interpolation 1295 1295 ! set of points and values, using a quadratic interpolation
! xdata and ydata must be increasingly ordered 1296 1296 ! xdata and ydata must be increasingly ordered
! the couple (x,y) must be as x within xdata(1) and xdata(nx) and 1297 1297 ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
! y within ydata(1) and ydata(ny) to actually do interpolation 1298 1298 ! y within ydata(1) and ydata(ny) to actually do interpolation
! otherwise extrapolation is done 1299 1299 ! otherwise extrapolation is done
integer(kind=4), intent(in) :: nx,ny 1300 1300 integer(kind=4), intent(in) :: nx,ny
real(kind=4), intent(in) :: x,y 1301 1301 real(kind=4), intent(in) :: x,y
real(kind=4), intent(in), dimension(nx) :: xdata 1302 1302 real(kind=4), intent(in), dimension(nx) :: xdata
real(kind=4), intent(in), dimension(ny) :: ydata 1303 1303 real(kind=4), intent(in), dimension(ny) :: ydata
real(kind=4), intent(in), dimension(nx,ny) :: zdata 1304 1304 real(kind=4), intent(in), dimension(nx,ny) :: zdata
real(kind=4) :: fvn_s_quad_2d_interpol 1305 1305 real(kind=4) :: fvn_s_quad_2d_interpol
1306 1306
integer(kind=4) :: ixinf,iyinf,basex,basey,i 1307 1307 integer(kind=4) :: ixinf,iyinf,basex,basey,i
real(kind=4),dimension(3) :: ztmp 1308 1308 real(kind=4),dimension(3) :: ztmp
!real(kind=4), external :: fvn_s_quad_interpol 1309 1309 !real(kind=4), external :: fvn_s_quad_interpol
1310 1310
call fvn_s_find_interval(x,ixinf,xdata,nx) 1311 1311 call fvn_s_find_interval(x,ixinf,xdata,nx)
call fvn_s_find_interval(y,iyinf,ydata,ny) 1312 1312 call fvn_s_find_interval(y,iyinf,ydata,ny)
1313 1313
! Settings for extrapolation 1314 1314 ! Settings for extrapolation
if (ixinf==0) then 1315 1315 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 1316 1316 ! TODO -> Lower x bound extrapolation warning
ixinf=1 1317 1317 ixinf=1
end if 1318 1318 end if
1319 1319
if (ixinf==nx) then 1320 1320 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 1321 1321 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 1322 1322 ixinf=nx-1
end if 1323 1323 end if
1324 1324
if (iyinf==0) then 1325 1325 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 1326 1326 ! TODO -> Lower y bound extrapolation warning
iyinf=1 1327 1327 iyinf=1
end if 1328 1328 end if
1329 1329
if (iyinf==ny) then 1330 1330 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 1331 1331 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 1332 1332 iyinf=ny-1
end if 1333 1333 end if
1334 1334
! The three points we will use are iinf-1,iinf and iinf+1 with the 1335 1335 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 1336 1336 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 1337 1337 if (ixinf==1) then
basex=0 1338 1338 basex=0
else 1339 1339 else
basex=ixinf-2 1340 1340 basex=ixinf-2
end if 1341 1341 end if
1342 1342
if (iyinf==1) then 1343 1343 if (iyinf==1) then
basey=0 1344 1344 basey=0
else 1345 1345 else
basey=iyinf-2 1346 1346 basey=iyinf-2
end if 1347 1347 end if
1348 1348
! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3) 1349 1349 ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
! stored in ztmp(1:3) 1350 1350 ! stored in ztmp(1:3)
do i=1,3 1351 1351 do i=1,3
ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i)) 1352 1352 ztmp(i)=fvn_s_quad_interpol(x,nx,xdata,zdata(:,basey+i))
end do 1353 1353 end do
1354 1354
! Then we make an interpolation for y using previous interpolations 1355 1355 ! Then we make an interpolation for y using previous interpolations
fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp) 1356 1356 fvn_s_quad_2d_interpol=fvn_s_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
end function 1357 1357 end function
1358 1358
1359 1359
function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata) 1360 1360 function fvn_d_quad_2d_interpol(x,y,nx,xdata,ny,ydata,zdata)
implicit none 1361 1361 implicit none
! This function evaluate the value of a two variable function defined by a 1362 1362 ! This function evaluate the value of a two variable function defined by a
! set of points and values, using a quadratic interpolation 1363 1363 ! set of points and values, using a quadratic interpolation
! xdata and ydata must be increasingly ordered 1364 1364 ! xdata and ydata must be increasingly ordered
! the couple (x,y) must be as x within xdata(1) and xdata(nx) and 1365 1365 ! the couple (x,y) must be as x within xdata(1) and xdata(nx) and
! y within ydata(1) and ydata(ny) to actually do interpolation 1366 1366 ! y within ydata(1) and ydata(ny) to actually do interpolation
! otherwise extrapolation is done 1367 1367 ! otherwise extrapolation is done
integer(kind=4), intent(in) :: nx,ny 1368 1368 integer(kind=4), intent(in) :: nx,ny
real(kind=8), intent(in) :: x,y 1369 1369 real(kind=8), intent(in) :: x,y
real(kind=8), intent(in), dimension(nx) :: xdata 1370 1370 real(kind=8), intent(in), dimension(nx) :: xdata
real(kind=8), intent(in), dimension(ny) :: ydata 1371 1371 real(kind=8), intent(in), dimension(ny) :: ydata
real(kind=8), intent(in), dimension(nx,ny) :: zdata 1372 1372 real(kind=8), intent(in), dimension(nx,ny) :: zdata
real(kind=8) :: fvn_d_quad_2d_interpol 1373 1373 real(kind=8) :: fvn_d_quad_2d_interpol
1374 1374
integer(kind=4) :: ixinf,iyinf,basex,basey,i 1375 1375 integer(kind=4) :: ixinf,iyinf,basex,basey,i
real(kind=8),dimension(3) :: ztmp 1376 1376 real(kind=8),dimension(3) :: ztmp
!real(kind=8), external :: fvn_d_quad_interpol 1377 1377 !real(kind=8), external :: fvn_d_quad_interpol
1378 1378
call fvn_d_find_interval(x,ixinf,xdata,nx) 1379 1379 call fvn_d_find_interval(x,ixinf,xdata,nx)
call fvn_d_find_interval(y,iyinf,ydata,ny) 1380 1380 call fvn_d_find_interval(y,iyinf,ydata,ny)
1381 1381
! Settings for extrapolation 1382 1382 ! Settings for extrapolation
if (ixinf==0) then 1383 1383 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 1384 1384 ! TODO -> Lower x bound extrapolation warning
ixinf=1 1385 1385 ixinf=1
end if 1386 1386 end if
1387 1387
if (ixinf==nx) then 1388 1388 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 1389 1389 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 1390 1390 ixinf=nx-1
end if 1391 1391 end if
1392 1392
if (iyinf==0) then 1393 1393 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 1394 1394 ! TODO -> Lower y bound extrapolation warning
iyinf=1 1395 1395 iyinf=1
end if 1396 1396 end if
1397 1397
if (iyinf==ny) then 1398 1398 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 1399 1399 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 1400 1400 iyinf=ny-1
end if 1401 1401 end if
1402 1402
! The three points we will use are iinf-1,iinf and iinf+1 with the 1403 1403 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 1404 1404 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 1405 1405 if (ixinf==1) then
basex=0 1406 1406 basex=0
else 1407 1407 else
basex=ixinf-2 1408 1408 basex=ixinf-2
end if 1409 1409 end if
1410 1410
if (iyinf==1) then 1411 1411 if (iyinf==1) then
basey=0 1412 1412 basey=0
else 1413 1413 else
basey=iyinf-2 1414 1414 basey=iyinf-2
end if 1415 1415 end if
1416 1416
! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3) 1417 1417 ! First we make 3 interpolations for x at y(base+1),y(base+2),y(base+3)
! stored in ztmp(1:3) 1418 1418 ! stored in ztmp(1:3)
do i=1,3 1419 1419 do i=1,3
ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i)) 1420 1420 ztmp(i)=fvn_d_quad_interpol(x,nx,xdata,zdata(:,basey+i))
end do 1421 1421 end do
1422 1422
! Then we make an interpolation for y using previous interpolations 1423 1423 ! Then we make an interpolation for y using previous interpolations
fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp) 1424 1424 fvn_d_quad_2d_interpol=fvn_d_quad_interpol(y,3,ydata(basey+1:basey+3),ztmp)
end function 1425 1425 end function
1426 1426
1427 1427
function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata) 1428 1428 function fvn_s_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
implicit none 1429 1429 implicit none
! This function evaluate the value of a 3 variables function defined by a 1430 1430 ! This function evaluate the value of a 3 variables function defined by a
! set of points and values, using a quadratic interpolation 1431 1431 ! set of points and values, using a quadratic interpolation
! xdata, ydata and zdata must be increasingly ordered 1432 1432 ! xdata, ydata and zdata must be increasingly ordered
! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 1433 1433 ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually
! perform an interpolation, otherwise extrapolation is done 1434 1434 ! perform an interpolation, otherwise extrapolation is done
integer(kind=4), intent(in) :: nx,ny,nz 1435 1435 integer(kind=4), intent(in) :: nx,ny,nz
real(kind=4), intent(in) :: x,y,z 1436 1436 real(kind=4), intent(in) :: x,y,z
real(kind=4), intent(in), dimension(nx) :: xdata 1437 1437 real(kind=4), intent(in), dimension(nx) :: xdata
real(kind=4), intent(in), dimension(ny) :: ydata 1438 1438 real(kind=4), intent(in), dimension(ny) :: ydata
real(kind=4), intent(in), dimension(nz) :: zdata 1439 1439 real(kind=4), intent(in), dimension(nz) :: zdata
real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata 1440 1440 real(kind=4), intent(in), dimension(nx,ny,nz) :: tdata
real(kind=4) :: fvn_s_quad_3d_interpol 1441 1441 real(kind=4) :: fvn_s_quad_3d_interpol
1442 1442
integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j 1443 1443 integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
!real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol 1444 1444 !real(kind=4), external :: fvn_s_quad_interpol,fvn_s_quad_2d_interpol
real(kind=4),dimension(3,3) :: ttmp 1445 1445 real(kind=4),dimension(3,3) :: ttmp
1446 1446
call fvn_s_find_interval(x,ixinf,xdata,nx) 1447 1447 call fvn_s_find_interval(x,ixinf,xdata,nx)
call fvn_s_find_interval(y,iyinf,ydata,ny) 1448 1448 call fvn_s_find_interval(y,iyinf,ydata,ny)
call fvn_s_find_interval(z,izinf,zdata,nz) 1449 1449 call fvn_s_find_interval(z,izinf,zdata,nz)
1450 1450
! Settings for extrapolation 1451 1451 ! Settings for extrapolation
if (ixinf==0) then 1452 1452 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 1453 1453 ! TODO -> Lower x bound extrapolation warning
ixinf=1 1454 1454 ixinf=1
end if 1455 1455 end if
1456 1456
if (ixinf==nx) then 1457 1457 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 1458 1458 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 1459 1459 ixinf=nx-1
end if 1460 1460 end if
1461 1461
if (iyinf==0) then 1462 1462 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 1463 1463 ! TODO -> Lower y bound extrapolation warning
iyinf=1 1464 1464 iyinf=1
end if 1465 1465 end if
1466 1466
if (iyinf==ny) then 1467 1467 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 1468 1468 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 1469 1469 iyinf=ny-1
end if 1470 1470 end if
1471 1471
if (izinf==0) then 1472 1472 if (izinf==0) then
! TODO -> Lower z bound extrapolation warning 1473 1473 ! TODO -> Lower z bound extrapolation warning
izinf=1 1474 1474 izinf=1
end if 1475 1475 end if
1476 1476
if (izinf==nz) then 1477 1477 if (izinf==nz) then
! TODO -> Higher z bound extrapolation warning 1478 1478 ! TODO -> Higher z bound extrapolation warning
izinf=nz-1 1479 1479 izinf=nz-1
end if 1480 1480 end if
1481 1481
! The three points we will use are iinf-1,iinf and iinf+1 with the 1482 1482 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 1483 1483 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 1484 1484 if (ixinf==1) then
basex=0 1485 1485 basex=0
else 1486 1486 else
basex=ixinf-2 1487 1487 basex=ixinf-2
end if 1488 1488 end if
1489 1489
if (iyinf==1) then 1490 1490 if (iyinf==1) then
basey=0 1491 1491 basey=0
else 1492 1492 else
basey=iyinf-2 1493 1493 basey=iyinf-2
end if 1494 1494 end if
1495 1495
if (izinf==1) then 1496 1496 if (izinf==1) then
basez=0 1497 1497 basez=0
else 1498 1498 else
basez=izinf-2 1499 1499 basez=izinf-2
end if 1500 1500 end if
1501 1501
! We first make 9 one dimensional interpolation on variable x. 1502 1502 ! We first make 9 one dimensional interpolation on variable x.
! results are stored in ttmp 1503 1503 ! results are stored in ttmp
do i=1,3 1504 1504 do i=1,3
do j=1,3 1505 1505 do j=1,3
ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j)) 1506 1506 ttmp(i,j)=fvn_s_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
end do 1507 1507 end do
end do 1508 1508 end do
1509 1509
! We then make a 2 dimensionnal interpolation on variables y and z 1510 1510 ! We then make a 2 dimensionnal interpolation on variables y and z
fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, & 1511 1511 fvn_s_quad_3d_interpol=fvn_s_quad_2d_interpol(y,z, &
3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp) 1512 1512 3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
end function 1513 1513 end function
1514 1514
1515 1515
function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata) 1516 1516 function fvn_d_quad_3d_interpol(x,y,z,nx,xdata,ny,ydata,nz,zdata,tdata)
implicit none 1517 1517 implicit none
! This function evaluate the value of a 3 variables function defined by a 1518 1518 ! This function evaluate the value of a 3 variables function defined by a
! set of points and values, using a quadratic interpolation 1519 1519 ! set of points and values, using a quadratic interpolation
! xdata, ydata and zdata must be increasingly ordered 1520 1520 ! xdata, ydata and zdata must be increasingly ordered
! The triplet (x,y,z) must be within xdata,ydata and zdata to actually 1521 1521 ! The triplet (x,y,z) must be within xdata,ydata and zdata to actually
! perform an interpolation, otherwise extrapolation is done 1522 1522 ! perform an interpolation, otherwise extrapolation is done
integer(kind=4), intent(in) :: nx,ny,nz 1523 1523 integer(kind=4), intent(in) :: nx,ny,nz
real(kind=8), intent(in) :: x,y,z 1524 1524 real(kind=8), intent(in) :: x,y,z
real(kind=8), intent(in), dimension(nx) :: xdata 1525 1525 real(kind=8), intent(in), dimension(nx) :: xdata
real(kind=8), intent(in), dimension(ny) :: ydata 1526 1526 real(kind=8), intent(in), dimension(ny) :: ydata
real(kind=8), intent(in), dimension(nz) :: zdata 1527 1527 real(kind=8), intent(in), dimension(nz) :: zdata
real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata 1528 1528 real(kind=8), intent(in), dimension(nx,ny,nz) :: tdata
real(kind=8) :: fvn_d_quad_3d_interpol 1529 1529 real(kind=8) :: fvn_d_quad_3d_interpol
1530 1530
integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j 1531 1531 integer(kind=4) :: ixinf,iyinf,izinf,basex,basey,basez,i,j
!real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol 1532 1532 !real(kind=8), external :: fvn_d_quad_interpol,fvn_d_quad_2d_interpol
real(kind=8),dimension(3,3) :: ttmp 1533 1533 real(kind=8),dimension(3,3) :: ttmp
1534 1534
call fvn_d_find_interval(x,ixinf,xdata,nx) 1535 1535 call fvn_d_find_interval(x,ixinf,xdata,nx)
call fvn_d_find_interval(y,iyinf,ydata,ny) 1536 1536 call fvn_d_find_interval(y,iyinf,ydata,ny)
call fvn_d_find_interval(z,izinf,zdata,nz) 1537 1537 call fvn_d_find_interval(z,izinf,zdata,nz)
1538 1538
! Settings for extrapolation 1539 1539 ! Settings for extrapolation
if (ixinf==0) then 1540 1540 if (ixinf==0) then
! TODO -> Lower x bound extrapolation warning 1541 1541 ! TODO -> Lower x bound extrapolation warning
ixinf=1 1542 1542 ixinf=1
end if 1543 1543 end if
1544 1544
if (ixinf==nx) then 1545 1545 if (ixinf==nx) then
! TODO -> Higher x bound extrapolation warning 1546 1546 ! TODO -> Higher x bound extrapolation warning
ixinf=nx-1 1547 1547 ixinf=nx-1
end if 1548 1548 end if
1549 1549
if (iyinf==0) then 1550 1550 if (iyinf==0) then
! TODO -> Lower y bound extrapolation warning 1551 1551 ! TODO -> Lower y bound extrapolation warning
iyinf=1 1552 1552 iyinf=1
end if 1553 1553 end if
1554 1554
if (iyinf==ny) then 1555 1555 if (iyinf==ny) then
! TODO -> Higher y bound extrapolation warning 1556 1556 ! TODO -> Higher y bound extrapolation warning
iyinf=ny-1 1557 1557 iyinf=ny-1
end if 1558 1558 end if
1559 1559
if (izinf==0) then 1560 1560 if (izinf==0) then
! TODO -> Lower z bound extrapolation warning 1561 1561 ! TODO -> Lower z bound extrapolation warning
izinf=1 1562 1562 izinf=1
end if 1563 1563 end if
1564 1564
if (izinf==nz) then 1565 1565 if (izinf==nz) then
! TODO -> Higher z bound extrapolation warning 1566 1566 ! TODO -> Higher z bound extrapolation warning
izinf=nz-1 1567 1567 izinf=nz-1
end if 1568 1568 end if
1569 1569
! The three points we will use are iinf-1,iinf and iinf+1 with the 1570 1570 ! The three points we will use are iinf-1,iinf and iinf+1 with the
! exception of the first interval, where iinf=1 we will use 1,2 and 3 1571 1571 ! exception of the first interval, where iinf=1 we will use 1,2 and 3
if (ixinf==1) then 1572 1572 if (ixinf==1) then
basex=0 1573 1573 basex=0
else 1574 1574 else
basex=ixinf-2 1575 1575 basex=ixinf-2
end if 1576 1576 end if
1577 1577
if (iyinf==1) then 1578 1578 if (iyinf==1) then
basey=0 1579 1579 basey=0
else 1580 1580 else
basey=iyinf-2 1581 1581 basey=iyinf-2
end if 1582 1582 end if
1583 1583
if (izinf==1) then 1584 1584 if (izinf==1) then
basez=0 1585 1585 basez=0
else 1586 1586 else
basez=izinf-2 1587 1587 basez=izinf-2
end if 1588 1588 end if
1589 1589
! We first make 9 one dimensional interpolation on variable x. 1590 1590 ! We first make 9 one dimensional interpolation on variable x.
! results are stored in ttmp 1591 1591 ! results are stored in ttmp
do i=1,3 1592 1592 do i=1,3
do j=1,3 1593 1593 do j=1,3
ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j)) 1594 1594 ttmp(i,j)=fvn_d_quad_interpol(x,nx,xdata,tdata(:,basey+i,basez+j))
end do 1595 1595 end do
end do 1596 1596 end do
1597 1597
! We then make a 2 dimensionnal interpolation on variables y and z 1598 1598 ! We then make a 2 dimensionnal interpolation on variables y and z
fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, & 1599 1599 fvn_d_quad_3d_interpol=fvn_d_quad_2d_interpol(y,z, &
3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp) 1600 1600 3,ydata(basey+1:basey+3),3,zdata(basez+1:basez+3),ttmp)
end function 1601 1601 end function
1602 1602
1603 1603
1604 1604
1605 1605
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1606 1606 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1607 1607 !
! Akima spline interpolation and spline evaluation 1608 1608 ! Akima spline interpolation and spline evaluation
! 1609 1609 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1610 1610 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1611 1611
! Single precision 1612 1612 ! Single precision
subroutine fvn_s_akima(n,x,y,br,co) 1613 1613 subroutine fvn_s_akima(n,x,y,br,co)
implicit none 1614 1614 implicit none
integer, intent(in) :: n 1615 1615 integer, intent(in) :: n
real, intent(in) :: x(n) 1616 1616 real, intent(in) :: x(n)
real, intent(in) :: y(n) 1617 1617 real, intent(in) :: y(n)
real, intent(out) :: br(n) 1618 1618 real, intent(out) :: br(n)
real, intent(out) :: co(4,n) 1619 1619 real, intent(out) :: co(4,n)
1620 1620
real, allocatable :: var(:),z(:) 1621 1621 real, allocatable :: var(:),z(:)
real :: wi_1,wi 1622 1622 real :: wi_1,wi
integer :: i 1623 1623 integer :: i
real :: dx,a,b 1624 1624 real :: dx,a,b
1625 1625
! br is just a copy of x 1626 1626 ! br is just a copy of x
br(:)=x(:) 1627 1627 br(:)=x(:)
1628 1628
allocate(var(n)) 1629 1629 allocate(var(n+3))
allocate(z(n)) 1630 1630 allocate(z(n))
! evaluate the variations 1631 1631 ! evaluate the variations
do i=1, n-1 1632 1632 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 1633 1633 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 1634 1634 end do
var(n+2)=2.e0*var(n+1)-var(n) 1635 1635 var(n+2)=2.e0*var(n+1)-var(n)
var(n+3)=2.e0*var(n+2)-var(n+1) 1636 1636 var(n+3)=2.e0*var(n+2)-var(n+1)
var(2)=2.e0*var(3)-var(4) 1637 1637 var(2)=2.e0*var(3)-var(4)
var(1)=2.e0*var(2)-var(3) 1638 1638 var(1)=2.e0*var(2)-var(3)
1639 1639
do i = 1, n 1640 1640 do i = 1, n
wi_1=abs(var(i+3)-var(i+2)) 1641 1641 wi_1=abs(var(i+3)-var(i+2))
wi=abs(var(i+1)-var(i)) 1642 1642 wi=abs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.e0) then 1643 1643 if ((wi_1+wi).eq.0.e0) then
z(i)=(var(i+2)+var(i+1))/2.e0 1644 1644 z(i)=(var(i+2)+var(i+1))/2.e0
else 1645 1645 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 1646 1646 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 1647 1647 end if
end do 1648 1648 end do
1649 1649
do i=1, n-1 1650 1650 do i=1, n-1
dx=x(i+1)-x(i) 1651 1651 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 1652 1652 a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd 1653 1653 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 1654 1654 co(1,i)=y(i)
co(2,i)=z(i) 1655 1655 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd 1656 1656 !co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd
!co(4,i)=(a-2.*b)/dx**3 ! mรฉthode wd 1657 1657 !co(4,i)=(a-2.*b)/dx**3 ! mรฉthode wd
co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx ! mรฉthode JP Moreau 1658 1658 co(3,i)=(3.e0*var(i+2)-2.e0*z(i)-z(i+1))/dx ! mรฉthode JP Moreau
co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2 ! 1659 1659 co(4,i)=(z(i)+z(i+1)-2.e0*var(i+2))/dx**2 !
! les coefficients donnรฉs par imsl sont co(3,i)*2 et co(4,i)*6 1660 1660 ! les coefficients donnรฉs par imsl sont co(3,i)*2 et co(4,i)*6
! etrangement la fonction csval corrige et donne la bonne valeur ... 1661 1661 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 1662 1662 end do
co(1,n)=y(n) 1663 1663 co(1,n)=y(n)
co(2,n)=z(n) 1664 1664 co(2,n)=z(n)
co(3,n)=0.e0 1665 1665 co(3,n)=0.e0
co(4,n)=0.e0 1666 1666 co(4,n)=0.e0
1667 1667
deallocate(z) 1668 1668 deallocate(z)
deallocate(var) 1669 1669 deallocate(var)
1670 1670
end subroutine 1671 1671 end subroutine
1672 1672
! Double precision 1673 1673 ! Double precision
subroutine fvn_d_akima(n,x,y,br,co) 1674 1674 subroutine fvn_d_akima(n,x,y,br,co)
1675 1675
implicit none 1676 1676 implicit none
integer, intent(in) :: n 1677 1677 integer, intent(in) :: n
double precision, intent(in) :: x(n) 1678 1678 double precision, intent(in) :: x(n)
double precision, intent(in) :: y(n) 1679 1679 double precision, intent(in) :: y(n)
double precision, intent(out) :: br(n) 1680 1680 double precision, intent(out) :: br(n)
double precision, intent(out) :: co(4,n) 1681 1681 double precision, intent(out) :: co(4,n)
1682 1682
double precision, allocatable :: var(:),z(:) 1683 1683 double precision, allocatable :: var(:),z(:)
double precision :: wi_1,wi 1684 1684 double precision :: wi_1,wi
integer :: i 1685 1685 integer :: i
double precision :: dx,a,b 1686 1686 double precision :: dx,a,b
1687 1687
! br is just a copy of x 1688 1688 ! br is just a copy of x
br(:)=x(:) 1689 1689 br(:)=x(:)
1690 1690
allocate(var(n)) 1691 1691 allocate(var(n+3))
allocate(z(n)) 1692 1692 allocate(z(n))
! evaluate the variations 1693 1693 ! evaluate the variations
do i=1, n-1 1694 1694 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 1695 1695 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 1696 1696 end do
var(n+2)=2.d0*var(n+1)-var(n) 1697 1697 var(n+2)=2.d0*var(n+1)-var(n)
var(n+3)=2.d0*var(n+2)-var(n+1) 1698 1698 var(n+3)=2.d0*var(n+2)-var(n+1)
var(2)=2.d0*var(3)-var(4) 1699 1699 var(2)=2.d0*var(3)-var(4)
var(1)=2.d0*var(2)-var(3) 1700 1700 var(1)=2.d0*var(2)-var(3)
1701 1701
do i = 1, n 1702 1702 do i = 1, n
wi_1=dabs(var(i+3)-var(i+2)) 1703 1703 wi_1=dabs(var(i+3)-var(i+2))
wi=dabs(var(i+1)-var(i)) 1704 1704 wi=dabs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.d0) then 1705 1705 if ((wi_1+wi).eq.0.d0) then
z(i)=(var(i+2)+var(i+1))/2.d0 1706 1706 z(i)=(var(i+2)+var(i+1))/2.d0
else 1707 1707 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 1708 1708 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 1709 1709 end if
end do 1710 1710 end do
1711 1711
do i=1, n-1 1712 1712 do i=1, n-1
dx=x(i+1)-x(i) 1713 1713 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 1714 1714 a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd
b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd 1715 1715 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 1716 1716 co(1,i)=y(i)
co(2,i)=z(i) 1717 1717 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd 1718 1718 !co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd
!co(4,i)=(a-2.*b)/dx**3 ! mรฉthode wd 1719 1719 !co(4,i)=(a-2.*b)/dx**3 ! mรฉthode wd
co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx ! mรฉthode JP Moreau 1720 1720 co(3,i)=(3.d0*var(i+2)-2.d0*z(i)-z(i+1))/dx ! mรฉthode JP Moreau
co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2 ! 1721 1721 co(4,i)=(z(i)+z(i+1)-2.d0*var(i+2))/dx**2 !
! les coefficients donnรฉs par imsl sont co(3,i)*2 et co(4,i)*6 1722 1722 ! les coefficients donnรฉs par imsl sont co(3,i)*2 et co(4,i)*6
! etrangement la fonction csval corrige et donne la bonne valeur ... 1723 1723 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 1724 1724 end do
co(1,n)=y(n) 1725 1725 co(1,n)=y(n)
co(2,n)=z(n) 1726 1726 co(2,n)=z(n)
co(3,n)=0.d0 1727 1727 co(3,n)=0.d0
co(4,n)=0.d0 1728 1728 co(4,n)=0.d0
1729 1729
deallocate(z) 1730 1730 deallocate(z)
deallocate(var) 1731 1731 deallocate(var)
1732 1732
end subroutine 1733 1733 end subroutine
1734 1734
! 1735 1735 !
! Single precision spline evaluation 1736 1736 ! Single precision spline evaluation
! 1737 1737 !
function fvn_s_spline_eval(x,n,br,co) 1738 1738 function fvn_s_spline_eval(x,n,br,co)
implicit none 1739 1739 implicit none
real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated 1740 1740 real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
integer, intent(in) :: n ! number of intervals 1741 1741 integer, intent(in) :: n ! number of intervals
real, intent(in) :: br(n+1) ! breakpoints 1742 1742 real, intent(in) :: br(n+1) ! breakpoints
real, intent(in) :: co(4,n+1) ! spline coeeficients 1743 1743 real, intent(in) :: co(4,n+1) ! spline coeeficients
real :: fvn_s_spline_eval 1744 1744 real :: fvn_s_spline_eval
1745 1745
integer :: i 1746 1746 integer :: i
real :: dx 1747 1747 real :: dx
1748 1748
if (x<=br(1)) then 1749 1749 if (x<=br(1)) then
i=1 1750 1750 i=1
else if (x>=br(n+1)) then 1751 1751 else if (x>=br(n+1)) then
i=n 1752 1752 i=n
else 1753 1753 else
i=1 1754 1754 i=1
do while(x>=br(i)) 1755 1755 do while(x>=br(i))
i=i+1 1756 1756 i=i+1
end do 1757 1757 end do
i=i-1 1758 1758 i=i-1
end if 1759 1759 end if
dx=x-br(i) 1760 1760 dx=x-br(i)
fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3 1761 1761 fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
1762 1762
end function 1763 1763 end function
1764 1764
! Double precision spline evaluation 1765 1765 ! Double precision spline evaluation
function fvn_d_spline_eval(x,n,br,co) 1766 1766 function fvn_d_spline_eval(x,n,br,co)
implicit none 1767 1767 implicit none
double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated 1768 1768 double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
integer, intent(in) :: n ! number of intervals 1769 1769 integer, intent(in) :: n ! number of intervals
double precision, intent(in) :: br(n+1) ! breakpoints 1770 1770 double precision, intent(in) :: br(n+1) ! breakpoints
double precision, intent(in) :: co(4,n+1) ! spline coeeficients 1771 1771 double precision, intent(in) :: co(4,n+1) ! spline coeeficients
double precision :: fvn_d_spline_eval 1772 1772 double precision :: fvn_d_spline_eval
1773 1773
integer :: i 1774 1774 integer :: i
double precision :: dx 1775 1775 double precision :: dx
1776 1776
1777 1777
if (x<=br(1)) then 1778 1778 if (x<=br(1)) then
i=1 1779 1779 i=1
else if (x>=br(n+1)) then 1780 1780 else if (x>=br(n+1)) then
i=n 1781 1781 i=n
else 1782 1782 else
i=1 1783 1783 i=1
do while(x>=br(i)) 1784 1784 do while(x>=br(i))
i=i+1 1785 1785 i=i+1
end do 1786 1786 end do
i=i-1 1787 1787 i=i-1
end if 1788 1788 end if
1789 1789
dx=x-br(i) 1790 1790 dx=x-br(i)
fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3 1791 1791 fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
1792 1792
end function 1793 1793 end function
1794 1794
1795 1795
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1796 1796 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1797 1797 !
! Least square problem 1798 1798 ! Least square problem
! 1799 1799 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1800 1800 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1801 1801 !
! 1802 1802 !
1803 1803
1804 1804
1805 1805
1806 1806
subroutine fvn_d_lspoly(np,x,y,deg,coeff,status) 1807 1807 subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
! 1808 1808 !
! Least square polynomial fitting 1809 1809 ! Least square polynomial fitting
! 1810 1810 !
! Find the coefficients of the least square polynomial of a given degree 1811 1811 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1812 1812 ! for a set of coordinates.
! 1813 1813 !
! The degree must be lower than the number of points 1814 1814 ! The degree must be lower than the number of points
! 1815 1815 !
! np (in) : number of points 1816 1816 ! np (in) : number of points
! x(np) (in) : x data 1817 1817 ! x(np) (in) : x data
! y(np) (in) : y data 1818 1818 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1819 1819 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1820 1820 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1821 1821 ! status (out) : =0 if a problem occurs
! 1822 1822 !
implicit none 1823 1823 implicit none
1824 1824
integer, intent(in) :: np,deg 1825 1825 integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y 1826 1826 real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff 1827 1827 real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status 1828 1828 integer, intent(out), optional :: status
1829 1829
real(kind=8), allocatable, dimension(:,:) :: mat,bmat 1830 1830 real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work 1831 1831 real(kind=8),dimension(:),allocatable :: work
real(kind=8),dimension(1) :: twork 1832 1832 real(kind=8),dimension(1) :: twork
integer :: lwork,info 1833 1833 integer :: lwork,info
1834 1834
integer :: i,j 1835 1835 integer :: i,j
1836 1836
if (present(status)) status=1 1837 1837 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1)) 1838 1838 allocate(mat(np,deg+1),bmat(np,1))
1839 1839
! Design matrix valorisation 1840 1840 ! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) ) 1841 1841 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1842 1842
! second member valorisation 1843 1843 ! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /)) 1844 1844 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1845 1845
! query workspace size 1846 1846 ! query workspace size
call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info) 1847 1847 call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
lwork=twork(1) 1848 1848 lwork=twork(1)
allocate(work(int(lwork))) 1849 1849 allocate(work(int(lwork)))
! real call 1850 1850 ! real call
call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info) 1851 1851 call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
1852 1852
if (info /= 0) then 1853 1853 if (info /= 0) then
if (present(status)) status=0 1854 1854 if (present(status)) status=0
end if 1855 1855 end if
1856 1856
coeff = (/ (bmat(i,1),i=1,deg+1) /) 1857 1857 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1858 1858
deallocate(work) 1859 1859 deallocate(work)
deallocate(mat,bmat) 1860 1860 deallocate(mat,bmat)
end subroutine 1861 1861 end subroutine
1862 1862
subroutine fvn_s_lspoly(np,x,y,deg,coeff,status) 1863 1863 subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)
! 1864 1864 !
! Least square polynomial fitting 1865 1865 ! Least square polynomial fitting
! 1866 1866 !
! Find the coefficients of the least square polynomial of a given degree 1867 1867 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1868 1868 ! for a set of coordinates.
! 1869 1869 !
! The degree must be lower than the number of points 1870 1870 ! The degree must be lower than the number of points
! 1871 1871 !
! np (in) : number of points 1872 1872 ! np (in) : number of points
! x(np) (in) : x data 1873 1873 ! x(np) (in) : x data
! y(np) (in) : y data 1874 1874 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1875 1875 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1876 1876 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1877 1877 ! status (out) : =0 if a problem occurs
! 1878 1878 !
implicit none 1879 1879 implicit none
1880 1880
integer, intent(in) :: np,deg 1881 1881 integer, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y 1882 1882 real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff 1883 1883 real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status 1884 1884 integer, intent(out), optional :: status
1885 1885
real(kind=4), allocatable, dimension(:,:) :: mat,bmat 1886 1886 real(kind=4), allocatable, dimension(:,:) :: mat,bmat
real(kind=4),dimension(:),allocatable :: work 1887 1887 real(kind=4),dimension(:),allocatable :: work
real(kind=4),dimension(1) :: twork 1888 1888 real(kind=4),dimension(1) :: twork
integer :: lwork,info 1889 1889 integer :: lwork,info
1890 1890
integer :: i,j 1891 1891 integer :: i,j
1892 1892
if (present(status)) status=1 1893 1893 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1)) 1894 1894 allocate(mat(np,deg+1),bmat(np,1))
1895 1895
! Design matrix valorisation 1896 1896 ! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) ) 1897 1897 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1898 1898
! second member valorisation 1899 1899 ! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /)) 1900 1900 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1901 1901
! query workspace size 1902 1902 ! query workspace size
call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info) 1903 1903 call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
lwork=twork(1) 1904 1904 lwork=twork(1)
allocate(work(int(lwork))) 1905 1905 allocate(work(int(lwork)))
! real call 1906 1906 ! real call
call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info) 1907 1907 call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
1908 1908
if (info /= 0) then 1909 1909 if (info /= 0) then
if (present(status)) status=0 1910 1910 if (present(status)) status=0
end if 1911 1911 end if
1912 1912
coeff = (/ (bmat(i,1),i=1,deg+1) /) 1913 1913 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1914 1914
deallocate(work) 1915 1915 deallocate(work)
deallocate(mat,bmat) 1916 1916 deallocate(mat,bmat)
end subroutine 1917 1917 end subroutine
1918 1918
1919 1919
1920 1920
1921 1921
1922 1922
1923 1923
1924 1924
1925 1925
subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status) 1926 1926 subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)
! 1927 1927 !
! Least square polynomial fitting using singular value decomposition 1928 1928 ! Least square polynomial fitting using singular value decomposition
! 1929 1929 !
! Find the coefficients of the least square polynomial of a given degree 1930 1930 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1931 1931 ! for a set of coordinates.
! 1932 1932 !
! The degree must be lower than the number of points 1933 1933 ! The degree must be lower than the number of points
! 1934 1934 !
! np (in) : number of points 1935 1935 ! np (in) : number of points
! x(np) (in) : x data 1936 1936 ! x(np) (in) : x data
! y(np) (in) : y data 1937 1937 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1938 1938 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1939 1939 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1940 1940 ! status (out) : =0 if a problem occurs
! 1941 1941 !
implicit none 1942 1942 implicit none
1943 1943
integer, intent(in) :: np,deg 1944 1944 integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y 1945 1945 real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff 1946 1946 real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status 1947 1947 integer, intent(out), optional :: status
1948 1948
real(kind=8), allocatable, dimension(:,:) :: mat,bmat 1949 1949 real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work,singval 1950 1950 real(kind=8),dimension(:),allocatable :: work,singval
real(kind=8),dimension(1) :: twork 1951 1951 real(kind=8),dimension(1) :: twork
integer :: lwork,info,rank 1952 1952 integer :: lwork,info,rank
1953 1953
integer :: i,j 1954 1954 integer :: i,j
1955 1955
if (present(status)) status=1 1956 1956 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1)) 1957 1957 allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
1958 1958
! Design matrix valorisation 1959 1959 ! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) ) 1960 1960 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1961 1961
! second member valorisation 1962 1962 ! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /)) 1963 1963 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1964 1964
! query workspace size 1965 1965 ! query workspace size
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info) 1966 1966 call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1) 1967 1967 lwork=twork(1)
allocate(work(int(lwork))) 1968 1968 allocate(work(int(lwork)))
! real call 1969 1969 ! real call
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info) 1970 1970 call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
1971 1971
if (info /= 0) then 1972 1972 if (info /= 0) then
if (present(status)) status=0 1973 1973 if (present(status)) status=0
end if 1974 1974 end if
1975 1975
coeff = (/ (bmat(i,1),i=1,deg+1) /) 1976 1976 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1977 1977
deallocate(work) 1978 1978 deallocate(work)
deallocate(mat,bmat,singval) 1979 1979 deallocate(mat,bmat,singval)
end subroutine 1980 1980 end subroutine
1981 1981
subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status) 1982 1982 subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)
! 1983 1983 !
! Least square polynomial fitting using singular value decomposition 1984 1984 ! Least square polynomial fitting using singular value decomposition
! 1985 1985 !
! Find the coefficients of the least square polynomial of a given degree 1986 1986 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1987 1987 ! for a set of coordinates.
! 1988 1988 !
! The degree must be lower than the number of points 1989 1989 ! The degree must be lower than the number of points
! 1990 1990 !
! np (in) : number of points 1991 1991 ! np (in) : number of points
! x(np) (in) : x data 1992 1992 ! x(np) (in) : x data
! y(np) (in) : y data 1993 1993 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1994 1994 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1995 1995 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1996 1996 ! status (out) : =0 if a problem occurs
! 1997 1997 !
implicit none 1998 1998 implicit none
1999 1999
integer, intent(in) :: np,deg 2000 2000 integer, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y 2001 2001 real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff 2002 2002 real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out), optional :: status 2003 2003 integer, intent(out), optional :: status
2004 2004
real(kind=4), allocatable, dimension(:,:) :: mat,bmat 2005 2005 real(kind=4), allocatable, dimension(:,:) :: mat,bmat
real(kind=4),dimension(:),allocatable :: work,singval 2006 2006 real(kind=4),dimension(:),allocatable :: work,singval
real(kind=4),dimension(1) :: twork 2007 2007 real(kind=4),dimension(1) :: twork
integer :: lwork,info,rank 2008 2008 integer :: lwork,info,rank
2009 2009
integer :: i,j 2010 2010 integer :: i,j
2011 2011
if (present(status)) status=1 2012 2012 if (present(status)) status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1)) 2013 2013 allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
2014 2014
! Design matrix valorisation 2015 2015 ! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) ) 2016 2016 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
2017 2017
! second member valorisation 2018 2018 ! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /)) 2019 2019 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
2020 2020
! query workspace size 2021 2021 ! query workspace size
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info) 2022 2022 call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1) 2023 2023 lwork=twork(1)
allocate(work(int(lwork))) 2024 2024 allocate(work(int(lwork)))
! real call 2025 2025 ! real call
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info) 2026 2026 call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
2027 2027
if (info /= 0) then 2028 2028 if (info /= 0) then
if (present(status)) status=0 2029 2029 if (present(status)) status=0
end if 2030 2030 end if
2031 2031
coeff = (/ (bmat(i,1),i=1,deg+1) /) 2032 2032 coeff = (/ (bmat(i,1),i=1,deg+1) /)
2033 2033
deallocate(work) 2034 2034 deallocate(work)
deallocate(mat,bmat,singval) 2035 2035 deallocate(mat,bmat,singval)
end subroutine 2036 2036 end subroutine
2037 2037
2038 2038
! 2039 2039 !
! Muller 2040 2040 ! Muller
! 2041 2041 !
! 2042 2042 !
! 2043 2043 !
! William Daniau 2007 2044 2044 ! William Daniau 2007
! 2045 2045 !
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f 2046 2046 ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f 2047 2047 ! http://plato.asu.edu/ftp/other_software/muller.f
! 2048 2048 !
! it can be used as a replacement for imsl routine dzanly with minor changes 2049 2049 ! it can be used as a replacement for imsl routine dzanly with minor changes
! 2050 2050 !
!----------------------------------------------------------------------- 2051 2051 !-----------------------------------------------------------------------
! 2052 2052 !
! purpose - zeros of an analytic complex function 2053 2053 ! purpose - zeros of an analytic complex function
! using the muller method with deflation 2054 2054 ! using the muller method with deflation
! 2055 2055 !
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, 2056 2056 ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier) 2057 2057 ! infer,ier)
! 2058 2058 !
! arguments f - a complex function subprogram, f(z), written 2059 2059 ! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose 2060 2060 ! by the user specifying the equation whose
! roots are to be found. f must appear in 2061 2061 ! roots are to be found. f must appear in
! an external statement in the calling pro- 2062 2062 ! an external statement in the calling pro-
! gram. 2063 2063 ! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p 2064 2064 ! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) 2065 2065 ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found 2066 2066 ! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and. 2067 2067 ! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted 2068 2068 ! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input) 2069 2069 ! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted 2070 2070 ! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given 2071 2071 ! if two successive approximations to a given
! root agree within eps1. (input) 2072 2072 ! root agree within eps1. (input)
! note. if either or both of the stopping 2073 2073 ! note. if either or both of the stopping
! criteria are fulfilled, the root is 2074 2074 ! criteria are fulfilled, the root is
! accepted. 2075 2075 ! accepted.
! kn - the number of known roots which must be stored 2076 2076 ! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller 2077 2077 ! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these 2078 2078 ! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),..., 2079 2079 ! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal 2080 2080 ! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input) 2081 2081 ! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by 2082 2082 ! n - the number of new roots to be found by
! muller (input) 2083 2083 ! muller (input)
! x - a complex vector of length kn+n. x(1),..., 2084 2084 ! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known 2085 2085 ! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may, 2086 2086 ! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for 2087 2087 ! on user option, contain initial guesses for
! the n new roots which are to be computed. 2088 2088 ! the n new roots which are to be computed.
! if the user does not provide an initial 2089 2089 ! if the user does not provide an initial
! guess, zero is used. 2090 2090 ! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the 2091 2091 ! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller. 2092 2092 ! approximate roots found by muller.
! itmax - the maximum allowable number of iterations 2093 2093 ! itmax - the maximum allowable number of iterations
! per root (input) 2094 2094 ! per root (input)
! infer - an integer vector of length kn+n. on 2095 2095 ! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of 2096 2096 ! output infer(j) contains the number of
! iterations used in finding the j-th root 2097 2097 ! iterations used in finding the j-th root
! when convergence was achieved. if 2098 2098 ! when convergence was achieved. if
! convergence was not obtained in itmax 2099 2099 ! convergence was not obtained in itmax
! iterations, infer(j) will be greater than 2100 2100 ! iterations, infer(j) will be greater than
! itmax (output). 2101 2101 ! itmax (output).
! ier - error parameter (output) 2102 2102 ! ier - error parameter (output)
! warning error 2103 2103 ! warning error
! ier = 33 indicates failure to converge with- 2104 2104 ! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of 2105 2105 ! in itmax iterations for at least one of
! the (n) new roots. 2106 2106 ! the (n) new roots.
! 2107 2107 !
! 2108 2108 !
! remarks muller always returns the last approximation for root j 2109 2109 ! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied, 2110 2110 ! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the 2111 2111 ! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j) 2112 2112 ! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater 2113 2113 ! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did 2114 2114 ! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter- 2115 2115 ! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax 2116 2116 ! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con- 2117 2117 ! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla- 2118 2118 ! vergence was obtained (on iteration k) for the defla-
! ted function 2119 2119 ! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) 2120 2120 ! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
! 2121 2121 !
! but failed for f(z). in this case, better initial 2122 2122 ! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax 2123 2123 ! guesses might help or, it might be necessary to relax
! the convergence criterion. 2124 2124 ! the convergence criterion.
! 2125 2125 !
!----------------------------------------------------------------------- 2126 2126 !-----------------------------------------------------------------------
! 2127 2127 !
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) 2128 2128 subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none 2129 2129 implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq 2130 2130 double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, & 2131 2131 double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & 2132 2132 tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5 2133 2133 zero,p1,one,four,p5
2134 2134
double complex, external :: f 2135 2135 double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & 2136 2136 integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode 2137 2137 knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n) 2138 2138 double complex :: x(kn+n)
integer :: infer(kn+n) 2139 2139 integer :: infer(kn+n)
2140 2140
2141 2141
data zero/(0.0,0.0)/,p1/(0.1,0.0)/, & 2142 2142 data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
one/(1.0,0.0)/,four/(4.0,0.0)/, & 2143 2143 one/(1.0,0.0)/,four/(4.0,0.0)/, &
p5/(0.5,0.0)/, & 2144 2144 p5/(0.5,0.0)/, &
rzero/0.0/,rten/10.0/,rhun/100.0/, & 2145 2145 rzero/0.0/,rten/10.0/,rhun/100.0/, &
ax/0.1/,ickmax/3/,rp01/0.01/ 2146 2146 ax/0.1/,ickmax/3/,rp01/0.01/
2147 2147
ier = 0 2148 2148 ier = 0
if (n .lt. 1) then ! What the hell are doing here then ... 2149 2149 if (n .lt. 1) then ! What the hell are doing here then ...
return 2150 2150 return
end if 2151 2151 end if
!eps1 = rten **(-nsig) 2152 2152 !eps1 = rten **(-nsig)
eps1 = min(eps1,rp01) 2153 2153 eps1 = min(eps1,rp01)
2154 2154
knp1 = kn+1 2155 2155 knp1 = kn+1
knpn = kn+n 2156 2156 knpn = kn+n
knpng = kn+nguess 2157 2157 knpng = kn+nguess
do i=1,knpn 2158 2158 do i=1,knpn
infer(i) = 0 2159 2159 infer(i) = 0
if (i .gt. knpng) x(i) = zero 2160 2160 if (i .gt. knpng) x(i) = zero
end do 2161 2161 end do
l= knp1 2162 2162 l= knp1
2163 2163
ic=0 2164 2164 ic=0
rloop: do while (l<=knpn) ! Main loop over new roots 2165 2165 rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0 2166 2166 jk = 0
ick = 0 2167 2167 ick = 0
xl = x(l) 2168 2168 xl = x(l)
icloop: do 2169 2169 icloop: do
ic = 0 2170 2170 ic = 0
h = ax 2171 2171 h = ax
h = p1*h 2172 2172 h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl 2173 2173 if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are 2174 2174 ! first three points are
! xl+h, xl-h, xl 2175 2175 ! xl+h, xl-h, xl
rt = xl+h 2176 2176 rt = xl+h
call deflated_work(errcode) 2177 2177 call deflated_work(errcode)
if (errcode == 1) then 2178 2178 if (errcode == 1) then
exit icloop 2179 2179 exit icloop
end if 2180 2180 end if
2181 2181
z0 = fprt 2182 2182 z0 = fprt
y0 = frt 2183 2183 y0 = frt
x0 = rt 2184 2184 x0 = rt
rt = xl-h 2185 2185 rt = xl-h
call deflated_work(errcode) 2186 2186 call deflated_work(errcode)
if (errcode == 1) then 2187 2187 if (errcode == 1) then
exit icloop 2188 2188 exit icloop
end if 2189 2189 end if
2190 2190
z1 = fprt 2191 2191 z1 = fprt
y1 = frt 2192 2192 y1 = frt
h = xl-rt 2193 2193 h = xl-rt
d = h/(rt-x0) 2194 2194 d = h/(rt-x0)
rt = xl 2195 2195 rt = xl
2196 2196
call deflated_work(errcode) 2197 2197 call deflated_work(errcode)
if (errcode == 1) then 2198 2198 if (errcode == 1) then
exit icloop 2199 2199 exit icloop
end if 2200 2200 end if
2201 2201
2202 2202
z2 = fprt 2203 2203 z2 = fprt
y2 = frt 2204 2204 y2 = frt
! begin main algorithm 2205 2205 ! begin main algorithm
iloop: do 2206 2206 iloop: do
dd = one + d 2207 2207 dd = one + d
t1 = z0*d*d 2208 2208 t1 = z0*d*d
t2 = z1*dd*dd 2209 2209 t2 = z1*dd*dd
xx = z2*dd 2210 2210 xx = z2*dd
t3 = z2*d 2211 2211 t3 = z2*d
bi = t1-t2+xx+t3 2212 2212 bi = t1-t2+xx+t3
den = bi*bi-four*(xx*t1-t3*(t2-xx)) 2213 2213 den = bi*bi-four*(xx*t1-t3*(t2-xx))
! use denominator of maximum amplitude 2214 2214 ! use denominator of maximum amplitude
t1 = cdsqrt(den) 2215 2215 t1 = cdsqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1)) 2216 2216 qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1 2217 2217 t2 = bi + t1
tpq = cdabs(t2)+qz 2218 2218 tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero 2219 2219 if (tpq .eq. qz) t2 = zero
t3 = bi - t1 2220 2220 t3 = bi - t1
tpq = cdabs(t3) + qz 2221 2221 tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero 2222 2222 if (tpq .eq. qz) t3 = zero
den = t2 2223 2223 den = t2
qz = cdabs(t3)-cdabs(t2) 2224 2224 qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3 2225 2225 if (qz .gt. rzero) den = t3
! test for zero denominator 2226 2226 ! test for zero denominator
if (cdabs(den) .eq. rzero) then 2227 2227 if (cdabs(den) .eq. rzero) then
call trans_rt() 2228 2228 call trans_rt()
call deflated_work(errcode) 2229 2229 call deflated_work(errcode)
if (errcode == 1) then 2230 2230 if (errcode == 1) then
exit icloop 2231 2231 exit icloop
end if 2232 2232 end if
z2 = fprt 2233 2233 z2 = fprt
y2 = frt 2234 2234 y2 = frt
cycle iloop 2235 2235 cycle iloop
end if 2236 2236 end if
2237 2237
2238 2238
d = -xx/den 2239 2239 d = -xx/den
d = d+d 2240 2240 d = d+d
h = d*h 2241 2241 h = d*h
rt = rt + h 2242 2242 rt = rt + h
! check convergence of the first kind 2243 2243 ! check convergence of the first kind
if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then 2244 2244 if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
if (ic .ne. 0) then 2245 2245 if (ic .ne. 0) then
exit icloop 2246 2246 exit icloop
end if 2247 2247 end if
ic = 1 2248 2248 ic = 1
z0 = y1 2249 2249 z0 = y1
z1 = y2 2250 2250 z1 = y2
z2 = f(rt) 2251 2251 z2 = f(rt)
xl = rt 2252 2252 xl = rt
ick = ick+1 2253 2253 ick = ick+1
if (ick .le. ickmax) then 2254 2254 if (ick .le. ickmax) then
cycle iloop 2255 2255 cycle iloop
end if 2256 2256 end if
! warning error, itmax = maximum 2257 2257 ! warning error, itmax = maximum
jk = itmax + jk 2258 2258 jk = itmax + jk
ier = 33 2259 2259 ier = 33
end if 2260 2260 end if
if (ic .ne. 0) then 2261 2261 if (ic .ne. 0) then
cycle icloop 2262 2262 cycle icloop
end if 2263 2263 end if
call deflated_work(errcode) 2264 2264 call deflated_work(errcode)
if (errcode == 1) then 2265 2265 if (errcode == 1) then
exit icloop 2266 2266 exit icloop
end if 2267 2267 end if
2268 2268
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) 2269 2269 do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce 2270 2270 ! take remedial action to induce
! convergence 2271 2271 ! convergence
d = d*p5 2272 2272 d = d*p5
h = h*p5 2273 2273 h = h*p5
rt = rt-h 2274 2274 rt = rt-h
call deflated_work(errcode) 2275 2275 call deflated_work(errcode)
if (errcode == 1) then 2276 2276 if (errcode == 1) then
exit icloop 2277 2277 exit icloop
end if 2278 2278 end if
end do 2279 2279 end do
z0 = z1 2280 2280 z0 = z1
z1 = z2 2281 2281 z1 = z2
z2 = fprt 2282 2282 z2 = fprt
y0 = y1 2283 2283 y0 = y1
y1 = y2 2284 2284 y1 = y2
y2 = frt 2285 2285 y2 = frt
end do iloop 2286 2286 end do iloop
end do icloop 2287 2287 end do icloop
x(l) = rt 2288 2288 x(l) = rt
infer(l) = jk 2289 2289 infer(l) = jk
l = l+1 2290 2290 l = l+1
end do rloop 2291 2291 end do rloop
2292 2292
contains 2293 2293 contains
subroutine trans_rt() 2294 2294 subroutine trans_rt()
tem = rten*eps1 2295 2295 tem = rten*eps1
if (cdabs(rt) .gt. ax) tem = tem*rt 2296 2296 if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem 2297 2297 rt = rt+tem
d = (h+tem)*d/h 2298 2298 d = (h+tem)*d/h
h = h+tem 2299 2299 h = h+tem
end subroutine trans_rt 2300 2300 end subroutine trans_rt
2301 2301
subroutine deflated_work(errcode) 2302 2302 subroutine deflated_work(errcode)
! errcode=0 => no errors 2303 2303 ! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved 2304 2304 ! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag 2305 2305 integer :: errcode,flag
2306 2306
flag=1 2307 2307 flag=1
loop1: do while(flag==1) 2308 2308 loop1: do while(flag==1)
errcode=0 2309 2309 errcode=0
jk = jk+1 2310 2310 jk = jk+1
if (jk .gt. itmax) then 2311 2311 if (jk .gt. itmax) then
ier=33 2312 2312 ier=33
errcode=1 2313 2313 errcode=1
return 2314 2314 return
end if 2315 2315 end if
frt = f(rt) 2316 2316 frt = f(rt)
fprt = frt 2317 2317 fprt = frt
if (l /= 1) then 2318 2318 if (l /= 1) then
lm1 = l-1 2319 2319 lm1 = l-1
do i=1,lm1 2320 2320 do i=1,lm1
tem = rt - x(i) 2321 2321 tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then 2322 2322 if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible? 2323 2323 !if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt() 2324 2324 call trans_rt()
cycle loop1 2325 2325 cycle loop1
end if 2326 2326 end if
fprt = fprt/tem 2327 2327 fprt = fprt/tem
end do 2328 2328 end do
end if 2329 2329 end if
flag=0 2330 2330 flag=0
end do loop1 2331 2331 end do loop1
2332 2332
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then 2333 2333 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1 2334 2334 errcode=1
return 2335 2335 return
end if 2336 2336 end if
2337 2337
end subroutine deflated_work 2338 2338 end subroutine deflated_work
2339 2339
end subroutine 2340 2340 end subroutine
2341 2341
2342 2342
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2343 2343 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 2344 2344 !
! Integration 2345 2345 ! Integration
! 2346 2346 !
! Only double precision coded atm 2347 2347 ! Only double precision coded atm
! 2348 2348 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2349 2349 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2350 2350
2351 2351
subroutine fvn_d_gauss_legendre(n,qx,qw) 2352 2352 subroutine fvn_d_gauss_legendre(n,qx,qw)
! 2353 2353 !
! This routine compute the n Gauss Legendre abscissas and weights 2354 2354 ! This routine compute the n Gauss Legendre abscissas and weights
! Adapted from Numerical Recipes routine gauleg 2355 2355 ! Adapted from Numerical Recipes routine gauleg
! 2356 2356 !
! n (in) : number of points 2357 2357 ! n (in) : number of points
! qx(out) : abscissas 2358 2358 ! qx(out) : abscissas
! qw(out) : weights 2359 2359 ! qw(out) : weights
! 2360 2360 !
implicit none 2361 2361 implicit none
double precision,parameter :: pi=3.141592653589793d0 2362 2362 double precision,parameter :: pi=3.141592653589793d0
integer, intent(in) :: n 2363 2363 integer, intent(in) :: n
double precision, intent(out) :: qx(n),qw(n) 2364 2364 double precision, intent(out) :: qx(n),qw(n)
2365 2365
integer :: m,i,j 2366 2366 integer :: m,i,j
double precision :: z,z1,p1,p2,p3,pp 2367 2367 double precision :: z,z1,p1,p2,p3,pp
m=(n+1)/2 2368 2368 m=(n+1)/2
do i=1,m 2369 2369 do i=1,m
z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0)) 2370 2370 z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
iloop: do 2371 2371 iloop: do
p1=1.d0 2372 2372 p1=1.d0
p2=0.d0 2373 2373 p2=0.d0
do j=1,n 2374 2374 do j=1,n
p3=p2 2375 2375 p3=p2
p2=p1 2376 2376 p2=p1
p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j) 2377 2377 p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
end do 2378 2378 end do
pp=dble(n)*(z*p1-p2)/(z*z-1.d0) 2379 2379 pp=dble(n)*(z*p1-p2)/(z*z-1.d0)
z1=z 2380 2380 z1=z
z=z1-p1/pp 2381 2381 z=z1-p1/pp
if (dabs(z-z1)<=epsilon(z)) then 2382 2382 if (dabs(z-z1)<=epsilon(z)) then
exit iloop 2383 2383 exit iloop
end if 2384 2384 end if
end do iloop 2385 2385 end do iloop
qx(i)=-z 2386 2386 qx(i)=-z
qx(n+1-i)=z 2387 2387 qx(n+1-i)=z
qw(i)=2.d0/((1.d0-z*z)*pp*pp) 2388 2388 qw(i)=2.d0/((1.d0-z*z)*pp*pp)
qw(n+1-i)=qw(i) 2389 2389 qw(n+1-i)=qw(i)
end do 2390 2390 end do
end subroutine 2391 2391 end subroutine
2392 2392
2393 2393
2394 2394
subroutine fvn_d_gl_integ(f,a,b,n,res) 2395 2395 subroutine fvn_d_gl_integ(f,a,b,n,res)
! 2396 2396 !
! This is a simple non adaptative integration routine 2397 2397 ! This is a simple non adaptative integration routine
! using n gauss legendre abscissas and weights 2398 2398 ! using n gauss legendre abscissas and weights
! 2399 2399 !
! f(in) : the function to integrate 2400 2400 ! f(in) : the function to integrate
! a(in) : lower bound 2401 2401 ! a(in) : lower bound
! b(in) : higher bound 2402 2402 ! b(in) : higher bound
! n(in) : number of gauss legendre pairs 2403 2403 ! n(in) : number of gauss legendre pairs
! res(out): the evaluation of the integral 2404 2404 ! res(out): the evaluation of the integral
! 2405 2405 !
double precision,external :: f 2406 2406 double precision,external :: f
double precision, intent(in) :: a,b 2407 2407 double precision, intent(in) :: a,b
integer, intent(in):: n 2408 2408 integer, intent(in):: n
double precision, intent(out) :: res 2409 2409 double precision, intent(out) :: res
2410 2410
double precision, allocatable :: qx(:),qw(:) 2411 2411 double precision, allocatable :: qx(:),qw(:)
double precision :: xm,xr 2412 2412 double precision :: xm,xr
integer :: i 2413 2413 integer :: i
2414 2414
! First compute n gauss legendre abs and weight 2415 2415 ! First compute n gauss legendre abs and weight
allocate(qx(n)) 2416 2416 allocate(qx(n))
allocate(qw(n)) 2417 2417 allocate(qw(n))
call fvn_d_gauss_legendre(n,qx,qw) 2418 2418 call fvn_d_gauss_legendre(n,qx,qw)
2419 2419
xm=0.5d0*(b+a) 2420 2420 xm=0.5d0*(b+a)
xr=0.5d0*(b-a) 2421 2421 xr=0.5d0*(b-a)
2422 2422
res=0.d0 2423 2423 res=0.d0
2424 2424
do i=1,n 2425 2425 do i=1,n
res=res+qw(i)*f(xm+xr*qx(i)) 2426 2426 res=res+qw(i)*f(xm+xr*qx(i))
end do 2427 2427 end do
2428 2428
res=xr*res 2429 2429 res=xr*res
2430 2430
deallocate(qw) 2431 2431 deallocate(qw)
deallocate(qx) 2432 2432 deallocate(qx)
2433 2433
end subroutine 2434 2434 end subroutine
2435 2435
!!!!!!!!!!!!!!!!!!!!!!!! 2436 2436 !!!!!!!!!!!!!!!!!!!!!!!!
! 2437 2437 !
! Simple and double adaptative Gauss Kronrod integration based on 2438 2438 ! Simple and double adaptative Gauss Kronrod integration based on
! a modified version of quadpack ( http://www.netlib.org/quadpack 2439 2439 ! a modified version of quadpack ( http://www.netlib.org/quadpack
! 2440 2440 !
! Common parameters : 2441 2441 ! Common parameters :
! 2442 2442 !
! key (in) 2443 2443 ! key (in)
! epsabs 2444 2444 ! epsabs
! epsrel 2445 2445 ! epsrel
! 2446 2446 !
! 2447 2447 !
!!!!!!!!!!!!!!!!!!!!!!!! 2448 2448 !!!!!!!!!!!!!!!!!!!!!!!!
2449 2449
subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit) 2450 2450 subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
! 2451 2451 !
! Evaluate the integral of function f(x) between a and b 2452 2452 ! Evaluate the integral of function f(x) between a and b
! 2453 2453 !
! f(in) : the function 2454 2454 ! f(in) : the function
! a(in) : lower bound 2455 2455 ! a(in) : lower bound
! b(in) : higher bound 2456 2456 ! b(in) : higher bound
! epsabs(in) : desired absolute error 2457 2457 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 2458 2458 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 2459 2459 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 2460 2460 ! 1: 7 - 15 points
! 2: 10 - 21 points 2461 2461 ! 2: 10 - 21 points
! 3: 15 - 31 points 2462 2462 ! 3: 15 - 31 points
! 4: 20 - 41 points 2463 2463 ! 4: 20 - 41 points
! 5: 25 - 51 points 2464 2464 ! 5: 25 - 51 points
! 6: 30 - 61 points 2465 2465 ! 6: 30 - 61 points
! 2466 2466 !
! limit(in) : maximum number of subintervals in the partition of the 2467 2467 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 2468 2468 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 2469 2469 ! behaviour as the imsl routine dqdag
! 2470 2470 !
! res(out) : estimated integral value 2471 2471 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 2472 2472 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 2473 2473 ! ier(out) : error flag from quadpack routines
! 0 : no error 2474 2474 ! 0 : no error
! 1 : maximum number of subdivisions allowed 2475 2475 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 2476 2476 ! has been achieved. one can allow more
! subdivisions by increasing the value of 2477 2477 ! subdivisions by increasing the value of
! limit (and taking the according dimension 2478 2478 ! limit (and taking the according dimension
! adjustments into account). however, if 2479 2479 ! adjustments into account). however, if
! this yield no improvement it is advised 2480 2480 ! this yield no improvement it is advised
! to analyze the integrand in order to 2481 2481 ! to analyze the integrand in order to
! determine the integration difficulaties. 2482 2482 ! determine the integration difficulaties.
! if the position of a local difficulty can 2483 2483 ! if the position of a local difficulty can
! be determined (i.e.singularity, 2484 2484 ! be determined (i.e.singularity,
! discontinuity within the interval) one 2485 2485 ! discontinuity within the interval) one
! will probably gain from splitting up the 2486 2486 ! will probably gain from splitting up the
! interval at this point and calling the 2487 2487 ! interval at this point and calling the
! integrator on the subranges. if possible, 2488 2488 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 2489 2489 ! an appropriate special-purpose integrator
! should be used which is designed for 2490 2490 ! should be used which is designed for
! handling the type of difficulty involved. 2491 2491 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 2492 2492 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 2493 2493 ! detected, which prevents the requested
! tolerance from being achieved. 2494 2494 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 2495 2495 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 2496 2496 ! at some points of the integration
! interval. 2497 2497 ! interval.
! 6 : the input is invalid, because 2498 2498 ! 6 : the input is invalid, because
! (epsabs.le.0 and 2499 2499 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 2500 2500 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 2501 2501 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 2502 2502 ! result, abserr, neval, last are set
! to zero. 2503 2503 ! to zero.
! except when lenw is invalid, iwork(1), 2504 2504 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 2505 2505 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 2506 2506 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 2507 2507 ! work(limit+1) to b.
2508 2508
implicit none 2509 2509 implicit none
double precision, external :: f 2510 2510 double precision, external :: f
double precision, intent(in) :: a,b,epsabs,epsrel 2511 2511 double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key 2512 2512 integer, intent(in) :: key
integer, intent(in) :: limit 2513 2513 integer, intent(in) :: limit
double precision, intent(out) :: res,abserr 2514 2514 double precision, intent(out) :: res,abserr
integer, intent(out) :: ier 2515 2515 integer, intent(out) :: ier
2516 2516
double precision, allocatable :: work(:) 2517 2517 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 2518 2518 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 2519 2519 integer :: lenw,neval,last
2520 2520
! imsl value for limit is 500 2521 2521 ! imsl value for limit is 500
lenw=limit*4 2522 2522 lenw=limit*4
2523 2523
allocate(iwork(limit)) 2524 2524 allocate(iwork(limit))
allocate(work(lenw)) 2525 2525 allocate(work(lenw))
2526 2526
call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 2527 2527 call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
2528 2528
deallocate(work) 2529 2529 deallocate(work)
deallocate(iwork) 2530 2530 deallocate(iwork)
2531 2531
end subroutine 2532 2532 end subroutine
2533 2533
2534 2534
2535 2535
subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit) 2536 2536 subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
! 2537 2537 !
! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x) 2538 2538 ! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)
! 2539 2539 !
! f(in) : the function 2540 2540 ! f(in) : the function
! a(in) : lower bound 2541 2541 ! a(in) : lower bound
! b(in) : higher bound 2542 2542 ! b(in) : higher bound
! g(in) : external function describing lower bound for y 2543 2543 ! g(in) : external function describing lower bound for y
! h(in) : external function describing higher bound for y 2544 2544 ! h(in) : external function describing higher bound for y
! epsabs(in) : desired absolute error 2545 2545 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 2546 2546 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 2547 2547 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 2548 2548 ! 1: 7 - 15 points
! 2: 10 - 21 points 2549 2549 ! 2: 10 - 21 points
! 3: 15 - 31 points 2550 2550 ! 3: 15 - 31 points
! 4: 20 - 41 points 2551 2551 ! 4: 20 - 41 points
! 5: 25 - 51 points 2552 2552 ! 5: 25 - 51 points
! 6: 30 - 61 points 2553 2553 ! 6: 30 - 61 points
! 2554 2554 !
! limit(in) : maximum number of subintervals in the partition of the 2555 2555 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 2556 2556 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 2557 2557 ! behaviour as the imsl routine dqdag
! 2558 2558 !
! res(out) : estimated integral value 2559 2559 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 2560 2560 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 2561 2561 ! ier(out) : error flag from quadpack routines
! 0 : no error 2562 2562 ! 0 : no error
! 1 : maximum number of subdivisions allowed 2563 2563 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 2564 2564 ! has been achieved. one can allow more
! subdivisions by increasing the value of 2565 2565 ! subdivisions by increasing the value of
! limit (and taking the according dimension 2566 2566 ! limit (and taking the according dimension
! adjustments into account). however, if 2567 2567 ! adjustments into account). however, if
! this yield no improvement it is advised 2568 2568 ! this yield no improvement it is advised
! to analyze the integrand in order to 2569 2569 ! to analyze the integrand in order to
! determine the integration difficulaties. 2570 2570 ! determine the integration difficulaties.
! if the position of a local difficulty can 2571 2571 ! if the position of a local difficulty can
! be determined (i.e.singularity, 2572 2572 ! be determined (i.e.singularity,
! discontinuity within the interval) one 2573 2573 ! discontinuity within the interval) one
! will probably gain from splitting up the 2574 2574 ! will probably gain from splitting up the
! interval at this point and calling the 2575 2575 ! interval at this point and calling the
! integrator on the subranges. if possible, 2576 2576 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 2577 2577 ! an appropriate special-purpose integrator
! should be used which is designed for 2578 2578 ! should be used which is designed for
! handling the type of difficulty involved. 2579 2579 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 2580 2580 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 2581 2581 ! detected, which prevents the requested
! tolerance from being achieved. 2582 2582 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 2583 2583 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 2584 2584 ! at some points of the integration
! interval. 2585 2585 ! interval.
! 6 : the input is invalid, because 2586 2586 ! 6 : the input is invalid, because
! (epsabs.le.0 and 2587 2587 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 2588 2588 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 2589 2589 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 2590 2590 ! result, abserr, neval, last are set
! to zero. 2591 2591 ! to zero.
! except when lenw is invalid, iwork(1), 2592 2592 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 2593 2593 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 2594 2594 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 2595 2595 ! work(limit+1) to b.
2596 2596
implicit none 2597 2597 implicit none
double precision, external:: f,g,h 2598 2598 double precision, external:: f,g,h
double precision, intent(in) :: a,b,epsabs,epsrel 2599 2599 double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key,limit 2600 2600 integer, intent(in) :: key,limit
integer, intent(out) :: ier 2601 2601 integer, intent(out) :: ier
double precision, intent(out) :: res,abserr 2602 2602 double precision, intent(out) :: res,abserr
2603 2603
2604 2604
double precision, allocatable :: work(:) 2605 2605 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 2606 2606 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 2607 2607 integer :: lenw,neval,last
2608 2608
! imsl value for limit is 500 2609 2609 ! imsl value for limit is 500
lenw=limit*4 2610 2610 lenw=limit*4
allocate(work(lenw)) 2611 2611 allocate(work(lenw))
allocate(iwork(limit)) 2612 2612 allocate(iwork(limit))
2613 2613
call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 2614 2614 call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
2615 2615
deallocate(iwork) 2616 2616 deallocate(iwork)
deallocate(work) 2617 2617 deallocate(work)
end subroutine 2618 2618 end subroutine
2619 2619
2620 2620
2621 2621
subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit) 2622 2622 subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
! 2623 2623 !
! Evaluate the single integral of function f(x,y) for y between a and b with a 2624 2624 ! Evaluate the single integral of function f(x,y) for y between a and b with a
! given x value 2625 2625 ! given x value
! 2626 2626 !
! This function is used for the evaluation of the double integral fvn_d_integ_2_gk 2627 2627 ! This function is used for the evaluation of the double integral fvn_d_integ_2_gk
! 2628 2628 !
! f(in) : the function 2629 2629 ! f(in) : the function
! x(in) : x 2630 2630 ! x(in) : x
! a(in) : lower bound 2631 2631 ! a(in) : lower bound
! b(in) : higher bound 2632 2632 ! b(in) : higher bound
! epsabs(in) : desired absolute error 2633 2633 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 2634 2634 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 2635 2635 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 2636 2636 ! 1: 7 - 15 points
! 2: 10 - 21 points 2637 2637 ! 2: 10 - 21 points
! 3: 15 - 31 points 2638 2638 ! 3: 15 - 31 points
! 4: 20 - 41 points 2639 2639 ! 4: 20 - 41 points
! 5: 25 - 51 points 2640 2640 ! 5: 25 - 51 points
! 6: 30 - 61 points 2641 2641 ! 6: 30 - 61 points
! 2642 2642 !
! limit(in) : maximum number of subintervals in the partition of the 2643 2643 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 2644 2644 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 2645 2645 ! behaviour as the imsl routine dqdag
! 2646 2646 !
! res(out) : estimated integral value 2647 2647 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 2648 2648 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 2649 2649 ! ier(out) : error flag from quadpack routines
! 0 : no error 2650 2650 ! 0 : no error
! 1 : maximum number of subdivisions allowed 2651 2651 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 2652 2652 ! has been achieved. one can allow more
! subdivisions by increasing the value of 2653 2653 ! subdivisions by increasing the value of
! limit (and taking the according dimension 2654 2654 ! limit (and taking the according dimension
! adjustments into account). however, if 2655 2655 ! adjustments into account). however, if
! this yield no improvement it is advised 2656 2656 ! this yield no improvement it is advised
! to analyze the integrand in order to 2657 2657 ! to analyze the integrand in order to
! determine the integration difficulaties. 2658 2658 ! determine the integration difficulaties.
! if the position of a local difficulty can 2659 2659 ! if the position of a local difficulty can
! be determined (i.e.singularity, 2660 2660 ! be determined (i.e.singularity,
! discontinuity within the interval) one 2661 2661 ! discontinuity within the interval) one
! will probably gain from splitting up the 2662 2662 ! will probably gain from splitting up the
! interval at this point and calling the 2663 2663 ! interval at this point and calling the
! integrator on the subranges. if possible, 2664 2664 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 2665 2665 ! an appropriate special-purpose integrator
! should be used which is designed for 2666 2666 ! should be used which is designed for
! handling the type of difficulty involved. 2667 2667 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 2668 2668 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 2669 2669 ! detected, which prevents the requested
! tolerance from being achieved. 2670 2670 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 2671 2671 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 2672 2672 ! at some points of the integration
! interval. 2673 2673 ! interval.
! 6 : the input is invalid, because 2674 2674 ! 6 : the input is invalid, because
! (epsabs.le.0 and 2675 2675 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 2676 2676 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 2677 2677 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 2678 2678 ! result, abserr, neval, last are set
! to zero. 2679 2679 ! to zero.
! except when lenw is invalid, iwork(1), 2680 2680 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 2681 2681 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 2682 2682 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 2683 2683 ! work(limit+1) to b.
2684 2684
implicit none 2685 2685 implicit none
double precision, external:: f 2686 2686 double precision, external:: f
double precision, intent(in) :: x,a,b,epsabs,epsrel 2687 2687 double precision, intent(in) :: x,a,b,epsabs,epsrel
integer, intent(in) :: key,limit 2688 2688 integer, intent(in) :: key,limit
integer, intent(out) :: ier 2689 2689 integer, intent(out) :: ier
double precision, intent(out) :: res,abserr 2690 2690 double precision, intent(out) :: res,abserr
2691 2691
2692 2692
double precision, allocatable :: work(:) 2693 2693 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 2694 2694 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 2695 2695 integer :: lenw,neval,last
2696 2696
! imsl value for limit is 500 2697 2697 ! imsl value for limit is 500
lenw=limit*4 2698 2698 lenw=limit*4
allocate(work(lenw)) 2699 2699 allocate(work(lenw))
allocate(iwork(limit)) 2700 2700 allocate(iwork(limit))
2701 2701
call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 2702 2702 call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
2703 2703
deallocate(iwork) 2704 2704 deallocate(iwork)
deallocate(work) 2705 2705 deallocate(work)
end subroutine 2706 2706 end subroutine
2707 2707
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2708 2708 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Include the modified quadpack files 2709 2709 ! Include the modified quadpack files
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2710 2710 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
include "fvn_quadpack/dqag_2d_inner.f" 2711 2711 include "fvn_quadpack/dqag_2d_inner.f"
include "fvn_quadpack/dqk15_2d_inner.f" 2712 2712 include "fvn_quadpack/dqk15_2d_inner.f"
include "fvn_quadpack/dqk31_2d_outer.f" 2713 2713 include "fvn_quadpack/dqk31_2d_outer.f"
include "fvn_quadpack/d1mach.f" 2714 2714 include "fvn_quadpack/d1mach.f"
include "fvn_quadpack/dqk31_2d_inner.f" 2715 2715 include "fvn_quadpack/dqk31_2d_inner.f"
include "fvn_quadpack/dqage.f" 2716 2716 include "fvn_quadpack/dqage.f"
include "fvn_quadpack/dqk15.f" 2717 2717 include "fvn_quadpack/dqk15.f"
include "fvn_quadpack/dqk21.f" 2718 2718 include "fvn_quadpack/dqk21.f"
include "fvn_quadpack/dqk31.f" 2719 2719 include "fvn_quadpack/dqk31.f"
include "fvn_quadpack/dqk41.f" 2720 2720 include "fvn_quadpack/dqk41.f"
include "fvn_quadpack/dqk51.f" 2721 2721 include "fvn_quadpack/dqk51.f"
include "fvn_quadpack/dqk61.f" 2722 2722 include "fvn_quadpack/dqk61.f"
include "fvn_quadpack/dqk41_2d_outer.f" 2723 2723 include "fvn_quadpack/dqk41_2d_outer.f"
include "fvn_quadpack/dqk41_2d_inner.f" 2724 2724 include "fvn_quadpack/dqk41_2d_inner.f"
include "fvn_quadpack/dqag_2d_outer.f" 2725 2725 include "fvn_quadpack/dqag_2d_outer.f"
include "fvn_quadpack/dqpsrt.f" 2726 2726 include "fvn_quadpack/dqpsrt.f"
include "fvn_quadpack/dqag.f" 2727 2727 include "fvn_quadpack/dqag.f"
include "fvn_quadpack/dqage_2d_outer.f" 2728 2728 include "fvn_quadpack/dqage_2d_outer.f"
include "fvn_quadpack/dqage_2d_inner.f" 2729 2729 include "fvn_quadpack/dqage_2d_inner.f"
include "fvn_quadpack/dqk51_2d_outer.f" 2730 2730 include "fvn_quadpack/dqk51_2d_outer.f"
include "fvn_quadpack/dqk51_2d_inner.f" 2731 2731 include "fvn_quadpack/dqk51_2d_inner.f"
include "fvn_quadpack/dqk61_2d_outer.f" 2732 2732 include "fvn_quadpack/dqk61_2d_outer.f"
include "fvn_quadpack/dqk21_2d_outer.f" 2733 2733 include "fvn_quadpack/dqk21_2d_outer.f"
include "fvn_quadpack/dqk61_2d_inner.f" 2734 2734 include "fvn_quadpack/dqk61_2d_inner.f"
include "fvn_quadpack/dqk21_2d_inner.f" 2735 2735 include "fvn_quadpack/dqk21_2d_inner.f"
include "fvn_quadpack/dqk15_2d_outer.f" 2736 2736 include "fvn_quadpack/dqk15_2d_outer.f"
2737 2737
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2738 2738 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 2739 2739 !
! Trigonometric functions 2740 2740 ! Trigonometric functions
! 2741 2741 !
! fvn_z_acos, fvn_z_asin : complex arc cosine and sine 2742 2742 ! fvn_z_acos, fvn_z_asin : complex arc cosine and sine
! fvn_d_acosh : arc cosinus hyperbolic 2743 2743 ! fvn_d_acosh : arc cosinus hyperbolic
! fvn_d_asinh : arc sinus hyperbolic 2744 2744 ! fvn_d_asinh : arc sinus hyperbolic
! 2745 2745 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2746 2746 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! February 2008 2747 2747 ! February 2008
! All Trigonometric functions removed due to implementation of fnlib 2748 2748 ! All Trigonometric functions removed due to implementation of fnlib
2749 2749
function fvn_z_acos(z) 2750 2750 function fvn_z_acos(z)
! double complex arccos function adapted from 2751 2751 ! double complex arccos function adapted from
! the c gsl library 2752 2752 ! the c gsl library
! http://www.gnu.org/software/gsl/ 2753 2753 ! http://www.gnu.org/software/gsl/
implicit none 2754 2754 implicit none
complex(kind=8) :: fvn_z_acos 2755 2755 complex(kind=8) :: fvn_z_acos
complex(kind=8) :: z 2756 2756 complex(kind=8) :: z
real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1 2757 2757 real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8 2758 2758 real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
complex(kind=8),parameter :: i=(0._8,1._8) 2759 2759 complex(kind=8),parameter :: i=(0._8,1._8)
real(kind=8) :: r_res,i_res 2760 2760 real(kind=8) :: r_res,i_res
2761 2761
rz=dreal(z) 2762 2762 rz=dreal(z)
iz=dimag(z) 2763 2763 iz=dimag(z)
if ( iz == 0._8 ) then 2764 2764 if ( iz == 0._8 ) then
fvn_z_acos=fvn_z_acos_real(rz) 2765 2765 fvn_z_acos=fvn_z_acos_real(rz)
return 2766 2766 return
end if 2767 2767 end if
2768 2768
x=dabs(rz) 2769 2769 x=dabs(rz)
y=dabs(iz) 2770 2770 y=dabs(iz)
r=fvn_d_hypot(x+1.,y) 2771 2771 r=fvn_d_hypot(x+1.,y)
s=fvn_d_hypot(x-1.,y) 2772 2772 s=fvn_d_hypot(x-1.,y)
a=0.5*(r + s) 2773 2773 a=0.5*(r + s)
b=x/a 2774 2774 b=x/a
y2=y*y 2775 2775 y2=y*y
2776 2776
if (b <= b_crossover) then 2777 2777 if (b <= b_crossover) then
r_res=dacos(b) 2778 2778 r_res=dacos(b)
else 2779 2779 else
if (x <= 1.) then 2780 2780 if (x <= 1.) then
d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x))) 2781 2781 d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
r_res=datan(dsqrt(d)/x) 2782 2782 r_res=datan(dsqrt(d)/x)
else 2783 2783 else
apx=a+x 2784 2784 apx=a+x
d=0.5*(apx/(r+x+1)+apx/(s + (x - 1))) 2785 2785 d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
r_res=datan((y*dsqrt(d))/x); 2786 2786 r_res=datan((y*dsqrt(d))/x);
end if 2787 2787 end if
end if 2788 2788 end if
2789 2789
if (a <= a_crossover) then 2790 2790 if (a <= a_crossover) then
if (x < 1.) then 2791 2791 if (x < 1.) then
am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x))) 2792 2792 am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
else 2793 2793 else
am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1))) 2794 2794 am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
end if 2795 2795 end if
i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1)))); 2796 2796 i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
else 2797 2797 else
i_res = dlog (a + dsqrt (a*a - 1.)); 2798 2798 i_res = dlog (a + dsqrt (a*a - 1.));
end if 2799 2799 end if
if (rz <0.) then 2800 2800 if (rz <0.) then
r_res=fvn_pi-r_res 2801 2801 r_res=fvn_pi-r_res
end if 2802 2802 end if
i_res=-sign(1._8,iz)*i_res 2803 2803 i_res=-sign(1._8,iz)*i_res
fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res) 2804 2804 fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)
2805 2805
end function fvn_z_acos 2806 2806 end function fvn_z_acos
2807 2807
function fvn_z_acos_real(r) 2808 2808 function fvn_z_acos_real(r)
! return the double complex arc cosinus for a 2809 2809 ! return the double complex arc cosinus for a
! double precision argument 2810 2810 ! double precision argument
implicit none 2811 2811 implicit none
real(kind=8) :: r 2812 2812 real(kind=8) :: r
complex(kind=8) :: fvn_z_acos_real 2813 2813 complex(kind=8) :: fvn_z_acos_real
2814 2814
if (dabs(r)<=1._8) then 2815 2815 if (dabs(r)<=1._8) then
fvn_z_acos_real=dcmplx(dacos(r)) 2816 2816 fvn_z_acos_real=dcmplx(dacos(r))
return 2817 2817 return
end if 2818 2818 end if
if (r < 0._8) then 2819 2819 if (r < 0._8) then
fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r)) 2820 2820 fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))
else 2821 2821 else
fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r)) 2822 2822 fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))
end if 2823 2823 end if
end function 2824 2824 end function
2825 2825
2826 2826
function fvn_z_asin(z) 2827 2827 function fvn_z_asin(z)
! double complex arcsin function derived from 2828 2828 ! double complex arcsin function derived from
! the c gsl library 2829 2829 ! the c gsl library
! http://www.gnu.org/software/gsl/ 2830 2830 ! http://www.gnu.org/software/gsl/
implicit none 2831 2831 implicit none
complex(kind=8) :: fvn_z_asin 2832 2832 complex(kind=8) :: fvn_z_asin
complex(kind=8) :: z 2833 2833 complex(kind=8) :: z
real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1 2834 2834 real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1
real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8 2835 2835 real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
real(kind=8) :: r_res,i_res 2836 2836 real(kind=8) :: r_res,i_res
2837 2837
rz=dreal(z) 2838 2838 rz=dreal(z)
iz=dimag(z) 2839 2839 iz=dimag(z)
if ( iz == 0._8 ) then 2840 2840 if ( iz == 0._8 ) then
! z is real 2841 2841 ! z is real
fvn_z_asin=fvn_z_asin_real(rz) 2842 2842 fvn_z_asin=fvn_z_asin_real(rz)
return 2843 2843 return
end if 2844 2844 end if
2845 2845
x=dabs(rz) 2846 2846 x=dabs(rz)
y=dabs(iz) 2847 2847 y=dabs(iz)
r=fvn_d_hypot(x+1.,y) 2848 2848 r=fvn_d_hypot(x+1.,y)
s=fvn_d_hypot(x-1.,y) 2849 2849 s=fvn_d_hypot(x-1.,y)
a=0.5*(r + s) 2850 2850 a=0.5*(r + s)
b=x/a 2851 2851 b=x/a
y2=y*y 2852 2852 y2=y*y
2853 2853
if (b <= b_crossover) then 2854 2854 if (b <= b_crossover) then
r_res=dasin(b) 2855 2855 r_res=dasin(b)
else 2856 2856 else
if (x <= 1.) then 2857 2857 if (x <= 1.) then
d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x))) 2858 2858 d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
r_res=datan(x/dsqrt(d)) 2859 2859 r_res=datan(x/dsqrt(d))
else 2860 2860 else
apx=a+x 2861 2861 apx=a+x
d=0.5*(apx/(r+x+1)+apx/(s + (x - 1))) 2862 2862 d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
r_res=datan(x/(y*dsqrt(d))); 2863 2863 r_res=datan(x/(y*dsqrt(d)));
end if 2864 2864 end if
end if 2865 2865 end if
2866 2866
if (a <= a_crossover) then 2867 2867 if (a <= a_crossover) then
if (x < 1.) then 2868 2868 if (x < 1.) then
am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x))) 2869 2869 am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
else 2870 2870 else
am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1))) 2871 2871 am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
end if 2872 2872 end if
i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1)))); 2873 2873 i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
else 2874 2874 else
i_res = dlog (a + dsqrt (a*a - 1.)); 2875 2875 i_res = dlog (a + dsqrt (a*a - 1.));
end if 2876 2876 end if
r_res=sign(1._8,rz)*r_res 2877 2877 r_res=sign(1._8,rz)*r_res
i_res=sign(1._8,iz)*i_res 2878 2878 i_res=sign(1._8,iz)*i_res
fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res) 2879 2879 fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)
2880 2880
end function fvn_z_asin 2881 2881 end function fvn_z_asin
2882 2882
function fvn_z_asin_real(r) 2883 2883 function fvn_z_asin_real(r)
! return the double complex arc sinus for a 2884 2884 ! return the double complex arc sinus for a
! double precision argument 2885 2885 ! double precision argument
implicit none 2886 2886 implicit none
real(kind=8) :: r 2887 2887 real(kind=8) :: r
complex(kind=8) :: fvn_z_asin_real 2888 2888 complex(kind=8) :: fvn_z_asin_real
2889 2889
if (dabs(r)<=1._8) then 2890 2890 if (dabs(r)<=1._8) then
fvn_z_asin_real=dcmplx(dasin(r)) 2891 2891 fvn_z_asin_real=dcmplx(dasin(r))
return 2892 2892 return
end if 2893 2893 end if
if (r < 0._8) then 2894 2894 if (r < 0._8) then
fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r)) 2895 2895 fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))
else 2896 2896 else
fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r)) 2897 2897 fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))
end if 2898 2898 end if
end function fvn_z_asin_real 2899 2899 end function fvn_z_asin_real
2900 2900
function fvn_d_acosh(r) 2901 2901 function fvn_d_acosh(r)
! return the arc hyperbolic cosine 2902 2902 ! return the arc hyperbolic cosine
implicit none 2903 2903 implicit none
real(kind=8) :: r 2904 2904 real(kind=8) :: r
real(kind=8) :: fvn_d_acosh 2905 2905 real(kind=8) :: fvn_d_acosh
if (r >=1) then 2906 2906 if (r >=1) then
fvn_d_acosh=dlog(r+dsqrt(r*r-1)) 2907 2907 fvn_d_acosh=dlog(r+dsqrt(r*r-1))
else 2908 2908 else
!! TODO : Better error handling!!!!!! 2909 2909 !! TODO : Better error handling!!!!!!