lxsd / fvn

Download as
Browse Code »

Commit 15adbf52baf8ed530882453ee63cbe9f4d1c2132

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

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

Showing 2 changed files with 539 additions and 0 deletions Inline Diff

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