lxsd / fvn

Download as
Browse Code ยป

Commit ef002e0aa8946d1a2c09d106b627a3093d1033d8

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

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

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