lxsd / fvn

Download as
Browse Code ยป

Commit ea0bc1fe357570ea079fa56adbe6566c7028b797

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

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

Showing 1 changed file with 2 additions and 1 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 9 9 ! William Daniau 2007
! 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 ;-) Nevertheless, you
23 ! may give credits to quadpack authors.
! 23 24 !
! Version 1.1 24 25 ! Version 1.1
! 25 26 !
! TO DO LIST : 26 27 ! TO DO LIST :
! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm 27 28 ! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm
! eigenvalues are given with no particular order. 28 29 ! eigenvalues are given with no particular order.
! + Generic interface for fvn_x_name family -> fvn_name 29 30 ! + Generic interface for fvn_x_name family -> fvn_name
! + Make some parameters optional, status for example 30 31 ! + Make some parameters optional, status for example
! + use f95 kinds "double complex" -> complex(kind=8) 31 32 ! + use f95 kinds "double complex" -> complex(kind=8)
! + unify quadpack routines 32 33 ! + unify quadpack routines
! + ... 33 34 ! + ...
! 34 35 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 35 36 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
36 37
implicit none 37 38 implicit none
! All quadpack routines are private to the module 38 39 ! All quadpack routines are private to the module
private :: d1mach,dqag,dqag_2d_inner,dqag_2d_outer,dqage,dqage_2d_inner, & 39 40 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, & 40 41 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, & 41 42 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 42 43 dqk51,dqk51_2d_inner,dqk51_2d_outer,dqk61,dqk61_2d_inner,dqk61_2d_outer,dqpsrt
43 44
44 45
contains 45 46 contains
46 47
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 47 48 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 48 49 !
! Matrix inversion subroutines 49 50 ! Matrix inversion subroutines
! 50 51 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 51 52 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine fvn_s_matinv(d,a,inva,status) 52 53 subroutine fvn_s_matinv(d,a,inva,status)
! 53 54 !
! Matrix inversion of a real matrix using BLAS and LAPACK 54 55 ! Matrix inversion of a real matrix using BLAS and LAPACK
! 55 56 !
! d (in) : matrix rank 56 57 ! d (in) : matrix rank
! a (in) : input matrix 57 58 ! a (in) : input matrix
! inva (out) : inversed matrix 58 59 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 59 60 ! status (ou) : =0 if something failed
! 60 61 !
integer, intent(in) :: d 61 62 integer, intent(in) :: d
real, intent(in) :: a(d,d) 62 63 real, intent(in) :: a(d,d)
real, intent(out) :: inva(d,d) 63 64 real, intent(out) :: inva(d,d)
integer, intent(out) :: status 64 65 integer, intent(out) :: status
65 66
integer, allocatable :: ipiv(:) 66 67 integer, allocatable :: ipiv(:)
real, allocatable :: work(:) 67 68 real, allocatable :: work(:)
real twork(1) 68 69 real twork(1)
integer :: info 69 70 integer :: info
integer :: lwork 70 71 integer :: lwork
71 72
status=1 72 73 status=1
73 74
allocate(ipiv(d)) 74 75 allocate(ipiv(d))
! copy a into inva using BLAS 75 76 ! copy a into inva using BLAS
!call scopy(d*d,a,1,inva,1) 76 77 !call scopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 77 78 inva(:,:)=a(:,:)
! LU factorization using LAPACK 78 79 ! LU factorization using LAPACK
call sgetrf(d,d,inva,d,ipiv,info) 79 80 call sgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 80 81 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 81 82 if (info /= 0) then
status=0 82 83 status=0
deallocate(ipiv) 83 84 deallocate(ipiv)
return 84 85 return
end if 85 86 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 86 87 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call sgetri(d,inva,d,ipiv,twork,-1,info) 87 88 call sgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 88 89 lwork=int(twork(1))
allocate(work(lwork)) 89 90 allocate(work(lwork))
! Matrix inversion using LAPACK 90 91 ! Matrix inversion using LAPACK
call sgetri(d,inva,d,ipiv,work,lwork,info) 91 92 call sgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 92 93 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 93 94 if (info /= 0) then
status=0 94 95 status=0
end if 95 96 end if
deallocate(work) 96 97 deallocate(work)
deallocate(ipiv) 97 98 deallocate(ipiv)
end subroutine 98 99 end subroutine
99 100
subroutine fvn_d_matinv(d,a,inva,status) 100 101 subroutine fvn_d_matinv(d,a,inva,status)
! 101 102 !
! Matrix inversion of a double precision matrix using BLAS and LAPACK 102 103 ! Matrix inversion of a double precision matrix using BLAS and LAPACK
! 103 104 !
! d (in) : matrix rank 104 105 ! d (in) : matrix rank
! a (in) : input matrix 105 106 ! a (in) : input matrix
! inva (out) : inversed matrix 106 107 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 107 108 ! status (ou) : =0 if something failed
! 108 109 !
integer, intent(in) :: d 109 110 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 110 111 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: inva(d,d) 111 112 double precision, intent(out) :: inva(d,d)
integer, intent(out) :: status 112 113 integer, intent(out) :: status
113 114
integer, allocatable :: ipiv(:) 114 115 integer, allocatable :: ipiv(:)
double precision, allocatable :: work(:) 115 116 double precision, allocatable :: work(:)
double precision :: twork(1) 116 117 double precision :: twork(1)
integer :: info 117 118 integer :: info
integer :: lwork 118 119 integer :: lwork
119 120
status=1 120 121 status=1
121 122
allocate(ipiv(d)) 122 123 allocate(ipiv(d))
! copy a into inva using BLAS 123 124 ! copy a into inva using BLAS
!call dcopy(d*d,a,1,inva,1) 124 125 !call dcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 125 126 inva(:,:)=a(:,:)
! LU factorization using LAPACK 126 127 ! LU factorization using LAPACK
call dgetrf(d,d,inva,d,ipiv,info) 127 128 call dgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 128 129 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 129 130 if (info /= 0) then
status=0 130 131 status=0
deallocate(ipiv) 131 132 deallocate(ipiv)
return 132 133 return
end if 133 134 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 134 135 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call dgetri(d,inva,d,ipiv,twork,-1,info) 135 136 call dgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 136 137 lwork=int(twork(1))
allocate(work(lwork)) 137 138 allocate(work(lwork))
! Matrix inversion using LAPACK 138 139 ! Matrix inversion using LAPACK
call dgetri(d,inva,d,ipiv,work,lwork,info) 139 140 call dgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 140 141 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 141 142 if (info /= 0) then
status=0 142 143 status=0
end if 143 144 end if
deallocate(work) 144 145 deallocate(work)
deallocate(ipiv) 145 146 deallocate(ipiv)
end subroutine 146 147 end subroutine
147 148
subroutine fvn_c_matinv(d,a,inva,status) 148 149 subroutine fvn_c_matinv(d,a,inva,status)
! 149 150 !
! Matrix inversion of a complex matrix using BLAS and LAPACK 150 151 ! Matrix inversion of a complex matrix using BLAS and LAPACK
! 151 152 !
! d (in) : matrix rank 152 153 ! d (in) : matrix rank
! a (in) : input matrix 153 154 ! a (in) : input matrix
! inva (out) : inversed matrix 154 155 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 155 156 ! status (ou) : =0 if something failed
! 156 157 !
integer, intent(in) :: d 157 158 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 158 159 complex, intent(in) :: a(d,d)
complex, intent(out) :: inva(d,d) 159 160 complex, intent(out) :: inva(d,d)
integer, intent(out) :: status 160 161 integer, intent(out) :: status
161 162
integer, allocatable :: ipiv(:) 162 163 integer, allocatable :: ipiv(:)
complex, allocatable :: work(:) 163 164 complex, allocatable :: work(:)
complex :: twork(1) 164 165 complex :: twork(1)
integer :: info 165 166 integer :: info
integer :: lwork 166 167 integer :: lwork
167 168
status=1 168 169 status=1
169 170
allocate(ipiv(d)) 170 171 allocate(ipiv(d))
! copy a into inva using BLAS 171 172 ! copy a into inva using BLAS
!call ccopy(d*d,a,1,inva,1) 172 173 !call ccopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 173 174 inva(:,:)=a(:,:)
174 175
! LU factorization using LAPACK 175 176 ! LU factorization using LAPACK
call cgetrf(d,d,inva,d,ipiv,info) 176 177 call cgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 177 178 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 178 179 if (info /= 0) then
status=0 179 180 status=0
deallocate(ipiv) 180 181 deallocate(ipiv)
return 181 182 return
end if 182 183 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 183 184 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call cgetri(d,inva,d,ipiv,twork,-1,info) 184 185 call cgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 185 186 lwork=int(twork(1))
allocate(work(lwork)) 186 187 allocate(work(lwork))
! Matrix inversion using LAPACK 187 188 ! Matrix inversion using LAPACK
call cgetri(d,inva,d,ipiv,work,lwork,info) 188 189 call cgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 189 190 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 190 191 if (info /= 0) then
status=0 191 192 status=0
end if 192 193 end if
deallocate(work) 193 194 deallocate(work)
deallocate(ipiv) 194 195 deallocate(ipiv)
end subroutine 195 196 end subroutine
196 197
subroutine fvn_z_matinv(d,a,inva,status) 197 198 subroutine fvn_z_matinv(d,a,inva,status)
! 198 199 !
! Matrix inversion of a double complex matrix using BLAS and LAPACK 199 200 ! Matrix inversion of a double complex matrix using BLAS and LAPACK
! 200 201 !
! d (in) : matrix rank 201 202 ! d (in) : matrix rank
! a (in) : input matrix 202 203 ! a (in) : input matrix
! inva (out) : inversed matrix 203 204 ! inva (out) : inversed matrix
! status (ou) : =0 if something failed 204 205 ! status (ou) : =0 if something failed
! 205 206 !
integer, intent(in) :: d 206 207 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 207 208 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: inva(d,d) 208 209 double complex, intent(out) :: inva(d,d)
integer, intent(out) :: status 209 210 integer, intent(out) :: status
210 211
integer, allocatable :: ipiv(:) 211 212 integer, allocatable :: ipiv(:)
double complex, allocatable :: work(:) 212 213 double complex, allocatable :: work(:)
double complex :: twork(1) 213 214 double complex :: twork(1)
integer :: info 214 215 integer :: info
integer :: lwork 215 216 integer :: lwork
216 217
status=1 217 218 status=1
218 219
allocate(ipiv(d)) 219 220 allocate(ipiv(d))
! copy a into inva using BLAS 220 221 ! copy a into inva using BLAS
!call zcopy(d*d,a,1,inva,1) 221 222 !call zcopy(d*d,a,1,inva,1)
inva(:,:)=a(:,:) 222 223 inva(:,:)=a(:,:)
223 224
! LU factorization using LAPACK 224 225 ! LU factorization using LAPACK
call zgetrf(d,d,inva,d,ipiv,info) 225 226 call zgetrf(d,d,inva,d,ipiv,info)
! if info is not equal to 0, something went wrong we exit setting status to 0 226 227 ! if info is not equal to 0, something went wrong we exit setting status to 0
if (info /= 0) then 227 228 if (info /= 0) then
status=0 228 229 status=0
deallocate(ipiv) 229 230 deallocate(ipiv)
return 230 231 return
end if 231 232 end if
! we use the query fonction of xxxtri to obtain the optimal workspace size 232 233 ! we use the query fonction of xxxtri to obtain the optimal workspace size
call zgetri(d,inva,d,ipiv,twork,-1,info) 233 234 call zgetri(d,inva,d,ipiv,twork,-1,info)
lwork=int(twork(1)) 234 235 lwork=int(twork(1))
allocate(work(lwork)) 235 236 allocate(work(lwork))
! Matrix inversion using LAPACK 236 237 ! Matrix inversion using LAPACK
call zgetri(d,inva,d,ipiv,work,lwork,info) 237 238 call zgetri(d,inva,d,ipiv,work,lwork,info)
! again if info is not equal to 0, we exit setting status to 0 238 239 ! again if info is not equal to 0, we exit setting status to 0
if (info /= 0) then 239 240 if (info /= 0) then
status=0 240 241 status=0
end if 241 242 end if
deallocate(work) 242 243 deallocate(work)
deallocate(ipiv) 243 244 deallocate(ipiv)
end subroutine 244 245 end subroutine
245 246
246 247
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 247 248 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 248 249 !
! Determinants 249 250 ! Determinants
! 250 251 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 251 252 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_s_det(d,a,status) 252 253 function fvn_s_det(d,a,status)
! 253 254 !
! Evaluate the determinant of a square matrix using lapack LU factorization 254 255 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 255 256 !
! d (in) : matrix rank 256 257 ! d (in) : matrix rank
! a (in) : The Matrix 257 258 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 258 259 ! status (out) : =0 if LU factorization failed
! 259 260 !
integer, intent(in) :: d 260 261 integer, intent(in) :: d
real, intent(in) :: a(d,d) 261 262 real, intent(in) :: a(d,d)
integer, intent(out) :: status 262 263 integer, intent(out) :: status
real :: fvn_s_det 263 264 real :: fvn_s_det
264 265
real, allocatable :: wc_a(:,:) 265 266 real, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 266 267 integer, allocatable :: ipiv(:)
integer :: info,i 267 268 integer :: info,i
268 269
status=1 269 270 status=1
allocate(wc_a(d,d)) 270 271 allocate(wc_a(d,d))
allocate(ipiv(d)) 271 272 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 272 273 wc_a(:,:)=a(:,:)
call sgetrf(d,d,wc_a,d,ipiv,info) 273 274 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 274 275 if (info/= 0) then
status=0 275 276 status=0
fvn_s_det=0.e0 276 277 fvn_s_det=0.e0
deallocate(ipiv) 277 278 deallocate(ipiv)
deallocate(wc_a) 278 279 deallocate(wc_a)
return 279 280 return
end if 280 281 end if
fvn_s_det=1.e0 281 282 fvn_s_det=1.e0
do i=1,d 282 283 do i=1,d
if (ipiv(i)==i) then 283 284 if (ipiv(i)==i) then
fvn_s_det=fvn_s_det*wc_a(i,i) 284 285 fvn_s_det=fvn_s_det*wc_a(i,i)
else 285 286 else
fvn_s_det=-fvn_s_det*wc_a(i,i) 286 287 fvn_s_det=-fvn_s_det*wc_a(i,i)
end if 287 288 end if
end do 288 289 end do
deallocate(ipiv) 289 290 deallocate(ipiv)
deallocate(wc_a) 290 291 deallocate(wc_a)
291 292
end function 292 293 end function
293 294
function fvn_d_det(d,a,status) 294 295 function fvn_d_det(d,a,status)
! 295 296 !
! Evaluate the determinant of a square matrix using lapack LU factorization 296 297 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 297 298 !
! d (in) : matrix rank 298 299 ! d (in) : matrix rank
! a (in) : The Matrix 299 300 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 300 301 ! status (out) : =0 if LU factorization failed
! 301 302 !
integer, intent(in) :: d 302 303 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 303 304 double precision, intent(in) :: a(d,d)
integer, intent(out) :: status 304 305 integer, intent(out) :: status
double precision :: fvn_d_det 305 306 double precision :: fvn_d_det
306 307
double precision, allocatable :: wc_a(:,:) 307 308 double precision, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 308 309 integer, allocatable :: ipiv(:)
integer :: info,i 309 310 integer :: info,i
310 311
status=1 311 312 status=1
allocate(wc_a(d,d)) 312 313 allocate(wc_a(d,d))
allocate(ipiv(d)) 313 314 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 314 315 wc_a(:,:)=a(:,:)
call dgetrf(d,d,wc_a,d,ipiv,info) 315 316 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 316 317 if (info/= 0) then
status=0 317 318 status=0
fvn_d_det=0.d0 318 319 fvn_d_det=0.d0
deallocate(ipiv) 319 320 deallocate(ipiv)
deallocate(wc_a) 320 321 deallocate(wc_a)
return 321 322 return
end if 322 323 end if
fvn_d_det=1.d0 323 324 fvn_d_det=1.d0
do i=1,d 324 325 do i=1,d
if (ipiv(i)==i) then 325 326 if (ipiv(i)==i) then
fvn_d_det=fvn_d_det*wc_a(i,i) 326 327 fvn_d_det=fvn_d_det*wc_a(i,i)
else 327 328 else
fvn_d_det=-fvn_d_det*wc_a(i,i) 328 329 fvn_d_det=-fvn_d_det*wc_a(i,i)
end if 329 330 end if
end do 330 331 end do
deallocate(ipiv) 331 332 deallocate(ipiv)
deallocate(wc_a) 332 333 deallocate(wc_a)
333 334
end function 334 335 end function
335 336
function fvn_c_det(d,a,status) ! 336 337 function fvn_c_det(d,a,status) !
! Evaluate the determinant of a square matrix using lapack LU factorization 337 338 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 338 339 !
! d (in) : matrix rank 339 340 ! d (in) : matrix rank
! a (in) : The Matrix 340 341 ! a (in) : The Matrix
! status (out) : =0 if LU factorization failed 341 342 ! status (out) : =0 if LU factorization failed
! 342 343 !
integer, intent(in) :: d 343 344 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 344 345 complex, intent(in) :: a(d,d)
integer, intent(out) :: status 345 346 integer, intent(out) :: status
complex :: fvn_c_det 346 347 complex :: fvn_c_det
347 348
complex, allocatable :: wc_a(:,:) 348 349 complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 349 350 integer, allocatable :: ipiv(:)
integer :: info,i 350 351 integer :: info,i
351 352
status=1 352 353 status=1
allocate(wc_a(d,d)) 353 354 allocate(wc_a(d,d))
allocate(ipiv(d)) 354 355 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 355 356 wc_a(:,:)=a(:,:)
call cgetrf(d,d,wc_a,d,ipiv,info) 356 357 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 357 358 if (info/= 0) then
status=0 358 359 status=0
fvn_c_det=(0.e0,0.e0) 359 360 fvn_c_det=(0.e0,0.e0)
deallocate(ipiv) 360 361 deallocate(ipiv)
deallocate(wc_a) 361 362 deallocate(wc_a)
return 362 363 return
end if 363 364 end if
fvn_c_det=(1.e0,0.e0) 364 365 fvn_c_det=(1.e0,0.e0)
do i=1,d 365 366 do i=1,d
if (ipiv(i)==i) then 366 367 if (ipiv(i)==i) then
fvn_c_det=fvn_c_det*wc_a(i,i) 367 368 fvn_c_det=fvn_c_det*wc_a(i,i)
else 368 369 else
fvn_c_det=-fvn_c_det*wc_a(i,i) 369 370 fvn_c_det=-fvn_c_det*wc_a(i,i)
end if 370 371 end if
end do 371 372 end do
deallocate(ipiv) 372 373 deallocate(ipiv)
deallocate(wc_a) 373 374 deallocate(wc_a)
374 375
end function 375 376 end function
376 377
function fvn_z_det(d,a,status) 377 378 function fvn_z_det(d,a,status)
! 378 379 !
! Evaluate the determinant of a square matrix using lapack LU factorization 379 380 ! Evaluate the determinant of a square matrix using lapack LU factorization
! 380 381 !
! d (in) : matrix rank 381 382 ! d (in) : matrix rank
! a (in) : The Matrix 382 383 ! a (in) : The Matrix
! det (out) : determinant 383 384 ! det (out) : determinant
! status (out) : =0 if LU factorization failed 384 385 ! status (out) : =0 if LU factorization failed
! 385 386 !
integer, intent(in) :: d 386 387 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 387 388 double complex, intent(in) :: a(d,d)
integer, intent(out) :: status 388 389 integer, intent(out) :: status
double complex :: fvn_z_det 389 390 double complex :: fvn_z_det
390 391
double complex, allocatable :: wc_a(:,:) 391 392 double complex, allocatable :: wc_a(:,:)
integer, allocatable :: ipiv(:) 392 393 integer, allocatable :: ipiv(:)
integer :: info,i 393 394 integer :: info,i
394 395
status=1 395 396 status=1
allocate(wc_a(d,d)) 396 397 allocate(wc_a(d,d))
allocate(ipiv(d)) 397 398 allocate(ipiv(d))
wc_a(:,:)=a(:,:) 398 399 wc_a(:,:)=a(:,:)
call zgetrf(d,d,wc_a,d,ipiv,info) 399 400 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info/= 0) then 400 401 if (info/= 0) then
status=0 401 402 status=0
fvn_z_det=(0.d0,0.d0) 402 403 fvn_z_det=(0.d0,0.d0)
deallocate(ipiv) 403 404 deallocate(ipiv)
deallocate(wc_a) 404 405 deallocate(wc_a)
return 405 406 return
end if 406 407 end if
fvn_z_det=(1.d0,0.d0) 407 408 fvn_z_det=(1.d0,0.d0)
do i=1,d 408 409 do i=1,d
if (ipiv(i)==i) then 409 410 if (ipiv(i)==i) then
fvn_z_det=fvn_z_det*wc_a(i,i) 410 411 fvn_z_det=fvn_z_det*wc_a(i,i)
else 411 412 else
fvn_z_det=-fvn_z_det*wc_a(i,i) 412 413 fvn_z_det=-fvn_z_det*wc_a(i,i)
end if 413 414 end if
end do 414 415 end do
deallocate(ipiv) 415 416 deallocate(ipiv)
deallocate(wc_a) 416 417 deallocate(wc_a)
417 418
end function 418 419 end function
419 420
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 420 421 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 421 422 !
! Condition test 422 423 ! Condition test
! 423 424 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 424 425 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1-norm 425 426 ! 1-norm
! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm 426 427 ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm
! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond 427 428 ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond
! 428 429 !
subroutine fvn_s_matcon(d,a,rcond,status) 429 430 subroutine fvn_s_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 430 431 ! Matrix condition (reciprocal of condition number)
! 431 432 !
! d (in) : matrix rank 432 433 ! d (in) : matrix rank
! a (in) : The Matrix 433 434 ! a (in) : The Matrix
! rcond (out) : guess what 434 435 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 435 436 ! status (out) : =0 if something went wrong
! 436 437 !
integer, intent(in) :: d 437 438 integer, intent(in) :: d
real, intent(in) :: a(d,d) 438 439 real, intent(in) :: a(d,d)
real, intent(out) :: rcond 439 440 real, intent(out) :: rcond
integer, intent(out) :: status 440 441 integer, intent(out) :: status
441 442
real, allocatable :: work(:) 442 443 real, allocatable :: work(:)
integer, allocatable :: iwork(:) 443 444 integer, allocatable :: iwork(:)
real :: anorm 444 445 real :: anorm
real, allocatable :: wc_a(:,:) ! working copy of a 445 446 real, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 446 447 integer :: info
integer, allocatable :: ipiv(:) 447 448 integer, allocatable :: ipiv(:)
448 449
real, external :: slange 449 450 real, external :: slange
450 451
451 452
status=1 452 453 status=1
453 454
anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 454 455 anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
455 456
allocate(wc_a(d,d)) 456 457 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 457 458 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 458 459 wc_a(:,:)=a(:,:)
459 460
allocate(ipiv(d)) 460 461 allocate(ipiv(d))
call sgetrf(d,d,wc_a,d,ipiv,info) 461 462 call sgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 462 463 if (info /= 0) then
status=0 463 464 status=0
deallocate(ipiv) 464 465 deallocate(ipiv)
deallocate(wc_a) 465 466 deallocate(wc_a)
return 466 467 return
end if 467 468 end if
allocate(work(4*d)) 468 469 allocate(work(4*d))
allocate(iwork(d)) 469 470 allocate(iwork(d))
call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 470 471 call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 471 472 if (info /= 0) then
status=0 472 473 status=0
end if 473 474 end if
deallocate(iwork) 474 475 deallocate(iwork)
deallocate(work) 475 476 deallocate(work)
deallocate(ipiv) 476 477 deallocate(ipiv)
deallocate(wc_a) 477 478 deallocate(wc_a)
478 479
end subroutine 479 480 end subroutine
480 481
subroutine fvn_d_matcon(d,a,rcond,status) 481 482 subroutine fvn_d_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 482 483 ! Matrix condition (reciprocal of condition number)
! 483 484 !
! d (in) : matrix rank 484 485 ! d (in) : matrix rank
! a (in) : The Matrix 485 486 ! a (in) : The Matrix
! rcond (out) : guess what 486 487 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 487 488 ! status (out) : =0 if something went wrong
! 488 489 !
integer, intent(in) :: d 489 490 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 490 491 double precision, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 491 492 double precision, intent(out) :: rcond
integer, intent(out) :: status 492 493 integer, intent(out) :: status
493 494
double precision, allocatable :: work(:) 494 495 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 495 496 integer, allocatable :: iwork(:)
double precision :: anorm 496 497 double precision :: anorm
double precision, allocatable :: wc_a(:,:) ! working copy of a 497 498 double precision, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 498 499 integer :: info
integer, allocatable :: ipiv(:) 499 500 integer, allocatable :: ipiv(:)
500 501
double precision, external :: dlange 501 502 double precision, external :: dlange
502 503
503 504
status=1 504 505 status=1
505 506
anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm 506 507 anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm
507 508
allocate(wc_a(d,d)) 508 509 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 509 510 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 510 511 wc_a(:,:)=a(:,:)
511 512
allocate(ipiv(d)) 512 513 allocate(ipiv(d))
call dgetrf(d,d,wc_a,d,ipiv,info) 513 514 call dgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 514 515 if (info /= 0) then
status=0 515 516 status=0
deallocate(ipiv) 516 517 deallocate(ipiv)
deallocate(wc_a) 517 518 deallocate(wc_a)
return 518 519 return
end if 519 520 end if
520 521
allocate(work(4*d)) 521 522 allocate(work(4*d))
allocate(iwork(d)) 522 523 allocate(iwork(d))
call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) 523 524 call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info)
if (info /= 0) then 524 525 if (info /= 0) then
status=0 525 526 status=0
end if 526 527 end if
deallocate(iwork) 527 528 deallocate(iwork)
deallocate(work) 528 529 deallocate(work)
deallocate(ipiv) 529 530 deallocate(ipiv)
deallocate(wc_a) 530 531 deallocate(wc_a)
531 532
end subroutine 532 533 end subroutine
533 534
subroutine fvn_c_matcon(d,a,rcond,status) 534 535 subroutine fvn_c_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 535 536 ! Matrix condition (reciprocal of condition number)
! 536 537 !
! d (in) : matrix rank 537 538 ! d (in) : matrix rank
! a (in) : The Matrix 538 539 ! a (in) : The Matrix
! rcond (out) : guess what 539 540 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 540 541 ! status (out) : =0 if something went wrong
! 541 542 !
integer, intent(in) :: d 542 543 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 543 544 complex, intent(in) :: a(d,d)
real, intent(out) :: rcond 544 545 real, intent(out) :: rcond
integer, intent(out) :: status 545 546 integer, intent(out) :: status
546 547
real, allocatable :: rwork(:) 547 548 real, allocatable :: rwork(:)
complex, allocatable :: work(:) 548 549 complex, allocatable :: work(:)
integer, allocatable :: iwork(:) 549 550 integer, allocatable :: iwork(:)
real :: anorm 550 551 real :: anorm
complex, allocatable :: wc_a(:,:) ! working copy of a 551 552 complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 552 553 integer :: info
integer, allocatable :: ipiv(:) 553 554 integer, allocatable :: ipiv(:)
554 555
real, external :: clange 555 556 real, external :: clange
556 557
557 558
status=1 558 559 status=1
559 560
anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 560 561 anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
561 562
allocate(wc_a(d,d)) 562 563 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 563 564 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 564 565 wc_a(:,:)=a(:,:)
565 566
allocate(ipiv(d)) 566 567 allocate(ipiv(d))
call cgetrf(d,d,wc_a,d,ipiv,info) 567 568 call cgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 568 569 if (info /= 0) then
status=0 569 570 status=0
deallocate(ipiv) 570 571 deallocate(ipiv)
deallocate(wc_a) 571 572 deallocate(wc_a)
return 572 573 return
end if 573 574 end if
allocate(work(2*d)) 574 575 allocate(work(2*d))
allocate(rwork(2*d)) 575 576 allocate(rwork(2*d))
call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 576 577 call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 577 578 if (info /= 0) then
status=0 578 579 status=0
end if 579 580 end if
deallocate(rwork) 580 581 deallocate(rwork)
deallocate(work) 581 582 deallocate(work)
deallocate(ipiv) 582 583 deallocate(ipiv)
deallocate(wc_a) 583 584 deallocate(wc_a)
end subroutine 584 585 end subroutine
585 586
subroutine fvn_z_matcon(d,a,rcond,status) 586 587 subroutine fvn_z_matcon(d,a,rcond,status)
! Matrix condition (reciprocal of condition number) 587 588 ! Matrix condition (reciprocal of condition number)
! 588 589 !
! d (in) : matrix rank 589 590 ! d (in) : matrix rank
! a (in) : The Matrix 590 591 ! a (in) : The Matrix
! rcond (out) : guess what 591 592 ! rcond (out) : guess what
! status (out) : =0 if something went wrong 592 593 ! status (out) : =0 if something went wrong
! 593 594 !
integer, intent(in) :: d 594 595 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 595 596 double complex, intent(in) :: a(d,d)
double precision, intent(out) :: rcond 596 597 double precision, intent(out) :: rcond
integer, intent(out) :: status 597 598 integer, intent(out) :: status
598 599
double complex, allocatable :: work(:) 599 600 double complex, allocatable :: work(:)
double precision, allocatable :: rwork(:) 600 601 double precision, allocatable :: rwork(:)
double precision :: anorm 601 602 double precision :: anorm
double complex, allocatable :: wc_a(:,:) ! working copy of a 602 603 double complex, allocatable :: wc_a(:,:) ! working copy of a
integer :: info 603 604 integer :: info
integer, allocatable :: ipiv(:) 604 605 integer, allocatable :: ipiv(:)
605 606
double precision, external :: zlange 606 607 double precision, external :: zlange
607 608
608 609
status=1 609 610 status=1
610 611
anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm 611 612 anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm
612 613
allocate(wc_a(d,d)) 613 614 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 614 615 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 615 616 wc_a(:,:)=a(:,:)
616 617
allocate(ipiv(d)) 617 618 allocate(ipiv(d))
call zgetrf(d,d,wc_a,d,ipiv,info) 618 619 call zgetrf(d,d,wc_a,d,ipiv,info)
if (info /= 0) then 619 620 if (info /= 0) then
status=0 620 621 status=0
deallocate(ipiv) 621 622 deallocate(ipiv)
deallocate(wc_a) 622 623 deallocate(wc_a)
return 623 624 return
end if 624 625 end if
625 626
allocate(work(2*d)) 626 627 allocate(work(2*d))
allocate(rwork(2*d)) 627 628 allocate(rwork(2*d))
call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) 628 629 call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info)
if (info /= 0) then 629 630 if (info /= 0) then
status=0 630 631 status=0
end if 631 632 end if
deallocate(rwork) 632 633 deallocate(rwork)
deallocate(work) 633 634 deallocate(work)
deallocate(ipiv) 634 635 deallocate(ipiv)
deallocate(wc_a) 635 636 deallocate(wc_a)
end subroutine 636 637 end subroutine
637 638
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 638 639 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 639 640 !
! Valeurs propres/ Vecteurs propre 640 641 ! Valeurs propres/ Vecteurs propre
! 641 642 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 642 643 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
643 644
subroutine fvn_s_matev(d,a,evala,eveca,status) 644 645 subroutine fvn_s_matev(d,a,evala,eveca,status)
! 645 646 !
! integer d (in) : matrice rank 646 647 ! integer d (in) : matrice rank
! real a(d,d) (in) : The Matrix 647 648 ! real a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 648 649 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 649 650 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 650 651 ! integer (out) : status =0 if something went wrong
! 651 652 !
! interfacing Lapack routine SGEEV 652 653 ! interfacing Lapack routine SGEEV
653 654
integer, intent(in) :: d 654 655 integer, intent(in) :: d
real, intent(in) :: a(d,d) 655 656 real, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 656 657 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 657 658 complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status 658 659 integer, intent(out) :: status
659 660
real, allocatable :: wc_a(:,:) ! a working copy of a 660 661 real, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 661 662 integer :: info
integer :: lwork 662 663 integer :: lwork
real, allocatable :: wr(:),wi(:) 663 664 real, allocatable :: wr(:),wi(:)
real :: vl ! unused but necessary for the call 664 665 real :: vl ! unused but necessary for the call
real, allocatable :: vr(:,:) 665 666 real, allocatable :: vr(:,:)
real, allocatable :: work(:) 666 667 real, allocatable :: work(:)
real :: twork(1) 667 668 real :: twork(1)
integer i 668 669 integer i
integer j 669 670 integer j
670 671
! making a working copy of a 671 672 ! making a working copy of a
allocate(wc_a(d,d)) 672 673 allocate(wc_a(d,d))
!call scopy(d*d,a,1,wc_a,1) 673 674 !call scopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 674 675 wc_a(:,:)=a(:,:)
675 676
allocate(wr(d)) 676 677 allocate(wr(d))
allocate(wi(d)) 677 678 allocate(wi(d))
allocate(vr(d,d)) 678 679 allocate(vr(d,d))
! query optimal work size 679 680 ! query optimal work size
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 680 681 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 681 682 lwork=int(twork(1))
allocate(work(lwork)) 682 683 allocate(work(lwork))
call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 683 684 call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
684 685
if (info /= 0) then 685 686 if (info /= 0) then
status=0 686 687 status=0
deallocate(work) 687 688 deallocate(work)
deallocate(vr) 688 689 deallocate(vr)
deallocate(wi) 689 690 deallocate(wi)
deallocate(wr) 690 691 deallocate(wr)
deallocate(wc_a) 691 692 deallocate(wc_a)
return 692 693 return
end if 693 694 end if
694 695
! now fill in the results 695 696 ! now fill in the results
i=1 696 697 i=1
do while(i<=d) 697 698 do while(i<=d)
evala(i)=cmplx(wr(i),wi(i)) 698 699 evala(i)=cmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 699 700 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=cmplx(vr(:,i),0.) 700 701 eveca(:,i)=cmplx(vr(:,i),0.)
else ! eigenvalue is complex 701 702 else ! eigenvalue is complex
evala(i+1)=cmplx(wr(i+1),wi(i+1)) 702 703 evala(i+1)=cmplx(wr(i+1),wi(i+1))
eveca(:,i)=cmplx(vr(:,i),vr(:,i+1)) 703 704 eveca(:,i)=cmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1)) 704 705 eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1))
i=i+1 705 706 i=i+1
end if 706 707 end if
i=i+1 707 708 i=i+1
enddo 708 709 enddo
deallocate(work) 709 710 deallocate(work)
deallocate(vr) 710 711 deallocate(vr)
deallocate(wi) 711 712 deallocate(wi)
deallocate(wr) 712 713 deallocate(wr)
deallocate(wc_a) 713 714 deallocate(wc_a)
714 715
end subroutine 715 716 end subroutine
716 717
subroutine fvn_d_matev(d,a,evala,eveca,status) 717 718 subroutine fvn_d_matev(d,a,evala,eveca,status)
! 718 719 !
! integer d (in) : matrice rank 719 720 ! integer d (in) : matrice rank
! double precision a(d,d) (in) : The Matrix 720 721 ! double precision a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 721 722 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 722 723 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 723 724 ! integer (out) : status =0 if something went wrong
! 724 725 !
! interfacing Lapack routine DGEEV 725 726 ! interfacing Lapack routine DGEEV
integer, intent(in) :: d 726 727 integer, intent(in) :: d
double precision, intent(in) :: a(d,d) 727 728 double precision, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 728 729 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 729 730 double complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status 730 731 integer, intent(out) :: status
731 732
double precision, allocatable :: wc_a(:,:) ! a working copy of a 732 733 double precision, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 733 734 integer :: info
integer :: lwork 734 735 integer :: lwork
double precision, allocatable :: wr(:),wi(:) 735 736 double precision, allocatable :: wr(:),wi(:)
double precision :: vl ! unused but necessary for the call 736 737 double precision :: vl ! unused but necessary for the call
double precision, allocatable :: vr(:,:) 737 738 double precision, allocatable :: vr(:,:)
double precision, allocatable :: work(:) 738 739 double precision, allocatable :: work(:)
double precision :: twork(1) 739 740 double precision :: twork(1)
integer i 740 741 integer i
integer j 741 742 integer j
742 743
! making a working copy of a 743 744 ! making a working copy of a
allocate(wc_a(d,d)) 744 745 allocate(wc_a(d,d))
!call dcopy(d*d,a,1,wc_a,1) 745 746 !call dcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 746 747 wc_a(:,:)=a(:,:)
747 748
allocate(wr(d)) 748 749 allocate(wr(d))
allocate(wi(d)) 749 750 allocate(wi(d))
allocate(vr(d,d)) 750 751 allocate(vr(d,d))
! query optimal work size 751 752 ! query optimal work size
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) 752 753 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info)
lwork=int(twork(1)) 753 754 lwork=int(twork(1))
allocate(work(lwork)) 754 755 allocate(work(lwork))
call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) 755 756 call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info)
756 757
if (info /= 0) then 757 758 if (info /= 0) then
status=0 758 759 status=0
deallocate(work) 759 760 deallocate(work)
deallocate(vr) 760 761 deallocate(vr)
deallocate(wi) 761 762 deallocate(wi)
deallocate(wr) 762 763 deallocate(wr)
deallocate(wc_a) 763 764 deallocate(wc_a)
return 764 765 return
end if 765 766 end if
766 767
! now fill in the results 767 768 ! now fill in the results
i=1 768 769 i=1
do while(i<=d) 769 770 do while(i<=d)
evala(i)=dcmplx(wr(i),wi(i)) 770 771 evala(i)=dcmplx(wr(i),wi(i))
if (wi(i) == 0.) then ! eigenvalue is real 771 772 if (wi(i) == 0.) then ! eigenvalue is real
eveca(:,i)=dcmplx(vr(:,i),0.) 772 773 eveca(:,i)=dcmplx(vr(:,i),0.)
else ! eigenvalue is complex 773 774 else ! eigenvalue is complex
evala(i+1)=dcmplx(wr(i+1),wi(i+1)) 774 775 evala(i+1)=dcmplx(wr(i+1),wi(i+1))
eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1)) 775 776 eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1))
eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1)) 776 777 eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1))
i=i+1 777 778 i=i+1
end if 778 779 end if
i=i+1 779 780 i=i+1
enddo 780 781 enddo
781 782
deallocate(work) 782 783 deallocate(work)
deallocate(vr) 783 784 deallocate(vr)
deallocate(wi) 784 785 deallocate(wi)
deallocate(wr) 785 786 deallocate(wr)
deallocate(wc_a) 786 787 deallocate(wc_a)
787 788
end subroutine 788 789 end subroutine
789 790
subroutine fvn_c_matev(d,a,evala,eveca,status) 790 791 subroutine fvn_c_matev(d,a,evala,eveca,status)
! 791 792 !
! integer d (in) : matrice rank 792 793 ! integer d (in) : matrice rank
! complex a(d,d) (in) : The Matrix 793 794 ! complex a(d,d) (in) : The Matrix
! complex evala(d) (out) : eigenvalues 794 795 ! complex evala(d) (out) : eigenvalues
! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 795 796 ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 796 797 ! integer (out) : status =0 if something went wrong
! 797 798 !
! interfacing Lapack routine CGEEV 798 799 ! interfacing Lapack routine CGEEV
799 800
integer, intent(in) :: d 800 801 integer, intent(in) :: d
complex, intent(in) :: a(d,d) 801 802 complex, intent(in) :: a(d,d)
complex, intent(out) :: evala(d) 802 803 complex, intent(out) :: evala(d)
complex, intent(out) :: eveca(d,d) 803 804 complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status 804 805 integer, intent(out) :: status
805 806
complex, allocatable :: wc_a(:,:) ! a working copy of a 806 807 complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 807 808 integer :: info
integer :: lwork 808 809 integer :: lwork
complex, allocatable :: work(:) 809 810 complex, allocatable :: work(:)
complex :: twork(1) 810 811 complex :: twork(1)
real, allocatable :: rwork(:) 811 812 real, allocatable :: rwork(:)
complex :: vl ! unused but necessary for the call 812 813 complex :: vl ! unused but necessary for the call
813 814
status=1 814 815 status=1
815 816
! making a working copy of a 816 817 ! making a working copy of a
allocate(wc_a(d,d)) 817 818 allocate(wc_a(d,d))
!call ccopy(d*d,a,1,wc_a,1) 818 819 !call ccopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 819 820 wc_a(:,:)=a(:,:)
820 821
821 822
! query optimal work size 822 823 ! query optimal work size
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 823 824 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 824 825 lwork=int(twork(1))
allocate(work(lwork)) 825 826 allocate(work(lwork))
allocate(rwork(2*d)) 826 827 allocate(rwork(2*d))
call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 827 828 call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
828 829
if (info /= 0) then 829 830 if (info /= 0) then
status=0 830 831 status=0
end if 831 832 end if
deallocate(rwork) 832 833 deallocate(rwork)
deallocate(work) 833 834 deallocate(work)
deallocate(wc_a) 834 835 deallocate(wc_a)
835 836
end subroutine 836 837 end subroutine
837 838
subroutine fvn_z_matev(d,a,evala,eveca,status) 838 839 subroutine fvn_z_matev(d,a,evala,eveca,status)
! 839 840 !
! integer d (in) : matrice rank 840 841 ! integer d (in) : matrice rank
! double complex a(d,d) (in) : The Matrix 841 842 ! double complex a(d,d) (in) : The Matrix
! double complex evala(d) (out) : eigenvalues 842 843 ! double complex evala(d) (out) : eigenvalues
! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector 843 844 ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector
! integer (out) : status =0 if something went wrong 844 845 ! integer (out) : status =0 if something went wrong
! 845 846 !
! interfacing Lapack routine ZGEEV 846 847 ! interfacing Lapack routine ZGEEV
847 848
integer, intent(in) :: d 848 849 integer, intent(in) :: d
double complex, intent(in) :: a(d,d) 849 850 double complex, intent(in) :: a(d,d)
double complex, intent(out) :: evala(d) 850 851 double complex, intent(out) :: evala(d)
double complex, intent(out) :: eveca(d,d) 851 852 double complex, intent(out) :: eveca(d,d)
integer, intent(out) :: status 852 853 integer, intent(out) :: status
853 854
double complex, allocatable :: wc_a(:,:) ! a working copy of a 854 855 double complex, allocatable :: wc_a(:,:) ! a working copy of a
integer :: info 855 856 integer :: info
integer :: lwork 856 857 integer :: lwork
double complex, allocatable :: work(:) 857 858 double complex, allocatable :: work(:)
double complex :: twork(1) 858 859 double complex :: twork(1)
double precision, allocatable :: rwork(:) 859 860 double precision, allocatable :: rwork(:)
double complex :: vl ! unused but necessary for the call 860 861 double complex :: vl ! unused but necessary for the call
861 862
status=1 862 863 status=1
863 864
! making a working copy of a 864 865 ! making a working copy of a
allocate(wc_a(d,d)) 865 866 allocate(wc_a(d,d))
!call zcopy(d*d,a,1,wc_a,1) 866 867 !call zcopy(d*d,a,1,wc_a,1)
wc_a(:,:)=a(:,:) 867 868 wc_a(:,:)=a(:,:)
868 869
! query optimal work size 869 870 ! query optimal work size
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) 870 871 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info)
lwork=int(twork(1)) 871 872 lwork=int(twork(1))
allocate(work(lwork)) 872 873 allocate(work(lwork))
allocate(rwork(2*d)) 873 874 allocate(rwork(2*d))
call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) 874 875 call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info)
875 876
if (info /= 0) then 876 877 if (info /= 0) then
status=0 877 878 status=0
end if 878 879 end if
deallocate(rwork) 879 880 deallocate(rwork)
deallocate(work) 880 881 deallocate(work)
deallocate(wc_a) 881 882 deallocate(wc_a)
882 883
end subroutine 883 884 end subroutine
884 885
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 885 886 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 886 887 !
! Akima spline interpolation and spline evaluation 887 888 ! Akima spline interpolation and spline evaluation
! 888 889 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 889 890 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
890 891
! Single precision 891 892 ! Single precision
subroutine fvn_s_akima(n,x,y,br,co) 892 893 subroutine fvn_s_akima(n,x,y,br,co)
implicit none 893 894 implicit none
integer, intent(in) :: n 894 895 integer, intent(in) :: n
real, intent(in) :: x(n) 895 896 real, intent(in) :: x(n)
real, intent(in) :: y(n) 896 897 real, intent(in) :: y(n)
real, intent(out) :: br(n) 897 898 real, intent(out) :: br(n)
real, intent(out) :: co(4,n) 898 899 real, intent(out) :: co(4,n)
899 900
real, allocatable :: var(:),z(:) 900 901 real, allocatable :: var(:),z(:)
real :: wi_1,wi 901 902 real :: wi_1,wi
integer :: i 902 903 integer :: i
real :: dx,a,b 903 904 real :: dx,a,b
904 905
! br is just a copy of x 905 906 ! br is just a copy of x
br(:)=x(:) 906 907 br(:)=x(:)
907 908
allocate(var(n)) 908 909 allocate(var(n))
allocate(z(n)) 909 910 allocate(z(n))
! evaluate the variations 910 911 ! evaluate the variations
do i=1, n-1 911 912 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 912 913 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 913 914 end do
var(n+2)=2.e0*var(n+1)-var(n) 914 915 var(n+2)=2.e0*var(n+1)-var(n)
var(n+3)=2.e0*var(n+2)-var(n+1) 915 916 var(n+3)=2.e0*var(n+2)-var(n+1)
var(2)=2.e0*var(3)-var(4) 916 917 var(2)=2.e0*var(3)-var(4)
var(1)=2.e0*var(2)-var(3) 917 918 var(1)=2.e0*var(2)-var(3)
918 919
do i = 1, n 919 920 do i = 1, n
wi_1=abs(var(i+3)-var(i+2)) 920 921 wi_1=abs(var(i+3)-var(i+2))
wi=abs(var(i+1)-var(i)) 921 922 wi=abs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.e0) then 922 923 if ((wi_1+wi).eq.0.e0) then
z(i)=(var(i+2)+var(i+1))/2.e0 923 924 z(i)=(var(i+2)+var(i+1))/2.e0
else 924 925 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 925 926 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 926 927 end if
end do 927 928 end do
928 929
do i=1, n-1 929 930 do i=1, n-1
dx=x(i+1)-x(i) 930 931 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 931 932 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 932 933 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 933 934 co(1,i)=y(i)
co(2,i)=z(i) 934 935 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd 935 936 !co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd
!co(4,i)=(a-2.*b)/dx**3 ! mรฉthode wd 936 937 !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 937 938 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 ! 938 939 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 939 940 ! 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 ... 940 941 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 941 942 end do
co(1,n)=y(n) 942 943 co(1,n)=y(n)
co(2,n)=z(n) 943 944 co(2,n)=z(n)
co(3,n)=0.e0 944 945 co(3,n)=0.e0
co(4,n)=0.e0 945 946 co(4,n)=0.e0
946 947
deallocate(z) 947 948 deallocate(z)
deallocate(var) 948 949 deallocate(var)
949 950
end subroutine 950 951 end subroutine
951 952
! Double precision 952 953 ! Double precision
subroutine fvn_d_akima(n,x,y,br,co) 953 954 subroutine fvn_d_akima(n,x,y,br,co)
954 955
implicit none 955 956 implicit none
integer, intent(in) :: n 956 957 integer, intent(in) :: n
double precision, intent(in) :: x(n) 957 958 double precision, intent(in) :: x(n)
double precision, intent(in) :: y(n) 958 959 double precision, intent(in) :: y(n)
double precision, intent(out) :: br(n) 959 960 double precision, intent(out) :: br(n)
double precision, intent(out) :: co(4,n) 960 961 double precision, intent(out) :: co(4,n)
961 962
double precision, allocatable :: var(:),z(:) 962 963 double precision, allocatable :: var(:),z(:)
double precision :: wi_1,wi 963 964 double precision :: wi_1,wi
integer :: i 964 965 integer :: i
double precision :: dx,a,b 965 966 double precision :: dx,a,b
966 967
! br is just a copy of x 967 968 ! br is just a copy of x
br(:)=x(:) 968 969 br(:)=x(:)
969 970
allocate(var(n)) 970 971 allocate(var(n))
allocate(z(n)) 971 972 allocate(z(n))
! evaluate the variations 972 973 ! evaluate the variations
do i=1, n-1 973 974 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 974 975 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 975 976 end do
var(n+2)=2.d0*var(n+1)-var(n) 976 977 var(n+2)=2.d0*var(n+1)-var(n)
var(n+3)=2.d0*var(n+2)-var(n+1) 977 978 var(n+3)=2.d0*var(n+2)-var(n+1)
var(2)=2.d0*var(3)-var(4) 978 979 var(2)=2.d0*var(3)-var(4)
var(1)=2.d0*var(2)-var(3) 979 980 var(1)=2.d0*var(2)-var(3)
980 981
do i = 1, n 981 982 do i = 1, n
wi_1=dabs(var(i+3)-var(i+2)) 982 983 wi_1=dabs(var(i+3)-var(i+2))
wi=dabs(var(i+1)-var(i)) 983 984 wi=dabs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.d0) then 984 985 if ((wi_1+wi).eq.0.d0) then
z(i)=(var(i+2)+var(i+1))/2.d0 985 986 z(i)=(var(i+2)+var(i+1))/2.d0
else 986 987 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 987 988 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 988 989 end if
end do 989 990 end do
990 991
do i=1, n-1 991 992 do i=1, n-1
dx=x(i+1)-x(i) 992 993 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 993 994 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 994 995 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 995 996 co(1,i)=y(i)
co(2,i)=z(i) 996 997 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd 997 998 !co(3,i)=-(a-3.*b)/dx**2 ! mรฉthode wd
!co(4,i)=(a-2.*b)/dx**3 ! mรฉthode wd 998 999 !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 999 1000 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 ! 1000 1001 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 1001 1002 ! 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 ... 1002 1003 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 1003 1004 end do
co(1,n)=y(n) 1004 1005 co(1,n)=y(n)
co(2,n)=z(n) 1005 1006 co(2,n)=z(n)
co(3,n)=0.d0 1006 1007 co(3,n)=0.d0
co(4,n)=0.d0 1007 1008 co(4,n)=0.d0
1008 1009
deallocate(z) 1009 1010 deallocate(z)
deallocate(var) 1010 1011 deallocate(var)
1011 1012
end subroutine 1012 1013 end subroutine
1013 1014
! 1014 1015 !
! Single precision spline evaluation 1015 1016 ! Single precision spline evaluation
! 1016 1017 !
function fvn_s_spline_eval(x,n,br,co) 1017 1018 function fvn_s_spline_eval(x,n,br,co)
implicit none 1018 1019 implicit none
real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated 1019 1020 real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated
integer, intent(in) :: n ! number of intervals 1020 1021 integer, intent(in) :: n ! number of intervals
real, intent(in) :: br(n+1) ! breakpoints 1021 1022 real, intent(in) :: br(n+1) ! breakpoints
real, intent(in) :: co(4,n+1) ! spline coeeficients 1022 1023 real, intent(in) :: co(4,n+1) ! spline coeeficients
real :: fvn_s_spline_eval 1023 1024 real :: fvn_s_spline_eval
1024 1025
integer :: i 1025 1026 integer :: i
real :: dx 1026 1027 real :: dx
1027 1028
if (x<=br(1)) then 1028 1029 if (x<=br(1)) then
i=1 1029 1030 i=1
else if (x>=br(n+1)) then 1030 1031 else if (x>=br(n+1)) then
i=n 1031 1032 i=n
else 1032 1033 else
i=1 1033 1034 i=1
do while(x>=br(i)) 1034 1035 do while(x>=br(i))
i=i+1 1035 1036 i=i+1
end do 1036 1037 end do
i=i-1 1037 1038 i=i-1
end if 1038 1039 end if
dx=x-br(i) 1039 1040 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 1040 1041 fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
1041 1042
end function 1042 1043 end function
1043 1044
! Double precision spline evaluation 1044 1045 ! Double precision spline evaluation
function fvn_d_spline_eval(x,n,br,co) 1045 1046 function fvn_d_spline_eval(x,n,br,co)
implicit none 1046 1047 implicit none
double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated 1047 1048 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 1048 1049 integer, intent(in) :: n ! number of intervals
double precision, intent(in) :: br(n+1) ! breakpoints 1049 1050 double precision, intent(in) :: br(n+1) ! breakpoints
double precision, intent(in) :: co(4,n+1) ! spline coeeficients 1050 1051 double precision, intent(in) :: co(4,n+1) ! spline coeeficients
double precision :: fvn_d_spline_eval 1051 1052 double precision :: fvn_d_spline_eval
1052 1053
integer :: i 1053 1054 integer :: i
double precision :: dx 1054 1055 double precision :: dx
1055 1056
1056 1057
if (x<=br(1)) then 1057 1058 if (x<=br(1)) then
i=1 1058 1059 i=1
else if (x>=br(n+1)) then 1059 1060 else if (x>=br(n+1)) then
i=n 1060 1061 i=n
else 1061 1062 else
i=1 1062 1063 i=1
do while(x>=br(i)) 1063 1064 do while(x>=br(i))
i=i+1 1064 1065 i=i+1
end do 1065 1066 end do
i=i-1 1066 1067 i=i-1
end if 1067 1068 end if
1068 1069
dx=x-br(i) 1069 1070 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 1070 1071 fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
1071 1072
end function 1072 1073 end function
1073 1074
1074 1075
! 1075 1076 !
! Muller 1076 1077 ! Muller
! 1077 1078 !
! 1078 1079 !
! 1079 1080 !
! William Daniau 2007 1080 1081 ! William Daniau 2007
! 1081 1082 !
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f 1082 1083 ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f 1083 1084 ! http://plato.asu.edu/ftp/other_software/muller.f
! 1084 1085 !
! it can be used as a replacement for imsl routine dzanly with minor changes 1085 1086 ! it can be used as a replacement for imsl routine dzanly with minor changes
! 1086 1087 !
!----------------------------------------------------------------------- 1087 1088 !-----------------------------------------------------------------------
! 1088 1089 !
! purpose - zeros of an analytic complex function 1089 1090 ! purpose - zeros of an analytic complex function
! using the muller method with deflation 1090 1091 ! using the muller method with deflation
! 1091 1092 !
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, 1092 1093 ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier) 1093 1094 ! infer,ier)
! 1094 1095 !
! arguments f - a complex function subprogram, f(z), written 1095 1096 ! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose 1096 1097 ! by the user specifying the equation whose
! roots are to be found. f must appear in 1097 1098 ! roots are to be found. f must appear in
! an external statement in the calling pro- 1098 1099 ! an external statement in the calling pro-
! gram. 1099 1100 ! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p 1100 1101 ! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) 1101 1102 ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found 1102 1103 ! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and. 1103 1104 ! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted 1104 1105 ! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input) 1105 1106 ! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted 1106 1107 ! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given 1107 1108 ! if two successive approximations to a given
! root agree within eps1. (input) 1108 1109 ! root agree within eps1. (input)
! note. if either or both of the stopping 1109 1110 ! note. if either or both of the stopping
! criteria are fulfilled, the root is 1110 1111 ! criteria are fulfilled, the root is
! accepted. 1111 1112 ! accepted.
! kn - the number of known roots which must be stored 1112 1113 ! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller 1113 1114 ! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these 1114 1115 ! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),..., 1115 1116 ! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal 1116 1117 ! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input) 1117 1118 ! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by 1118 1119 ! n - the number of new roots to be found by
! muller (input) 1119 1120 ! muller (input)
! x - a complex vector of length kn+n. x(1),..., 1120 1121 ! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known 1121 1122 ! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may, 1122 1123 ! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for 1123 1124 ! on user option, contain initial guesses for
! the n new roots which are to be computed. 1124 1125 ! the n new roots which are to be computed.
! if the user does not provide an initial 1125 1126 ! if the user does not provide an initial
! guess, zero is used. 1126 1127 ! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the 1127 1128 ! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller. 1128 1129 ! approximate roots found by muller.
! itmax - the maximum allowable number of iterations 1129 1130 ! itmax - the maximum allowable number of iterations
! per root (input) 1130 1131 ! per root (input)
! infer - an integer vector of length kn+n. on 1131 1132 ! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of 1132 1133 ! output infer(j) contains the number of
! iterations used in finding the j-th root 1133 1134 ! iterations used in finding the j-th root
! when convergence was achieved. if 1134 1135 ! when convergence was achieved. if
! convergence was not obtained in itmax 1135 1136 ! convergence was not obtained in itmax
! iterations, infer(j) will be greater than 1136 1137 ! iterations, infer(j) will be greater than
! itmax (output). 1137 1138 ! itmax (output).
! ier - error parameter (output) 1138 1139 ! ier - error parameter (output)
! warning error 1139 1140 ! warning error
! ier = 33 indicates failure to converge with- 1140 1141 ! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of 1141 1142 ! in itmax iterations for at least one of
! the (n) new roots. 1142 1143 ! the (n) new roots.
! 1143 1144 !
! 1144 1145 !
! remarks muller always returns the last approximation for root j 1145 1146 ! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied, 1146 1147 ! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the 1147 1148 ! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j) 1148 1149 ! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater 1149 1150 ! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did 1150 1151 ! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter- 1151 1152 ! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax 1152 1153 ! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con- 1153 1154 ! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla- 1154 1155 ! vergence was obtained (on iteration k) for the defla-
! ted function 1155 1156 ! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) 1156 1157 ! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
! 1157 1158 !
! but failed for f(z). in this case, better initial 1158 1159 ! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax 1159 1160 ! guesses might help or, it might be necessary to relax
! the convergence criterion. 1160 1161 ! the convergence criterion.
! 1161 1162 !
!----------------------------------------------------------------------- 1162 1163 !-----------------------------------------------------------------------
! 1163 1164 !
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) 1164 1165 subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none 1165 1166 implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq 1166 1167 double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, & 1167 1168 double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & 1168 1169 tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5 1169 1170 zero,p1,one,four,p5
1170 1171
double complex, external :: f 1171 1172 double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & 1172 1173 integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode 1173 1174 knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n) 1174 1175 double complex :: x(kn+n)
integer :: infer(kn+n) 1175 1176 integer :: infer(kn+n)
1176 1177
1177 1178
data zero/(0.0,0.0)/,p1/(0.1,0.0)/, & 1178 1179 data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
one/(1.0,0.0)/,four/(4.0,0.0)/, & 1179 1180 one/(1.0,0.0)/,four/(4.0,0.0)/, &
p5/(0.5,0.0)/, & 1180 1181 p5/(0.5,0.0)/, &
rzero/0.0/,rten/10.0/,rhun/100.0/, & 1181 1182 rzero/0.0/,rten/10.0/,rhun/100.0/, &
ax/0.1/,ickmax/3/,rp01/0.01/ 1182 1183 ax/0.1/,ickmax/3/,rp01/0.01/
1183 1184
ier = 0 1184 1185 ier = 0
if (n .lt. 1) then ! What the hell are doing here then ... 1185 1186 if (n .lt. 1) then ! What the hell are doing here then ...
return 1186 1187 return
end if 1187 1188 end if
!eps1 = rten **(-nsig) 1188 1189 !eps1 = rten **(-nsig)
eps1 = min(eps1,rp01) 1189 1190 eps1 = min(eps1,rp01)
1190 1191
knp1 = kn+1 1191 1192 knp1 = kn+1
knpn = kn+n 1192 1193 knpn = kn+n
knpng = kn+nguess 1193 1194 knpng = kn+nguess
do i=1,knpn 1194 1195 do i=1,knpn
infer(i) = 0 1195 1196 infer(i) = 0
if (i .gt. knpng) x(i) = zero 1196 1197 if (i .gt. knpng) x(i) = zero
end do 1197 1198 end do
l= knp1 1198 1199 l= knp1
1199 1200
ic=0 1200 1201 ic=0
rloop: do while (l<=knpn) ! Main loop over new roots 1201 1202 rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0 1202 1203 jk = 0
ick = 0 1203 1204 ick = 0
xl = x(l) 1204 1205 xl = x(l)
icloop: do 1205 1206 icloop: do
ic = 0 1206 1207 ic = 0
h = ax 1207 1208 h = ax
h = p1*h 1208 1209 h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl 1209 1210 if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are 1210 1211 ! first three points are
! xl+h, xl-h, xl 1211 1212 ! xl+h, xl-h, xl
rt = xl+h 1212 1213 rt = xl+h
call deflated_work(errcode) 1213 1214 call deflated_work(errcode)
if (errcode == 1) then 1214 1215 if (errcode == 1) then
exit icloop 1215 1216 exit icloop
end if 1216 1217 end if
1217 1218
z0 = fprt 1218 1219 z0 = fprt
y0 = frt 1219 1220 y0 = frt
x0 = rt 1220 1221 x0 = rt
rt = xl-h 1221 1222 rt = xl-h
call deflated_work(errcode) 1222 1223 call deflated_work(errcode)
if (errcode == 1) then 1223 1224 if (errcode == 1) then
exit icloop 1224 1225 exit icloop
end if 1225 1226 end if
1226 1227
z1 = fprt 1227 1228 z1 = fprt
y1 = frt 1228 1229 y1 = frt
h = xl-rt 1229 1230 h = xl-rt
d = h/(rt-x0) 1230 1231 d = h/(rt-x0)
rt = xl 1231 1232 rt = xl
1232 1233
call deflated_work(errcode) 1233 1234 call deflated_work(errcode)
if (errcode == 1) then 1234 1235 if (errcode == 1) then
exit icloop 1235 1236 exit icloop
end if 1236 1237 end if
1237 1238
1238 1239
z2 = fprt 1239 1240 z2 = fprt
y2 = frt 1240 1241 y2 = frt
! begin main algorithm 1241 1242 ! begin main algorithm
iloop: do 1242 1243 iloop: do
dd = one + d 1243 1244 dd = one + d
t1 = z0*d*d 1244 1245 t1 = z0*d*d
t2 = z1*dd*dd 1245 1246 t2 = z1*dd*dd
xx = z2*dd 1246 1247 xx = z2*dd
t3 = z2*d 1247 1248 t3 = z2*d
bi = t1-t2+xx+t3 1248 1249 bi = t1-t2+xx+t3
den = bi*bi-four*(xx*t1-t3*(t2-xx)) 1249 1250 den = bi*bi-four*(xx*t1-t3*(t2-xx))
! use denominator of maximum amplitude 1250 1251 ! use denominator of maximum amplitude
t1 = cdsqrt(den) 1251 1252 t1 = cdsqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1)) 1252 1253 qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1 1253 1254 t2 = bi + t1
tpq = cdabs(t2)+qz 1254 1255 tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero 1255 1256 if (tpq .eq. qz) t2 = zero
t3 = bi - t1 1256 1257 t3 = bi - t1
tpq = cdabs(t3) + qz 1257 1258 tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero 1258 1259 if (tpq .eq. qz) t3 = zero
den = t2 1259 1260 den = t2
qz = cdabs(t3)-cdabs(t2) 1260 1261 qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3 1261 1262 if (qz .gt. rzero) den = t3
! test for zero denominator 1262 1263 ! test for zero denominator
if (cdabs(den) .eq. rzero) then 1263 1264 if (cdabs(den) .eq. rzero) then
call trans_rt() 1264 1265 call trans_rt()
call deflated_work(errcode) 1265 1266 call deflated_work(errcode)
if (errcode == 1) then 1266 1267 if (errcode == 1) then
exit icloop 1267 1268 exit icloop
end if 1268 1269 end if
z2 = fprt 1269 1270 z2 = fprt
y2 = frt 1270 1271 y2 = frt
cycle iloop 1271 1272 cycle iloop
end if 1272 1273 end if
1273 1274
1274 1275
d = -xx/den 1275 1276 d = -xx/den
d = d+d 1276 1277 d = d+d
h = d*h 1277 1278 h = d*h
rt = rt + h 1278 1279 rt = rt + h
! check convergence of the first kind 1279 1280 ! check convergence of the first kind
if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then 1280 1281 if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
if (ic .ne. 0) then 1281 1282 if (ic .ne. 0) then
exit icloop 1282 1283 exit icloop
end if 1283 1284 end if
ic = 1 1284 1285 ic = 1
z0 = y1 1285 1286 z0 = y1
z1 = y2 1286 1287 z1 = y2
z2 = f(rt) 1287 1288 z2 = f(rt)
xl = rt 1288 1289 xl = rt
ick = ick+1 1289 1290 ick = ick+1
if (ick .le. ickmax) then 1290 1291 if (ick .le. ickmax) then
cycle iloop 1291 1292 cycle iloop
end if 1292 1293 end if
! warning error, itmax = maximum 1293 1294 ! warning error, itmax = maximum
jk = itmax + jk 1294 1295 jk = itmax + jk
ier = 33 1295 1296 ier = 33
end if 1296 1297 end if
if (ic .ne. 0) then 1297 1298 if (ic .ne. 0) then
cycle icloop 1298 1299 cycle icloop
end if 1299 1300 end if
call deflated_work(errcode) 1300 1301 call deflated_work(errcode)
if (errcode == 1) then 1301 1302 if (errcode == 1) then
exit icloop 1302 1303 exit icloop
end if 1303 1304 end if
1304 1305
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) 1305 1306 do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce 1306 1307 ! take remedial action to induce
! convergence 1307 1308 ! convergence
d = d*p5 1308 1309 d = d*p5
h = h*p5 1309 1310 h = h*p5
rt = rt-h 1310 1311 rt = rt-h
call deflated_work(errcode) 1311 1312 call deflated_work(errcode)
if (errcode == 1) then 1312 1313 if (errcode == 1) then
exit icloop 1313 1314 exit icloop
end if 1314 1315 end if
end do 1315 1316 end do
z0 = z1 1316 1317 z0 = z1
z1 = z2 1317 1318 z1 = z2
z2 = fprt 1318 1319 z2 = fprt
y0 = y1 1319 1320 y0 = y1
y1 = y2 1320 1321 y1 = y2
y2 = frt 1321 1322 y2 = frt
end do iloop 1322 1323 end do iloop
end do icloop 1323 1324 end do icloop
x(l) = rt 1324 1325 x(l) = rt
infer(l) = jk 1325 1326 infer(l) = jk
l = l+1 1326 1327 l = l+1
end do rloop 1327 1328 end do rloop
1328 1329
contains 1329 1330 contains
subroutine trans_rt() 1330 1331 subroutine trans_rt()
tem = rten*eps1 1331 1332 tem = rten*eps1
if (cdabs(rt) .gt. ax) tem = tem*rt 1332 1333 if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem 1333 1334 rt = rt+tem
d = (h+tem)*d/h 1334 1335 d = (h+tem)*d/h
h = h+tem 1335 1336 h = h+tem
end subroutine trans_rt 1336 1337 end subroutine trans_rt
1337 1338
subroutine deflated_work(errcode) 1338 1339 subroutine deflated_work(errcode)
! errcode=0 => no errors 1339 1340 ! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved 1340 1341 ! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag 1341 1342 integer :: errcode,flag
1342 1343
flag=1 1343 1344 flag=1
loop1: do while(flag==1) 1344 1345 loop1: do while(flag==1)
errcode=0 1345 1346 errcode=0
jk = jk+1 1346 1347 jk = jk+1
if (jk .gt. itmax) then 1347 1348 if (jk .gt. itmax) then
ier=33 1348 1349 ier=33
errcode=1 1349 1350 errcode=1
return 1350 1351 return
end if 1351 1352 end if
frt = f(rt) 1352 1353 frt = f(rt)
fprt = frt 1353 1354 fprt = frt
if (l /= 1) then 1354 1355 if (l /= 1) then
lm1 = l-1 1355 1356 lm1 = l-1
do i=1,lm1 1356 1357 do i=1,lm1
tem = rt - x(i) 1357 1358 tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then 1358 1359 if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible? 1359 1360 !if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt() 1360 1361 call trans_rt()
cycle loop1 1361 1362 cycle loop1
end if 1362 1363 end if
fprt = fprt/tem 1363 1364 fprt = fprt/tem
end do 1364 1365 end do
end if 1365 1366 end if
flag=0 1366 1367 flag=0
end do loop1 1367 1368 end do loop1
1368 1369
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then 1369 1370 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1 1370 1371 errcode=1
return 1371 1372 return
end if 1372 1373 end if
1373 1374
end subroutine deflated_work 1374 1375 end subroutine deflated_work
1375 1376
end subroutine 1376 1377 end subroutine
1377 1378
1378 1379
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1379 1380 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1380 1381 !
! Integration 1381 1382 ! Integration
! 1382 1383 !
! Only double precision coded atm 1383 1384 ! Only double precision coded atm
! 1384 1385 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1385 1386 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1386 1387
1387 1388
subroutine fvn_d_gauss_legendre(n,qx,qw) 1388 1389 subroutine fvn_d_gauss_legendre(n,qx,qw)
! 1389 1390 !
! This routine compute the n Gauss Legendre abscissas and weights 1390 1391 ! This routine compute the n Gauss Legendre abscissas and weights
! Adapted from Numerical Recipes routine gauleg 1391 1392 ! Adapted from Numerical Recipes routine gauleg
! 1392 1393 !
! n (in) : number of points 1393 1394 ! n (in) : number of points
! qx(out) : abscissas 1394 1395 ! qx(out) : abscissas
! qw(out) : weights 1395 1396 ! qw(out) : weights
! 1396 1397 !
implicit none 1397 1398 implicit none
double precision,parameter :: pi=3.141592653589793d0 1398 1399 double precision,parameter :: pi=3.141592653589793d0
integer, intent(in) :: n 1399 1400 integer, intent(in) :: n
double precision, intent(out) :: qx(n),qw(n) 1400 1401 double precision, intent(out) :: qx(n),qw(n)
1401 1402
integer :: m,i,j 1402 1403 integer :: m,i,j
double precision :: z,z1,p1,p2,p3,pp 1403 1404 double precision :: z,z1,p1,p2,p3,pp
m=(n+1)/2 1404 1405 m=(n+1)/2
do i=1,m 1405 1406 do i=1,m
z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0)) 1406 1407 z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
iloop: do 1407 1408 iloop: do
p1=1.d0 1408 1409 p1=1.d0
p2=0.d0 1409 1410 p2=0.d0
do j=1,n 1410 1411 do j=1,n
p3=p2 1411 1412 p3=p2
p2=p1 1412 1413 p2=p1
p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j) 1413 1414 p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
end do 1414 1415 end do
pp=dble(n)*(z*p1-p2)/(z*z-1.d0) 1415 1416 pp=dble(n)*(z*p1-p2)/(z*z-1.d0)
z1=z 1416 1417 z1=z
z=z1-p1/pp 1417 1418 z=z1-p1/pp
if (dabs(z-z1)<=epsilon(z)) then 1418 1419 if (dabs(z-z1)<=epsilon(z)) then
exit iloop 1419 1420 exit iloop
end if 1420 1421 end if
end do iloop 1421 1422 end do iloop
qx(i)=-z 1422 1423 qx(i)=-z
qx(n+1-i)=z 1423 1424 qx(n+1-i)=z
qw(i)=2.d0/((1.d0-z*z)*pp*pp) 1424 1425 qw(i)=2.d0/((1.d0-z*z)*pp*pp)
qw(n+1-i)=qw(i) 1425 1426 qw(n+1-i)=qw(i)
end do 1426 1427 end do
end subroutine 1427 1428 end subroutine
1428 1429
1429 1430
1430 1431
subroutine fvn_d_gl_integ(f,a,b,n,res) 1431 1432 subroutine fvn_d_gl_integ(f,a,b,n,res)
! 1432 1433 !
! This is a simple non adaptative integration routine 1433 1434 ! This is a simple non adaptative integration routine
! using n gauss legendre abscissas and weights 1434 1435 ! using n gauss legendre abscissas and weights
! 1435 1436 !
! f(in) : the function to integrate 1436 1437 ! f(in) : the function to integrate
! a(in) : lower bound 1437 1438 ! a(in) : lower bound
! b(in) : higher bound 1438 1439 ! b(in) : higher bound
! n(in) : number of gauss legendre pairs 1439 1440 ! n(in) : number of gauss legendre pairs
! res(out): the evaluation of the integral 1440 1441 ! res(out): the evaluation of the integral
! 1441 1442 !
double precision,external :: f 1442 1443 double precision,external :: f
double precision, intent(in) :: a,b 1443 1444 double precision, intent(in) :: a,b
integer, intent(in):: n 1444 1445 integer, intent(in):: n
double precision, intent(out) :: res 1445 1446 double precision, intent(out) :: res
1446 1447
double precision, allocatable :: qx(:),qw(:) 1447 1448 double precision, allocatable :: qx(:),qw(:)
double precision :: xm,xr 1448 1449 double precision :: xm,xr
integer :: i 1449 1450 integer :: i
1450 1451
! First compute n gauss legendre abs and weight 1451 1452 ! First compute n gauss legendre abs and weight
allocate(qx(n)) 1452 1453 allocate(qx(n))
allocate(qw(n)) 1453 1454 allocate(qw(n))
call fvn_d_gauss_legendre(n,qx,qw) 1454 1455 call fvn_d_gauss_legendre(n,qx,qw)
1455 1456
xm=0.5d0*(b+a) 1456 1457 xm=0.5d0*(b+a)
xr=0.5d0*(b-a) 1457 1458 xr=0.5d0*(b-a)
1458 1459
res=0.d0 1459 1460 res=0.d0
1460 1461
do i=1,n 1461 1462 do i=1,n
res=res+qw(i)*f(xm+xr*qx(i)) 1462 1463 res=res+qw(i)*f(xm+xr*qx(i))
end do 1463 1464 end do
1464 1465
res=xr*res 1465 1466 res=xr*res
1466 1467
deallocate(qw) 1467 1468 deallocate(qw)
deallocate(qx) 1468 1469 deallocate(qx)
1469 1470
end subroutine 1470 1471 end subroutine
1471 1472
!!!!!!!!!!!!!!!!!!!!!!!! 1472 1473 !!!!!!!!!!!!!!!!!!!!!!!!
! 1473 1474 !
! Simple and double adaptative Gauss Kronrod integration based on 1474 1475 ! Simple and double adaptative Gauss Kronrod integration based on
! a modified version of quadpack ( http://www.netlib.org/quadpack 1475 1476 ! a modified version of quadpack ( http://www.netlib.org/quadpack
! 1476 1477 !
! Common parameters : 1477 1478 ! Common parameters :
! 1478 1479 !
! key (in) 1479 1480 ! key (in)
! epsabs 1480 1481 ! epsabs
! epsrel 1481 1482 ! epsrel
! 1482 1483 !
! 1483 1484 !
!!!!!!!!!!!!!!!!!!!!!!!! 1484 1485 !!!!!!!!!!!!!!!!!!!!!!!!
1485 1486
subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit) 1486 1487 subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
! 1487 1488 !
! Evaluate the integral of function f(x) between a and b 1488 1489 ! Evaluate the integral of function f(x) between a and b
! 1489 1490 !
! f(in) : the function 1490 1491 ! f(in) : the function
! a(in) : lower bound 1491 1492 ! a(in) : lower bound
! b(in) : higher bound 1492 1493 ! b(in) : higher bound
! epsabs(in) : desired absolute error 1493 1494 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 1494 1495 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 1495 1496 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 1496 1497 ! 1: 7 - 15 points
! 2: 10 - 21 points 1497 1498 ! 2: 10 - 21 points
! 3: 15 - 31 points 1498 1499 ! 3: 15 - 31 points
! 4: 20 - 41 points 1499 1500 ! 4: 20 - 41 points
! 5: 25 - 51 points 1500 1501 ! 5: 25 - 51 points
! 6: 30 - 61 points 1501 1502 ! 6: 30 - 61 points
! 1502 1503 !
! limit(in) : maximum number of subintervals in the partition of the 1503 1504 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 1504 1505 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 1505 1506 ! behaviour as the imsl routine dqdag
! 1506 1507 !
! res(out) : estimated integral value 1507 1508 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 1508 1509 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 1509 1510 ! ier(out) : error flag from quadpack routines
! 0 : no error 1510 1511 ! 0 : no error
! 1 : maximum number of subdivisions allowed 1511 1512 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 1512 1513 ! has been achieved. one can allow more
! subdivisions by increasing the value of 1513 1514 ! subdivisions by increasing the value of
! limit (and taking the according dimension 1514 1515 ! limit (and taking the according dimension
! adjustments into account). however, if 1515 1516 ! adjustments into account). however, if
! this yield no improvement it is advised 1516 1517 ! this yield no improvement it is advised
! to analyze the integrand in order to 1517 1518 ! to analyze the integrand in order to
! determine the integration difficulaties. 1518 1519 ! determine the integration difficulaties.
! if the position of a local difficulty can 1519 1520 ! if the position of a local difficulty can
! be determined (i.e.singularity, 1520 1521 ! be determined (i.e.singularity,
! discontinuity within the interval) one 1521 1522 ! discontinuity within the interval) one
! will probably gain from splitting up the 1522 1523 ! will probably gain from splitting up the
! interval at this point and calling the 1523 1524 ! interval at this point and calling the
! integrator on the subranges. if possible, 1524 1525 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 1525 1526 ! an appropriate special-purpose integrator
! should be used which is designed for 1526 1527 ! should be used which is designed for
! handling the type of difficulty involved. 1527 1528 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 1528 1529 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 1529 1530 ! detected, which prevents the requested
! tolerance from being achieved. 1530 1531 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 1531 1532 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 1532 1533 ! at some points of the integration
! interval. 1533 1534 ! interval.
! 6 : the input is invalid, because 1534 1535 ! 6 : the input is invalid, because
! (epsabs.le.0 and 1535 1536 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 1536 1537 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 1537 1538 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 1538 1539 ! result, abserr, neval, last are set
! to zero. 1539 1540 ! to zero.
! except when lenw is invalid, iwork(1), 1540 1541 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 1541 1542 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 1542 1543 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 1543 1544 ! work(limit+1) to b.
1544 1545
implicit none 1545 1546 implicit none
double precision, external :: f 1546 1547 double precision, external :: f
double precision, intent(in) :: a,b,epsabs,epsrel 1547 1548 double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key 1548 1549 integer, intent(in) :: key
integer, intent(in) :: limit 1549 1550 integer, intent(in) :: limit
double precision, intent(out) :: res,abserr 1550 1551 double precision, intent(out) :: res,abserr
integer, intent(out) :: ier 1551 1552 integer, intent(out) :: ier
1552 1553
double precision, allocatable :: work(:) 1553 1554 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 1554 1555 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 1555 1556 integer :: lenw,neval,last
1556 1557
! imsl value for limit is 500 1557 1558 ! imsl value for limit is 500
lenw=limit*4 1558 1559 lenw=limit*4
1559 1560
allocate(iwork(limit)) 1560 1561 allocate(iwork(limit))
allocate(work(lenw)) 1561 1562 allocate(work(lenw))
1562 1563
call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 1563 1564 call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
1564 1565
deallocate(work) 1565 1566 deallocate(work)
deallocate(iwork) 1566 1567 deallocate(iwork)
1567 1568
end subroutine 1568 1569 end subroutine
1569 1570
1570 1571
1571 1572