lxsd / fvn

Download as
Browse Code »

Commit 126a3aed07763d1b0975159db844a6348597f38e

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

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

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

1 1
module fvn 2 2 module fvn
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 3 3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 4 4 !
! fvn : a f95 module replacement for some imsl routines 5 5 ! fvn : a f95 module replacement for some imsl routines
! it uses lapack for linear algebra 6 6 ! it uses lapack for linear algebra
! it uses modified quadpack for integration 7 7 ! it uses modified quadpack for integration
! 8 8 !
! William Daniau 2007 9 9 ! William Daniau 2007
! william.daniau@femto-st.fr 10 10 ! william.daniau@femto-st.fr
! 11 11 !
! Routines naming scheme : 12 12 ! Routines naming scheme :
! 13 13 !
! fvn_x_name 14 14 ! fvn_x_name
! where x can be s : real 15 15 ! where x can be s : real
! d : real double precision 16 16 ! d : real double precision
! c : complex 17 17 ! c : complex
! z : double complex 18 18 ! z : double complex
! 19 19 !
! 20 20 !
! This piece of code is totally free! Do whatever you want with it. However 21 21 ! This piece of code is totally free! Do whatever you want with it. However
! if you find it usefull it would be kind to give credits ;-) 22 22 ! if you find it usefull it would be kind to give credits ;-)
! 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 905 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 906 906 !
! Akima spline interpolation and spline evaluation 907 907 ! Akima spline interpolation and spline evaluation
! 908 908 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 909 909 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
910 910
! Single precision 911 911 ! Single precision
subroutine fvn_s_akima(n,x,y,br,co) 912 912 subroutine fvn_s_akima(n,x,y,br,co)
implicit none 913 913 implicit none
integer, intent(in) :: n 914 914 integer, intent(in) :: n
real, intent(in) :: x(n) 915 915 real, intent(in) :: x(n)
real, intent(in) :: y(n) 916 916 real, intent(in) :: y(n)
real, intent(out) :: br(n) 917 917 real, intent(out) :: br(n)
real, intent(out) :: co(4,n) 918 918 real, intent(out) :: co(4,n)
919 919
real, allocatable :: var(:),z(:) 920 920 real, allocatable :: var(:),z(:)
real :: wi_1,wi 921 921 real :: wi_1,wi
integer :: i 922 922 integer :: i
real :: dx,a,b 923 923 real :: dx,a,b
924 924
! br is just a copy of x 925 925 ! br is just a copy of x
br(:)=x(:) 926 926 br(:)=x(:)
927 927
allocate(var(n)) 928 928 allocate(var(n))
allocate(z(n)) 929 929 allocate(z(n))
! evaluate the variations 930 930 ! evaluate the variations
do i=1, n-1 931 931 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 932 932 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 933 933 end do
var(n+2)=2.e0*var(n+1)-var(n) 934 934 var(n+2)=2.e0*var(n+1)-var(n)
var(n+3)=2.e0*var(n+2)-var(n+1) 935 935 var(n+3)=2.e0*var(n+2)-var(n+1)
var(2)=2.e0*var(3)-var(4) 936 936 var(2)=2.e0*var(3)-var(4)
var(1)=2.e0*var(2)-var(3) 937 937 var(1)=2.e0*var(2)-var(3)
938 938
do i = 1, n 939 939 do i = 1, n
wi_1=abs(var(i+3)-var(i+2)) 940 940 wi_1=abs(var(i+3)-var(i+2))
wi=abs(var(i+1)-var(i)) 941 941 wi=abs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.e0) then 942 942 if ((wi_1+wi).eq.0.e0) then
z(i)=(var(i+2)+var(i+1))/2.e0 943 943 z(i)=(var(i+2)+var(i+1))/2.e0
else 944 944 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 945 945 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 946 946 end if
end do 947 947 end do
948 948
do i=1, n-1 949 949 do i=1, n-1
dx=x(i+1)-x(i) 950 950 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 951 951 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 952 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 953 953 co(1,i)=y(i)
co(2,i)=z(i) 954 954 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd 955 955 !co(3,i)=-(a-3.*b)/dx**2 ! méthode wd
!co(4,i)=(a-2.*b)/dx**3 ! méthode wd 956 956 !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 957 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 958 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 959 ! 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 960 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 961 961 end do
co(1,n)=y(n) 962 962 co(1,n)=y(n)
co(2,n)=z(n) 963 963 co(2,n)=z(n)
co(3,n)=0.e0 964 964 co(3,n)=0.e0
co(4,n)=0.e0 965 965 co(4,n)=0.e0
966 966
deallocate(z) 967 967 deallocate(z)
deallocate(var) 968 968 deallocate(var)
969 969
end subroutine 970 970 end subroutine
971 971
! Double precision 972 972 ! Double precision
subroutine fvn_d_akima(n,x,y,br,co) 973 973 subroutine fvn_d_akima(n,x,y,br,co)
974 974
implicit none 975 975 implicit none
integer, intent(in) :: n 976 976 integer, intent(in) :: n
double precision, intent(in) :: x(n) 977 977 double precision, intent(in) :: x(n)
double precision, intent(in) :: y(n) 978 978 double precision, intent(in) :: y(n)
double precision, intent(out) :: br(n) 979 979 double precision, intent(out) :: br(n)
double precision, intent(out) :: co(4,n) 980 980 double precision, intent(out) :: co(4,n)
981 981
double precision, allocatable :: var(:),z(:) 982 982 double precision, allocatable :: var(:),z(:)
double precision :: wi_1,wi 983 983 double precision :: wi_1,wi
integer :: i 984 984 integer :: i
double precision :: dx,a,b 985 985 double precision :: dx,a,b
986 986
! br is just a copy of x 987 987 ! br is just a copy of x
br(:)=x(:) 988 988 br(:)=x(:)
989 989
allocate(var(n)) 990 990 allocate(var(n))
allocate(z(n)) 991 991 allocate(z(n))
! evaluate the variations 992 992 ! evaluate the variations
do i=1, n-1 993 993 do i=1, n-1
var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) 994 994 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i))
end do 995 995 end do
var(n+2)=2.d0*var(n+1)-var(n) 996 996 var(n+2)=2.d0*var(n+1)-var(n)
var(n+3)=2.d0*var(n+2)-var(n+1) 997 997 var(n+3)=2.d0*var(n+2)-var(n+1)
var(2)=2.d0*var(3)-var(4) 998 998 var(2)=2.d0*var(3)-var(4)
var(1)=2.d0*var(2)-var(3) 999 999 var(1)=2.d0*var(2)-var(3)
1000 1000
do i = 1, n 1001 1001 do i = 1, n
wi_1=dabs(var(i+3)-var(i+2)) 1002 1002 wi_1=dabs(var(i+3)-var(i+2))
wi=dabs(var(i+1)-var(i)) 1003 1003 wi=dabs(var(i+1)-var(i))
if ((wi_1+wi).eq.0.d0) then 1004 1004 if ((wi_1+wi).eq.0.d0) then
z(i)=(var(i+2)+var(i+1))/2.d0 1005 1005 z(i)=(var(i+2)+var(i+1))/2.d0
else 1006 1006 else
z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) 1007 1007 z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi)
end if 1008 1008 end if
end do 1009 1009 end do
1010 1010
do i=1, n-1 1011 1011 do i=1, n-1
dx=x(i+1)-x(i) 1012 1012 dx=x(i+1)-x(i)
a=(z(i+1)-z(i))*dx ! coeff intermediaires pour calcul wd 1013 1013 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 1014 b=y(i+1)-y(i)-z(i)*dx ! coeff intermediaires pour calcul wd
co(1,i)=y(i) 1015 1015 co(1,i)=y(i)
co(2,i)=z(i) 1016 1016 co(2,i)=z(i)
!co(3,i)=-(a-3.*b)/dx**2 ! méthode wd 1017 1017 !co(3,i)=-(a-3.*b)/dx**2 ! méthode wd
!co(4,i)=(a-2.*b)/dx**3 ! méthode wd 1018 1018 !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 1019 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 1020 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 1021 ! 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 1022 ! etrangement la fonction csval corrige et donne la bonne valeur ...
end do 1023 1023 end do
co(1,n)=y(n) 1024 1024 co(1,n)=y(n)
co(2,n)=z(n) 1025 1025 co(2,n)=z(n)
co(3,n)=0.d0 1026 1026 co(3,n)=0.d0
co(4,n)=0.d0 1027 1027 co(4,n)=0.d0
1028 1028
deallocate(z) 1029 1029 deallocate(z)
deallocate(var) 1030 1030 deallocate(var)
1031 1031
end subroutine 1032 1032 end subroutine
1033 1033
! 1034 1034 !
! Single precision spline evaluation 1035 1035 ! Single precision spline evaluation
! 1036 1036 !
function fvn_s_spline_eval(x,n,br,co) 1037 1037 function fvn_s_spline_eval(x,n,br,co)
implicit none 1038 1038 implicit none
real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated 1039 1039 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 1040 integer, intent(in) :: n ! number of intervals
real, intent(in) :: br(n+1) ! breakpoints 1041 1041 real, intent(in) :: br(n+1) ! breakpoints
real, intent(in) :: co(4,n+1) ! spline coeeficients 1042 1042 real, intent(in) :: co(4,n+1) ! spline coeeficients
real :: fvn_s_spline_eval 1043 1043 real :: fvn_s_spline_eval
1044 1044
integer :: i 1045 1045 integer :: i
real :: dx 1046 1046 real :: dx
1047 1047
if (x<=br(1)) then 1048 1048 if (x<=br(1)) then
i=1 1049 1049 i=1
else if (x>=br(n+1)) then 1050 1050 else if (x>=br(n+1)) then
i=n 1051 1051 i=n
else 1052 1052 else
i=1 1053 1053 i=1
do while(x>=br(i)) 1054 1054 do while(x>=br(i))
i=i+1 1055 1055 i=i+1
end do 1056 1056 end do
i=i-1 1057 1057 i=i-1
end if 1058 1058 end if
dx=x-br(i) 1059 1059 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 1060 fvn_s_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
1061 1061
end function 1062 1062 end function
1063 1063
! Double precision spline evaluation 1064 1064 ! Double precision spline evaluation
function fvn_d_spline_eval(x,n,br,co) 1065 1065 function fvn_d_spline_eval(x,n,br,co)
implicit none 1066 1066 implicit none
double precision, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated 1067 1067 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 1068 integer, intent(in) :: n ! number of intervals
double precision, intent(in) :: br(n+1) ! breakpoints 1069 1069 double precision, intent(in) :: br(n+1) ! breakpoints
double precision, intent(in) :: co(4,n+1) ! spline coeeficients 1070 1070 double precision, intent(in) :: co(4,n+1) ! spline coeeficients
double precision :: fvn_d_spline_eval 1071 1071 double precision :: fvn_d_spline_eval
1072 1072
integer :: i 1073 1073 integer :: i
double precision :: dx 1074 1074 double precision :: dx
1075 1075
1076 1076
if (x<=br(1)) then 1077 1077 if (x<=br(1)) then
i=1 1078 1078 i=1
else if (x>=br(n+1)) then 1079 1079 else if (x>=br(n+1)) then
i=n 1080 1080 i=n
else 1081 1081 else
i=1 1082 1082 i=1
do while(x>=br(i)) 1083 1083 do while(x>=br(i))
i=i+1 1084 1084 i=i+1
end do 1085 1085 end do
i=i-1 1086 1086 i=i-1
end if 1087 1087 end if
1088 1088
dx=x-br(i) 1089 1089 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 1090 fvn_d_spline_eval=co(1,i)+co(2,i)*dx+co(3,i)*dx**2+co(4,i)*dx**3
1091 1091
end function 1092 1092 end function
1093 1093
1094 1094
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1095 1095 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1096 1096 !
! Least square problem 1097 1097 ! Least square problem
! 1098 1098 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1099 1099 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1100 1100 !
! 1101 1101 !
1102
1103
1104
1105
subroutine fvn_d_lspoly(np,x,y,deg,coeff,status) 1102 1106 subroutine fvn_d_lspoly(np,x,y,deg,coeff,status)
! 1103 1107 !
! Least square polynomial fitting 1104 1108 ! Least square polynomial fitting
! 1105 1109 !
! Find the coefficients of the least square polynomial of a given degree 1106 1110 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1107 1111 ! for a set of coordinates.
! 1108 1112 !
! The degree must be lower than the number of points 1109 1113 ! The degree must be lower than the number of points
! 1110 1114 !
! np (in) : number of points 1111 1115 ! np (in) : number of points
! x(np) (in) : x data 1112 1116 ! x(np) (in) : x data
! y(np) (in) : y data 1113 1117 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1114 1118 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1115 1119 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1116 1120 ! status (out) : =0 if a problem occurs
! 1117 1121 !
implicit none 1118 1122 implicit none
1119 1123
integer, intent(in) :: np,deg 1120 1124 integer, intent(in) :: np,deg
real(kind=8), intent(in), dimension(np) :: x,y 1121 1125 real(kind=8), intent(in), dimension(np) :: x,y
real(kind=8), intent(out), dimension(deg+1) :: coeff 1122 1126 real(kind=8), intent(out), dimension(deg+1) :: coeff
integer, intent(out) :: status 1123 1127 integer, intent(out) :: status
1124 1128
real(kind=8), allocatable, dimension(:,:) :: mat,bmat 1125 1129 real(kind=8), allocatable, dimension(:,:) :: mat,bmat
1130 real(kind=8),dimension(:),allocatable :: work
1131 real(kind=8),dimension(1) :: twork
1132 integer :: lwork,info
1133
1134 integer :: i,j
1135
1136 status=1
1137 allocate(mat(np,deg+1),bmat(np,1))
1138
1139 ! Design matrix valorisation
1140 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1141
1142 ! second member valorisation
1143 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1144
1145 ! query workspace size
1146 call dgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
1147 lwork=twork(1)
1148 allocate(work(int(lwork)))
1149 ! real call
1150 call dgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
1151
1152 if (info /= 0) then
1153 status=0
1154 end if
1155
1156 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1157
1158 deallocate(work)
1159 deallocate(mat,bmat)
1160 end subroutine
1161
1162 subroutine fvn_s_lspoly(np,x,y,deg,coeff,status)
1163 !
1164 ! Least square polynomial fitting
1165 !
1166 ! Find the coefficients of the least square polynomial of a given degree
1167 ! for a set of coordinates.
1168 !
1169 ! The degree must be lower than the number of points
1170 !
1171 ! np (in) : number of points
1172 ! x(np) (in) : x data
1173 ! y(np) (in) : y data
1174 ! deg (in) : polynomial's degree
1175 ! coeff(deg+1) (out) : polynomial coefficients
1176 ! status (out) : =0 if a problem occurs
1177 !
1178 implicit none
1179
1180 integer, intent(in) :: np,deg
1181 real(kind=4), intent(in), dimension(np) :: x,y
1182 real(kind=4), intent(out), dimension(deg+1) :: coeff
1183 integer, intent(out) :: status
1184
1185 real(kind=4), allocatable, dimension(:,:) :: mat,bmat
1186 real(kind=4),dimension(:),allocatable :: work
1187 real(kind=4),dimension(1) :: twork
1188 integer :: lwork,info
1189
1190 integer :: i,j
1191
1192 status=1
1193 allocate(mat(np,deg+1),bmat(np,1))
1194
1195 ! Design matrix valorisation
1196 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1197
1198 ! second member valorisation
1199 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1200
1201 ! query workspace size
1202 call sgels('N',np,deg+1,1,mat,np,bmat,np,twork,-1,info)
1203 lwork=twork(1)
1204 allocate(work(int(lwork)))
1205 ! real call
1206 call sgels('N',np,deg+1,1,mat,np,bmat,np,work,lwork,info)
1207
1208 if (info /= 0) then
1209 status=0
1210 end if
1211
1212 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1213
1214 deallocate(work)
1215 deallocate(mat,bmat)
1216 end subroutine
1217
1218
1219
1220
1221
1222
1223
1224
1225 subroutine fvn_d_lspoly_svd(np,x,y,deg,coeff,status)
1226 !
1227 ! Least square polynomial fitting
1228 !
1229 ! Find the coefficients of the least square polynomial of a given degree
1230 ! for a set of coordinates.
1231 !
1232 ! The degree must be lower than the number of points
1233 !
1234 ! np (in) : number of points
1235 ! x(np) (in) : x data
1236 ! y(np) (in) : y data
1237 ! deg (in) : polynomial's degree
1238 ! coeff(deg+1) (out) : polynomial coefficients
1239 ! status (out) : =0 if a problem occurs
1240 !
1241 implicit none
1242
1243 integer, intent(in) :: np,deg
1244 real(kind=8), intent(in), dimension(np) :: x,y
1245 real(kind=8), intent(out), dimension(deg+1) :: coeff
1246 integer, intent(out) :: status
1247
1248 real(kind=8), allocatable, dimension(:,:) :: mat,bmat
real(kind=8),dimension(:),allocatable :: work,singval 1126 1249 real(kind=8),dimension(:),allocatable :: work,singval
real(kind=8),dimension(1) :: twork 1127 1250 real(kind=8),dimension(1) :: twork
integer :: lwork,info,rank 1128 1251 integer :: lwork,info,rank
1129 1252
integer :: i,j 1130 1253 integer :: i,j
1131 1254
status=1 1132 1255 status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1)) 1133 1256 allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
1134 1257
! Design matrix valorisation 1135 1258 ! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) ) 1136 1259 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1137 1260
! second member valorisation 1138 1261 ! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /)) 1139 1262 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1140 1263
! query workspace size 1141 1264 ! query workspace size
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info) 1142 1265 call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1) 1143 1266 lwork=twork(1)
allocate(work(int(lwork))) 1144 1267 allocate(work(int(lwork)))
! real call 1145 1268 ! real call
call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info) 1146 1269 call dgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
1147 1270
if (info /= 0) then 1148 1271 if (info /= 0) then
status=0 1149 1272 status=0
end if 1150 1273 end if
1151 1274
coeff = (/ (bmat(i,1),i=1,deg+1) /) 1152 1275 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1153 1276
deallocate(work) 1154 1277 deallocate(work)
deallocate(mat,bmat,singval) 1155 1278 deallocate(mat,bmat,singval)
end subroutine 1156 1279 end subroutine
1157 1280
subroutine fvn_s_lspoly(np,x,y,deg,coeff,status) 1158 1281 subroutine fvn_s_lspoly_svd(np,x,y,deg,coeff,status)
! 1159 1282 !
! Least square polynomial fitting 1160 1283 ! Least square polynomial fitting
! 1161 1284 !
! Find the coefficients of the least square polynomial of a given degree 1162 1285 ! Find the coefficients of the least square polynomial of a given degree
! for a set of coordinates. 1163 1286 ! for a set of coordinates.
! 1164 1287 !
! The degree must be lower than the number of points 1165 1288 ! The degree must be lower than the number of points
! 1166 1289 !
! np (in) : number of points 1167 1290 ! np (in) : number of points
! x(np) (in) : x data 1168 1291 ! x(np) (in) : x data
! y(np) (in) : y data 1169 1292 ! y(np) (in) : y data
! deg (in) : polynomial's degree 1170 1293 ! deg (in) : polynomial's degree
! coeff(deg+1) (out) : polynomial coefficients 1171 1294 ! coeff(deg+1) (out) : polynomial coefficients
! status (out) : =0 if a problem occurs 1172 1295 ! status (out) : =0 if a problem occurs
! 1173 1296 !
implicit none 1174 1297 implicit none
1175 1298
integer, intent(in) :: np,deg 1176 1299 integer, intent(in) :: np,deg
real(kind=4), intent(in), dimension(np) :: x,y 1177 1300 real(kind=4), intent(in), dimension(np) :: x,y
real(kind=4), intent(out), dimension(deg+1) :: coeff 1178 1301 real(kind=4), intent(out), dimension(deg+1) :: coeff
integer, intent(out) :: status 1179 1302 integer, intent(out) :: status
1180 1303
real(kind=4), allocatable, dimension(:,:) :: mat,bmat 1181 1304 real(kind=4), allocatable, dimension(:,:) :: mat,bmat
real(kind=4),dimension(:),allocatable :: work,singval 1182 1305 real(kind=4),dimension(:),allocatable :: work,singval
real(kind=4),dimension(1) :: twork 1183 1306 real(kind=4),dimension(1) :: twork
integer :: lwork,info,rank 1184 1307 integer :: lwork,info,rank
1185 1308
integer :: i,j 1186 1309 integer :: i,j
1187 1310
status=1 1188 1311 status=1
allocate(mat(np,deg+1),bmat(np,1),singval(deg+1)) 1189 1312 allocate(mat(np,deg+1),bmat(np,1),singval(deg+1))
1190 1313
! Design matrix valorisation 1191 1314 ! Design matrix valorisation
mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) ) 1192 1315 mat=reshape( (/ ((x(i)**(j-1),i=1,np),j=1,deg+1) /),shape=(/ np,deg+1 /) )
1193 1316
! second member valorisation 1194 1317 ! second member valorisation
bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /)) 1195 1318 bmat=reshape ( (/ (y(i),i=1,np) /) ,shape = (/ np,1 /))
1196 1319
! query workspace size 1197 1320 ! query workspace size
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info) 1198 1321 call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,twork,-1,info)
lwork=twork(1) 1199 1322 lwork=twork(1)
allocate(work(int(lwork))) 1200 1323 allocate(work(int(lwork)))
! real call 1201 1324 ! real call
call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info) 1202 1325 call sgelss(np,deg+1,1,mat,np,bmat,np,singval,-1.,rank,work,lwork,info)
1203 1326
if (info /= 0) then 1204 1327 if (info /= 0) then
status=0 1205 1328 status=0
end if 1206 1329 end if
1207 1330
coeff = (/ (bmat(i,1),i=1,deg+1) /) 1208 1331 coeff = (/ (bmat(i,1),i=1,deg+1) /)
1209 1332
deallocate(work) 1210 1333 deallocate(work)
deallocate(mat,bmat,singval) 1211 1334 deallocate(mat,bmat,singval)
end subroutine 1212 1335 end subroutine
1213 1336
1214 1337
! 1215 1338 !
! Muller 1216 1339 ! Muller
! 1217 1340 !
! 1218 1341 !
! 1219 1342 !
! William Daniau 2007 1220 1343 ! William Daniau 2007
! 1221 1344 !
! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f 1222 1345 ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f
! http://plato.asu.edu/ftp/other_software/muller.f 1223 1346 ! http://plato.asu.edu/ftp/other_software/muller.f
! 1224 1347 !
! it can be used as a replacement for imsl routine dzanly with minor changes 1225 1348 ! it can be used as a replacement for imsl routine dzanly with minor changes
! 1226 1349 !
!----------------------------------------------------------------------- 1227 1350 !-----------------------------------------------------------------------
! 1228 1351 !
! purpose - zeros of an analytic complex function 1229 1352 ! purpose - zeros of an analytic complex function
! using the muller method with deflation 1230 1353 ! using the muller method with deflation
! 1231 1354 !
! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, 1232 1355 ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax,
! infer,ier) 1233 1356 ! infer,ier)
! 1234 1357 !
! arguments f - a complex function subprogram, f(z), written 1235 1358 ! arguments f - a complex function subprogram, f(z), written
! by the user specifying the equation whose 1236 1359 ! by the user specifying the equation whose
! roots are to be found. f must appear in 1237 1360 ! roots are to be found. f must appear in
! an external statement in the calling pro- 1238 1361 ! an external statement in the calling pro-
! gram. 1239 1362 ! gram.
! eps - 1st stopping criterion. let fp(z)=f(z)/p 1240 1363 ! eps - 1st stopping criterion. let fp(z)=f(z)/p
! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) 1241 1364 ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1))
! and z(1),...,z(k-1) are previously found 1242 1365 ! and z(1),...,z(k-1) are previously found
! roots. if ((cdabs(f(z)).le.eps) .and. 1243 1366 ! roots. if ((cdabs(f(z)).le.eps) .and.
! (cdabs(fp(z)).le.eps)), then z is accepted 1244 1367 ! (cdabs(fp(z)).le.eps)), then z is accepted
! as a root. (input) 1245 1368 ! as a root. (input)
! eps1 - 2nd stopping criterion. a root is accepted 1246 1369 ! eps1 - 2nd stopping criterion. a root is accepted
! if two successive approximations to a given 1247 1370 ! if two successive approximations to a given
! root agree within eps1. (input) 1248 1371 ! root agree within eps1. (input)
! note. if either or both of the stopping 1249 1372 ! note. if either or both of the stopping
! criteria are fulfilled, the root is 1250 1373 ! criteria are fulfilled, the root is
! accepted. 1251 1374 ! accepted.
! kn - the number of known roots which must be stored 1252 1375 ! kn - the number of known roots which must be stored
! in x(1),...,x(kn), prior to entry to muller 1253 1376 ! in x(1),...,x(kn), prior to entry to muller
! nguess - the number of initial guesses provided. these 1254 1377 ! nguess - the number of initial guesses provided. these
! guesses must be stored in x(kn+1),..., 1255 1378 ! guesses must be stored in x(kn+1),...,
! x(kn+nguess). nguess must be set equal 1256 1379 ! x(kn+nguess). nguess must be set equal
! to zero if no guesses are provided. (input) 1257 1380 ! to zero if no guesses are provided. (input)
! n - the number of new roots to be found by 1258 1381 ! n - the number of new roots to be found by
! muller (input) 1259 1382 ! muller (input)
! x - a complex vector of length kn+n. x(1),..., 1260 1383 ! x - a complex vector of length kn+n. x(1),...,
! x(kn) on input must contain any known 1261 1384 ! x(kn) on input must contain any known
! roots. x(kn+1),..., x(kn+n) on input may, 1262 1385 ! roots. x(kn+1),..., x(kn+n) on input may,
! on user option, contain initial guesses for 1263 1386 ! on user option, contain initial guesses for
! the n new roots which are to be computed. 1264 1387 ! the n new roots which are to be computed.
! if the user does not provide an initial 1265 1388 ! if the user does not provide an initial
! guess, zero is used. 1266 1389 ! guess, zero is used.
! on output, x(kn+1),...,x(kn+n) contain the 1267 1390 ! on output, x(kn+1),...,x(kn+n) contain the
! approximate roots found by muller. 1268 1391 ! approximate roots found by muller.
! itmax - the maximum allowable number of iterations 1269 1392 ! itmax - the maximum allowable number of iterations
! per root (input) 1270 1393 ! per root (input)
! infer - an integer vector of length kn+n. on 1271 1394 ! infer - an integer vector of length kn+n. on
! output infer(j) contains the number of 1272 1395 ! output infer(j) contains the number of
! iterations used in finding the j-th root 1273 1396 ! iterations used in finding the j-th root
! when convergence was achieved. if 1274 1397 ! when convergence was achieved. if
! convergence was not obtained in itmax 1275 1398 ! convergence was not obtained in itmax
! iterations, infer(j) will be greater than 1276 1399 ! iterations, infer(j) will be greater than
! itmax (output). 1277 1400 ! itmax (output).
! ier - error parameter (output) 1278 1401 ! ier - error parameter (output)
! warning error 1279 1402 ! warning error
! ier = 33 indicates failure to converge with- 1280 1403 ! ier = 33 indicates failure to converge with-
! in itmax iterations for at least one of 1281 1404 ! in itmax iterations for at least one of
! the (n) new roots. 1282 1405 ! the (n) new roots.
! 1283 1406 !
! 1284 1407 !
! remarks muller always returns the last approximation for root j 1285 1408 ! remarks muller always returns the last approximation for root j
! in x(j). if the convergence criterion is satisfied, 1286 1409 ! in x(j). if the convergence criterion is satisfied,
! then infer(j) is less than or equal to itmax. if the 1287 1410 ! then infer(j) is less than or equal to itmax. if the
! convergence criterion is not satisified, then infer(j) 1288 1411 ! convergence criterion is not satisified, then infer(j)
! is set to either itmax+1 or itmax+k, with k greater 1289 1412 ! is set to either itmax+1 or itmax+k, with k greater
! than 1. infer(j) = itmax+1 indicates that muller did 1290 1413 ! than 1. infer(j) = itmax+1 indicates that muller did
! not obtain convergence in the allowed number of iter- 1291 1414 ! not obtain convergence in the allowed number of iter-
! ations. in this case, the user may wish to set itmax 1292 1415 ! ations. in this case, the user may wish to set itmax
! to a larger value. infer(j) = itmax+k means that con- 1293 1416 ! to a larger value. infer(j) = itmax+k means that con-
! vergence was obtained (on iteration k) for the defla- 1294 1417 ! vergence was obtained (on iteration k) for the defla-
! ted function 1295 1418 ! ted function
! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) 1296 1419 ! fp(z) = f(z)/((z-z(1)...(z-z(j-1)))
! 1297 1420 !
! but failed for f(z). in this case, better initial 1298 1421 ! but failed for f(z). in this case, better initial
! guesses might help or, it might be necessary to relax 1299 1422 ! guesses might help or, it might be necessary to relax
! the convergence criterion. 1300 1423 ! the convergence criterion.
! 1301 1424 !
!----------------------------------------------------------------------- 1302 1425 !-----------------------------------------------------------------------
! 1303 1426 !
subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) 1304 1427 subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier)
implicit none 1305 1428 implicit none
double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq 1306 1429 double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq
double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, & 1307 1430 double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, &
tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & 1308 1431 tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, &
zero,p1,one,four,p5 1309 1432 zero,p1,one,four,p5
1310 1433
double complex, external :: f 1311 1434 double complex, external :: f
integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & 1312 1435 integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, &
knpng,jk,ick,nn,lm1,errcode 1313 1436 knpng,jk,ick,nn,lm1,errcode
double complex :: x(kn+n) 1314 1437 double complex :: x(kn+n)
integer :: infer(kn+n) 1315 1438 integer :: infer(kn+n)
1316 1439
1317 1440
data zero/(0.0,0.0)/,p1/(0.1,0.0)/, & 1318 1441 data zero/(0.0,0.0)/,p1/(0.1,0.0)/, &
one/(1.0,0.0)/,four/(4.0,0.0)/, & 1319 1442 one/(1.0,0.0)/,four/(4.0,0.0)/, &
p5/(0.5,0.0)/, & 1320 1443 p5/(0.5,0.0)/, &
rzero/0.0/,rten/10.0/,rhun/100.0/, & 1321 1444 rzero/0.0/,rten/10.0/,rhun/100.0/, &
ax/0.1/,ickmax/3/,rp01/0.01/ 1322 1445 ax/0.1/,ickmax/3/,rp01/0.01/
1323 1446
ier = 0 1324 1447 ier = 0
if (n .lt. 1) then ! What the hell are doing here then ... 1325 1448 if (n .lt. 1) then ! What the hell are doing here then ...
return 1326 1449 return
end if 1327 1450 end if
!eps1 = rten **(-nsig) 1328 1451 !eps1 = rten **(-nsig)
eps1 = min(eps1,rp01) 1329 1452 eps1 = min(eps1,rp01)
1330 1453
knp1 = kn+1 1331 1454 knp1 = kn+1
knpn = kn+n 1332 1455 knpn = kn+n
knpng = kn+nguess 1333 1456 knpng = kn+nguess
do i=1,knpn 1334 1457 do i=1,knpn
infer(i) = 0 1335 1458 infer(i) = 0
if (i .gt. knpng) x(i) = zero 1336 1459 if (i .gt. knpng) x(i) = zero
end do 1337 1460 end do
l= knp1 1338 1461 l= knp1
1339 1462
ic=0 1340 1463 ic=0
rloop: do while (l<=knpn) ! Main loop over new roots 1341 1464 rloop: do while (l<=knpn) ! Main loop over new roots
jk = 0 1342 1465 jk = 0
ick = 0 1343 1466 ick = 0
xl = x(l) 1344 1467 xl = x(l)
icloop: do 1345 1468 icloop: do
ic = 0 1346 1469 ic = 0
h = ax 1347 1470 h = ax
h = p1*h 1348 1471 h = p1*h
if (cdabs(xl) .gt. ax) h = p1*xl 1349 1472 if (cdabs(xl) .gt. ax) h = p1*xl
! first three points are 1350 1473 ! first three points are
! xl+h, xl-h, xl 1351 1474 ! xl+h, xl-h, xl
rt = xl+h 1352 1475 rt = xl+h
call deflated_work(errcode) 1353 1476 call deflated_work(errcode)
if (errcode == 1) then 1354 1477 if (errcode == 1) then
exit icloop 1355 1478 exit icloop
end if 1356 1479 end if
1357 1480
z0 = fprt 1358 1481 z0 = fprt
y0 = frt 1359 1482 y0 = frt
x0 = rt 1360 1483 x0 = rt
rt = xl-h 1361 1484 rt = xl-h
call deflated_work(errcode) 1362 1485 call deflated_work(errcode)
if (errcode == 1) then 1363 1486 if (errcode == 1) then
exit icloop 1364 1487 exit icloop
end if 1365 1488 end if
1366 1489
z1 = fprt 1367 1490 z1 = fprt
y1 = frt 1368 1491 y1 = frt
h = xl-rt 1369 1492 h = xl-rt
d = h/(rt-x0) 1370 1493 d = h/(rt-x0)
rt = xl 1371 1494 rt = xl
1372 1495
call deflated_work(errcode) 1373 1496 call deflated_work(errcode)
if (errcode == 1) then 1374 1497 if (errcode == 1) then
exit icloop 1375 1498 exit icloop
end if 1376 1499 end if
1377 1500
1378 1501
z2 = fprt 1379 1502 z2 = fprt
y2 = frt 1380 1503 y2 = frt
! begin main algorithm 1381 1504 ! begin main algorithm
iloop: do 1382 1505 iloop: do
dd = one + d 1383 1506 dd = one + d
t1 = z0*d*d 1384 1507 t1 = z0*d*d
t2 = z1*dd*dd 1385 1508 t2 = z1*dd*dd
xx = z2*dd 1386 1509 xx = z2*dd
t3 = z2*d 1387 1510 t3 = z2*d
bi = t1-t2+xx+t3 1388 1511 bi = t1-t2+xx+t3
den = bi*bi-four*(xx*t1-t3*(t2-xx)) 1389 1512 den = bi*bi-four*(xx*t1-t3*(t2-xx))
! use denominator of maximum amplitude 1390 1513 ! use denominator of maximum amplitude
t1 = cdsqrt(den) 1391 1514 t1 = cdsqrt(den)
qz = rhun*max(cdabs(bi),cdabs(t1)) 1392 1515 qz = rhun*max(cdabs(bi),cdabs(t1))
t2 = bi + t1 1393 1516 t2 = bi + t1
tpq = cdabs(t2)+qz 1394 1517 tpq = cdabs(t2)+qz
if (tpq .eq. qz) t2 = zero 1395 1518 if (tpq .eq. qz) t2 = zero
t3 = bi - t1 1396 1519 t3 = bi - t1
tpq = cdabs(t3) + qz 1397 1520 tpq = cdabs(t3) + qz
if (tpq .eq. qz) t3 = zero 1398 1521 if (tpq .eq. qz) t3 = zero
den = t2 1399 1522 den = t2
qz = cdabs(t3)-cdabs(t2) 1400 1523 qz = cdabs(t3)-cdabs(t2)
if (qz .gt. rzero) den = t3 1401 1524 if (qz .gt. rzero) den = t3
! test for zero denominator 1402 1525 ! test for zero denominator
if (cdabs(den) .eq. rzero) then 1403 1526 if (cdabs(den) .eq. rzero) then
call trans_rt() 1404 1527 call trans_rt()
call deflated_work(errcode) 1405 1528 call deflated_work(errcode)
if (errcode == 1) then 1406 1529 if (errcode == 1) then
exit icloop 1407 1530 exit icloop
end if 1408 1531 end if
z2 = fprt 1409 1532 z2 = fprt
y2 = frt 1410 1533 y2 = frt
cycle iloop 1411 1534 cycle iloop
end if 1412 1535 end if
1413 1536
1414 1537
d = -xx/den 1415 1538 d = -xx/den
d = d+d 1416 1539 d = d+d
h = d*h 1417 1540 h = d*h
rt = rt + h 1418 1541 rt = rt + h
! check convergence of the first kind 1419 1542 ! check convergence of the first kind
if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then 1420 1543 if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then
if (ic .ne. 0) then 1421 1544 if (ic .ne. 0) then
exit icloop 1422 1545 exit icloop
end if 1423 1546 end if
ic = 1 1424 1547 ic = 1
z0 = y1 1425 1548 z0 = y1
z1 = y2 1426 1549 z1 = y2
z2 = f(rt) 1427 1550 z2 = f(rt)
xl = rt 1428 1551 xl = rt
ick = ick+1 1429 1552 ick = ick+1
if (ick .le. ickmax) then 1430 1553 if (ick .le. ickmax) then
cycle iloop 1431 1554 cycle iloop
end if 1432 1555 end if
! warning error, itmax = maximum 1433 1556 ! warning error, itmax = maximum
jk = itmax + jk 1434 1557 jk = itmax + jk
ier = 33 1435 1558 ier = 33
end if 1436 1559 end if
if (ic .ne. 0) then 1437 1560 if (ic .ne. 0) then
cycle icloop 1438 1561 cycle icloop
end if 1439 1562 end if
call deflated_work(errcode) 1440 1563 call deflated_work(errcode)
if (errcode == 1) then 1441 1564 if (errcode == 1) then
exit icloop 1442 1565 exit icloop
end if 1443 1566 end if
1444 1567
do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) 1445 1568 do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero)
! take remedial action to induce 1446 1569 ! take remedial action to induce
! convergence 1447 1570 ! convergence
d = d*p5 1448 1571 d = d*p5
h = h*p5 1449 1572 h = h*p5
rt = rt-h 1450 1573 rt = rt-h
call deflated_work(errcode) 1451 1574 call deflated_work(errcode)
if (errcode == 1) then 1452 1575 if (errcode == 1) then
exit icloop 1453 1576 exit icloop
end if 1454 1577 end if
end do 1455 1578 end do
z0 = z1 1456 1579 z0 = z1
z1 = z2 1457 1580 z1 = z2
z2 = fprt 1458 1581 z2 = fprt
y0 = y1 1459 1582 y0 = y1
y1 = y2 1460 1583 y1 = y2
y2 = frt 1461 1584 y2 = frt
end do iloop 1462 1585 end do iloop
end do icloop 1463 1586 end do icloop
x(l) = rt 1464 1587 x(l) = rt
infer(l) = jk 1465 1588 infer(l) = jk
l = l+1 1466 1589 l = l+1
end do rloop 1467 1590 end do rloop
1468 1591
contains 1469 1592 contains
subroutine trans_rt() 1470 1593 subroutine trans_rt()
tem = rten*eps1 1471 1594 tem = rten*eps1
if (cdabs(rt) .gt. ax) tem = tem*rt 1472 1595 if (cdabs(rt) .gt. ax) tem = tem*rt
rt = rt+tem 1473 1596 rt = rt+tem
d = (h+tem)*d/h 1474 1597 d = (h+tem)*d/h
h = h+tem 1475 1598 h = h+tem
end subroutine trans_rt 1476 1599 end subroutine trans_rt
1477 1600
subroutine deflated_work(errcode) 1478 1601 subroutine deflated_work(errcode)
! errcode=0 => no errors 1479 1602 ! errcode=0 => no errors
! errcode=1 => jk>itmax or convergence of second kind achieved 1480 1603 ! errcode=1 => jk>itmax or convergence of second kind achieved
integer :: errcode,flag 1481 1604 integer :: errcode,flag
1482 1605
flag=1 1483 1606 flag=1
loop1: do while(flag==1) 1484 1607 loop1: do while(flag==1)
errcode=0 1485 1608 errcode=0
jk = jk+1 1486 1609 jk = jk+1
if (jk .gt. itmax) then 1487 1610 if (jk .gt. itmax) then
ier=33 1488 1611 ier=33
errcode=1 1489 1612 errcode=1
return 1490 1613 return
end if 1491 1614 end if
frt = f(rt) 1492 1615 frt = f(rt)
fprt = frt 1493 1616 fprt = frt
if (l /= 1) then 1494 1617 if (l /= 1) then
lm1 = l-1 1495 1618 lm1 = l-1
do i=1,lm1 1496 1619 do i=1,lm1
tem = rt - x(i) 1497 1620 tem = rt - x(i)
if (cdabs(tem) .eq. rzero) then 1498 1621 if (cdabs(tem) .eq. rzero) then
!if (ic .ne. 0) go to 15 !! ?? possible? 1499 1622 !if (ic .ne. 0) go to 15 !! ?? possible?
call trans_rt() 1500 1623 call trans_rt()
cycle loop1 1501 1624 cycle loop1
end if 1502 1625 end if
fprt = fprt/tem 1503 1626 fprt = fprt/tem
end do 1504 1627 end do
end if 1505 1628 end if
flag=0 1506 1629 flag=0
end do loop1 1507 1630 end do loop1
1508 1631
if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then 1509 1632 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then
errcode=1 1510 1633 errcode=1
return 1511 1634 return
end if 1512 1635 end if
1513 1636
end subroutine deflated_work 1514 1637 end subroutine deflated_work
1515 1638
end subroutine 1516 1639 end subroutine
1517 1640
1518 1641
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1519 1642 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1520 1643 !
! Integration 1521 1644 ! Integration
! 1522 1645 !
! Only double precision coded atm 1523 1646 ! Only double precision coded atm
! 1524 1647 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1525 1648 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1526 1649
1527 1650
subroutine fvn_d_gauss_legendre(n,qx,qw) 1528 1651 subroutine fvn_d_gauss_legendre(n,qx,qw)
! 1529 1652 !
! This routine compute the n Gauss Legendre abscissas and weights 1530 1653 ! This routine compute the n Gauss Legendre abscissas and weights
! Adapted from Numerical Recipes routine gauleg 1531 1654 ! Adapted from Numerical Recipes routine gauleg
! 1532 1655 !
! n (in) : number of points 1533 1656 ! n (in) : number of points
! qx(out) : abscissas 1534 1657 ! qx(out) : abscissas
! qw(out) : weights 1535 1658 ! qw(out) : weights
! 1536 1659 !
implicit none 1537 1660 implicit none
double precision,parameter :: pi=3.141592653589793d0 1538 1661 double precision,parameter :: pi=3.141592653589793d0
integer, intent(in) :: n 1539 1662 integer, intent(in) :: n
double precision, intent(out) :: qx(n),qw(n) 1540 1663 double precision, intent(out) :: qx(n),qw(n)
1541 1664
integer :: m,i,j 1542 1665 integer :: m,i,j
double precision :: z,z1,p1,p2,p3,pp 1543 1666 double precision :: z,z1,p1,p2,p3,pp
m=(n+1)/2 1544 1667 m=(n+1)/2
do i=1,m 1545 1668 do i=1,m
z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0)) 1546 1669 z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0))
iloop: do 1547 1670 iloop: do
p1=1.d0 1548 1671 p1=1.d0
p2=0.d0 1549 1672 p2=0.d0
do j=1,n 1550 1673 do j=1,n
p3=p2 1551 1674 p3=p2
p2=p1 1552 1675 p2=p1
p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j) 1553 1676 p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j)
end do 1554 1677 end do
pp=dble(n)*(z*p1-p2)/(z*z-1.d0) 1555 1678 pp=dble(n)*(z*p1-p2)/(z*z-1.d0)
z1=z 1556 1679 z1=z
z=z1-p1/pp 1557 1680 z=z1-p1/pp
if (dabs(z-z1)<=epsilon(z)) then 1558 1681 if (dabs(z-z1)<=epsilon(z)) then
exit iloop 1559 1682 exit iloop
end if 1560 1683 end if
end do iloop 1561 1684 end do iloop
qx(i)=-z 1562 1685 qx(i)=-z
qx(n+1-i)=z 1563 1686 qx(n+1-i)=z
qw(i)=2.d0/((1.d0-z*z)*pp*pp) 1564 1687 qw(i)=2.d0/((1.d0-z*z)*pp*pp)
qw(n+1-i)=qw(i) 1565 1688 qw(n+1-i)=qw(i)
end do 1566 1689 end do
end subroutine 1567 1690 end subroutine
1568 1691
1569 1692
1570 1693
subroutine fvn_d_gl_integ(f,a,b,n,res) 1571 1694 subroutine fvn_d_gl_integ(f,a,b,n,res)
! 1572 1695 !
! This is a simple non adaptative integration routine 1573 1696 ! This is a simple non adaptative integration routine
! using n gauss legendre abscissas and weights 1574 1697 ! using n gauss legendre abscissas and weights
! 1575 1698 !
! f(in) : the function to integrate 1576 1699 ! f(in) : the function to integrate
! a(in) : lower bound 1577 1700 ! a(in) : lower bound
! b(in) : higher bound 1578 1701 ! b(in) : higher bound
! n(in) : number of gauss legendre pairs 1579 1702 ! n(in) : number of gauss legendre pairs
! res(out): the evaluation of the integral 1580 1703 ! res(out): the evaluation of the integral
! 1581 1704 !
double precision,external :: f 1582 1705 double precision,external :: f
double precision, intent(in) :: a,b 1583 1706 double precision, intent(in) :: a,b
integer, intent(in):: n 1584 1707 integer, intent(in):: n
double precision, intent(out) :: res 1585 1708 double precision, intent(out) :: res
1586 1709
double precision, allocatable :: qx(:),qw(:) 1587 1710 double precision, allocatable :: qx(:),qw(:)
double precision :: xm,xr 1588 1711 double precision :: xm,xr
integer :: i 1589 1712 integer :: i
1590 1713
! First compute n gauss legendre abs and weight 1591 1714 ! First compute n gauss legendre abs and weight
allocate(qx(n)) 1592 1715 allocate(qx(n))
allocate(qw(n)) 1593 1716 allocate(qw(n))
call fvn_d_gauss_legendre(n,qx,qw) 1594 1717 call fvn_d_gauss_legendre(n,qx,qw)
1595 1718
xm=0.5d0*(b+a) 1596 1719 xm=0.5d0*(b+a)
xr=0.5d0*(b-a) 1597 1720 xr=0.5d0*(b-a)
1598 1721
res=0.d0 1599 1722 res=0.d0
1600 1723
do i=1,n 1601 1724 do i=1,n
res=res+qw(i)*f(xm+xr*qx(i)) 1602 1725 res=res+qw(i)*f(xm+xr*qx(i))
end do 1603 1726 end do
1604 1727
res=xr*res 1605 1728 res=xr*res
1606 1729
deallocate(qw) 1607 1730 deallocate(qw)
deallocate(qx) 1608 1731 deallocate(qx)
1609 1732
end subroutine 1610 1733 end subroutine
1611 1734
!!!!!!!!!!!!!!!!!!!!!!!! 1612 1735 !!!!!!!!!!!!!!!!!!!!!!!!
! 1613 1736 !
! Simple and double adaptative Gauss Kronrod integration based on 1614 1737 ! Simple and double adaptative Gauss Kronrod integration based on
! a modified version of quadpack ( http://www.netlib.org/quadpack 1615 1738 ! a modified version of quadpack ( http://www.netlib.org/quadpack
! 1616 1739 !
! Common parameters : 1617 1740 ! Common parameters :
! 1618 1741 !
! key (in) 1619 1742 ! key (in)
! epsabs 1620 1743 ! epsabs
! epsrel 1621 1744 ! epsrel
! 1622 1745 !
! 1623 1746 !
!!!!!!!!!!!!!!!!!!!!!!!! 1624 1747 !!!!!!!!!!!!!!!!!!!!!!!!
1625 1748
subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit) 1626 1749 subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
! 1627 1750 !
! Evaluate the integral of function f(x) between a and b 1628 1751 ! Evaluate the integral of function f(x) between a and b
! 1629 1752 !
! f(in) : the function 1630 1753 ! f(in) : the function
! a(in) : lower bound 1631 1754 ! a(in) : lower bound
! b(in) : higher bound 1632 1755 ! b(in) : higher bound
! epsabs(in) : desired absolute error 1633 1756 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 1634 1757 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 1635 1758 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 1636 1759 ! 1: 7 - 15 points
! 2: 10 - 21 points 1637 1760 ! 2: 10 - 21 points
! 3: 15 - 31 points 1638 1761 ! 3: 15 - 31 points
! 4: 20 - 41 points 1639 1762 ! 4: 20 - 41 points
! 5: 25 - 51 points 1640 1763 ! 5: 25 - 51 points
! 6: 30 - 61 points 1641 1764 ! 6: 30 - 61 points
! 1642 1765 !
! limit(in) : maximum number of subintervals in the partition of the 1643 1766 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 1644 1767 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 1645 1768 ! behaviour as the imsl routine dqdag
! 1646 1769 !
! res(out) : estimated integral value 1647 1770 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 1648 1771 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 1649 1772 ! ier(out) : error flag from quadpack routines
! 0 : no error 1650 1773 ! 0 : no error
! 1 : maximum number of subdivisions allowed 1651 1774 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 1652 1775 ! has been achieved. one can allow more
! subdivisions by increasing the value of 1653 1776 ! subdivisions by increasing the value of
! limit (and taking the according dimension 1654 1777 ! limit (and taking the according dimension
! adjustments into account). however, if 1655 1778 ! adjustments into account). however, if
! this yield no improvement it is advised 1656 1779 ! this yield no improvement it is advised
! to analyze the integrand in order to 1657 1780 ! to analyze the integrand in order to
! determine the integration difficulaties. 1658 1781 ! determine the integration difficulaties.
! if the position of a local difficulty can 1659 1782 ! if the position of a local difficulty can
! be determined (i.e.singularity, 1660 1783 ! be determined (i.e.singularity,
! discontinuity within the interval) one 1661 1784 ! discontinuity within the interval) one
! will probably gain from splitting up the 1662 1785 ! will probably gain from splitting up the
! interval at this point and calling the 1663 1786 ! interval at this point and calling the
! integrator on the subranges. if possible, 1664 1787 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 1665 1788 ! an appropriate special-purpose integrator
! should be used which is designed for 1666 1789 ! should be used which is designed for
! handling the type of difficulty involved. 1667 1790 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 1668 1791 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 1669 1792 ! detected, which prevents the requested
! tolerance from being achieved. 1670 1793 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 1671 1794 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 1672 1795 ! at some points of the integration
! interval. 1673 1796 ! interval.
! 6 : the input is invalid, because 1674 1797 ! 6 : the input is invalid, because
! (epsabs.le.0 and 1675 1798 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 1676 1799 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 1677 1800 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 1678 1801 ! result, abserr, neval, last are set
! to zero. 1679 1802 ! to zero.
! except when lenw is invalid, iwork(1), 1680 1803 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 1681 1804 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 1682 1805 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 1683 1806 ! work(limit+1) to b.
1684 1807
implicit none 1685 1808 implicit none
double precision, external :: f 1686 1809 double precision, external :: f
double precision, intent(in) :: a,b,epsabs,epsrel 1687 1810 double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key 1688 1811 integer, intent(in) :: key
integer, intent(in) :: limit 1689 1812 integer, intent(in) :: limit
double precision, intent(out) :: res,abserr 1690 1813 double precision, intent(out) :: res,abserr
integer, intent(out) :: ier 1691 1814 integer, intent(out) :: ier
1692 1815
double precision, allocatable :: work(:) 1693 1816 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 1694 1817 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 1695 1818 integer :: lenw,neval,last
1696 1819
! imsl value for limit is 500 1697 1820 ! imsl value for limit is 500
lenw=limit*4 1698 1821 lenw=limit*4
1699 1822
allocate(iwork(limit)) 1700 1823 allocate(iwork(limit))
allocate(work(lenw)) 1701 1824 allocate(work(lenw))
1702 1825
call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 1703 1826 call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
1704 1827
deallocate(work) 1705 1828 deallocate(work)
deallocate(iwork) 1706 1829 deallocate(iwork)
1707 1830
end subroutine 1708 1831 end subroutine
1709 1832
1710 1833
1711 1834
subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit) 1712 1835 subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit)
! 1713 1836 !
! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x) 1714 1837 ! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x)
! 1715 1838 !
! f(in) : the function 1716 1839 ! f(in) : the function
! a(in) : lower bound 1717 1840 ! a(in) : lower bound
! b(in) : higher bound 1718 1841 ! b(in) : higher bound
! g(in) : external function describing lower bound for y 1719 1842 ! g(in) : external function describing lower bound for y
! h(in) : external function describing higher bound for y 1720 1843 ! h(in) : external function describing higher bound for y
! epsabs(in) : desired absolute error 1721 1844 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 1722 1845 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 1723 1846 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 1724 1847 ! 1: 7 - 15 points
! 2: 10 - 21 points 1725 1848 ! 2: 10 - 21 points
! 3: 15 - 31 points 1726 1849 ! 3: 15 - 31 points
! 4: 20 - 41 points 1727 1850 ! 4: 20 - 41 points
! 5: 25 - 51 points 1728 1851 ! 5: 25 - 51 points
! 6: 30 - 61 points 1729 1852 ! 6: 30 - 61 points
! 1730 1853 !
! limit(in) : maximum number of subintervals in the partition of the 1731 1854 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 1732 1855 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 1733 1856 ! behaviour as the imsl routine dqdag
! 1734 1857 !
! res(out) : estimated integral value 1735 1858 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 1736 1859 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 1737 1860 ! ier(out) : error flag from quadpack routines
! 0 : no error 1738 1861 ! 0 : no error
! 1 : maximum number of subdivisions allowed 1739 1862 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 1740 1863 ! has been achieved. one can allow more
! subdivisions by increasing the value of 1741 1864 ! subdivisions by increasing the value of
! limit (and taking the according dimension 1742 1865 ! limit (and taking the according dimension
! adjustments into account). however, if 1743 1866 ! adjustments into account). however, if
! this yield no improvement it is advised 1744 1867 ! this yield no improvement it is advised
! to analyze the integrand in order to 1745 1868 ! to analyze the integrand in order to
! determine the integration difficulaties. 1746 1869 ! determine the integration difficulaties.
! if the position of a local difficulty can 1747 1870 ! if the position of a local difficulty can
! be determined (i.e.singularity, 1748 1871 ! be determined (i.e.singularity,
! discontinuity within the interval) one 1749 1872 ! discontinuity within the interval) one
! will probably gain from splitting up the 1750 1873 ! will probably gain from splitting up the
! interval at this point and calling the 1751 1874 ! interval at this point and calling the
! integrator on the subranges. if possible, 1752 1875 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 1753 1876 ! an appropriate special-purpose integrator
! should be used which is designed for 1754 1877 ! should be used which is designed for
! handling the type of difficulty involved. 1755 1878 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 1756 1879 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 1757 1880 ! detected, which prevents the requested
! tolerance from being achieved. 1758 1881 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 1759 1882 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 1760 1883 ! at some points of the integration
! interval. 1761 1884 ! interval.
! 6 : the input is invalid, because 1762 1885 ! 6 : the input is invalid, because
! (epsabs.le.0 and 1763 1886 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 1764 1887 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 1765 1888 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 1766 1889 ! result, abserr, neval, last are set
! to zero. 1767 1890 ! to zero.
! except when lenw is invalid, iwork(1), 1768 1891 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 1769 1892 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 1770 1893 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 1771 1894 ! work(limit+1) to b.
1772 1895
implicit none 1773 1896 implicit none
double precision, external:: f,g,h 1774 1897 double precision, external:: f,g,h
double precision, intent(in) :: a,b,epsabs,epsrel 1775 1898 double precision, intent(in) :: a,b,epsabs,epsrel
integer, intent(in) :: key,limit 1776 1899 integer, intent(in) :: key,limit
integer, intent(out) :: ier 1777 1900 integer, intent(out) :: ier
double precision, intent(out) :: res,abserr 1778 1901 double precision, intent(out) :: res,abserr
1779 1902
1780 1903
double precision, allocatable :: work(:) 1781 1904 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 1782 1905 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 1783 1906 integer :: lenw,neval,last
1784 1907
! imsl value for limit is 500 1785 1908 ! imsl value for limit is 500
lenw=limit*4 1786 1909 lenw=limit*4
allocate(work(lenw)) 1787 1910 allocate(work(lenw))
allocate(iwork(limit)) 1788 1911 allocate(iwork(limit))
1789 1912
call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 1790 1913 call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
1791 1914
deallocate(iwork) 1792 1915 deallocate(iwork)
deallocate(work) 1793 1916 deallocate(work)
end subroutine 1794 1917 end subroutine
1795 1918
1796 1919
1797 1920
subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit) 1798 1921 subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit)
! 1799 1922 !
! Evaluate the single integral of function f(x,y) for y between a and b with a 1800 1923 ! Evaluate the single integral of function f(x,y) for y between a and b with a
! given x value 1801 1924 ! given x value
! 1802 1925 !
! This function is used for the evaluation of the double integral fvn_d_integ_2_gk 1803 1926 ! This function is used for the evaluation of the double integral fvn_d_integ_2_gk
! 1804 1927 !
! f(in) : the function 1805 1928 ! f(in) : the function
! x(in) : x 1806 1929 ! x(in) : x
! a(in) : lower bound 1807 1930 ! a(in) : lower bound
! b(in) : higher bound 1808 1931 ! b(in) : higher bound
! epsabs(in) : desired absolute error 1809 1932 ! epsabs(in) : desired absolute error
! epsrel(in) : desired relative error 1810 1933 ! epsrel(in) : desired relative error
! key(in) : gauss kronrod rule 1811 1934 ! key(in) : gauss kronrod rule
! 1: 7 - 15 points 1812 1935 ! 1: 7 - 15 points
! 2: 10 - 21 points 1813 1936 ! 2: 10 - 21 points
! 3: 15 - 31 points 1814 1937 ! 3: 15 - 31 points
! 4: 20 - 41 points 1815 1938 ! 4: 20 - 41 points
! 5: 25 - 51 points 1816 1939 ! 5: 25 - 51 points
! 6: 30 - 61 points 1817 1940 ! 6: 30 - 61 points
! 1818 1941 !
! limit(in) : maximum number of subintervals in the partition of the 1819 1942 ! limit(in) : maximum number of subintervals in the partition of the
! given integration interval (a,b). A value of 500 will give the same 1820 1943 ! given integration interval (a,b). A value of 500 will give the same
! behaviour as the imsl routine dqdag 1821 1944 ! behaviour as the imsl routine dqdag
! 1822 1945 !
! res(out) : estimated integral value 1823 1946 ! res(out) : estimated integral value
! abserr(out) : estimated absolute error 1824 1947 ! abserr(out) : estimated absolute error
! ier(out) : error flag from quadpack routines 1825 1948 ! ier(out) : error flag from quadpack routines
! 0 : no error 1826 1949 ! 0 : no error
! 1 : maximum number of subdivisions allowed 1827 1950 ! 1 : maximum number of subdivisions allowed
! has been achieved. one can allow more 1828 1951 ! has been achieved. one can allow more
! subdivisions by increasing the value of 1829 1952 ! subdivisions by increasing the value of
! limit (and taking the according dimension 1830 1953 ! limit (and taking the according dimension
! adjustments into account). however, if 1831 1954 ! adjustments into account). however, if
! this yield no improvement it is advised 1832 1955 ! this yield no improvement it is advised
! to analyze the integrand in order to 1833 1956 ! to analyze the integrand in order to
! determine the integration difficulaties. 1834 1957 ! determine the integration difficulaties.
! if the position of a local difficulty can 1835 1958 ! if the position of a local difficulty can
! be determined (i.e.singularity, 1836 1959 ! be determined (i.e.singularity,
! discontinuity within the interval) one 1837 1960 ! discontinuity within the interval) one
! will probably gain from splitting up the 1838 1961 ! will probably gain from splitting up the
! interval at this point and calling the 1839 1962 ! interval at this point and calling the
! integrator on the subranges. if possible, 1840 1963 ! integrator on the subranges. if possible,
! an appropriate special-purpose integrator 1841 1964 ! an appropriate special-purpose integrator
! should be used which is designed for 1842 1965 ! should be used which is designed for
! handling the type of difficulty involved. 1843 1966 ! handling the type of difficulty involved.
! 2 : the occurrence of roundoff error is 1844 1967 ! 2 : the occurrence of roundoff error is
! detected, which prevents the requested 1845 1968 ! detected, which prevents the requested
! tolerance from being achieved. 1846 1969 ! tolerance from being achieved.
! 3 : extremely bad integrand behaviour occurs 1847 1970 ! 3 : extremely bad integrand behaviour occurs
! at some points of the integration 1848 1971 ! at some points of the integration
! interval. 1849 1972 ! interval.
! 6 : the input is invalid, because 1850 1973 ! 6 : the input is invalid, because
! (epsabs.le.0 and 1851 1974 ! (epsabs.le.0 and
! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) 1852 1975 ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28))
! or limit.lt.1 or lenw.lt.limit*4. 1853 1976 ! or limit.lt.1 or lenw.lt.limit*4.
! result, abserr, neval, last are set 1854 1977 ! result, abserr, neval, last are set
! to zero. 1855 1978 ! to zero.
! except when lenw is invalid, iwork(1), 1856 1979 ! except when lenw is invalid, iwork(1),
! work(limit*2+1) and work(limit*3+1) are 1857 1980 ! work(limit*2+1) and work(limit*3+1) are
! set to zero, work(1) is set to a and 1858 1981 ! set to zero, work(1) is set to a and
! work(limit+1) to b. 1859 1982 ! work(limit+1) to b.
1860 1983
implicit none 1861 1984 implicit none
double precision, external:: f 1862 1985 double precision, external:: f
double precision, intent(in) :: x,a,b,epsabs,epsrel 1863 1986 double precision, intent(in) :: x,a,b,epsabs,epsrel
integer, intent(in) :: key,limit 1864 1987 integer, intent(in) :: key,limit
integer, intent(out) :: ier 1865 1988 integer, intent(out) :: ier
double precision, intent(out) :: res,abserr 1866 1989 double precision, intent(out) :: res,abserr
1867 1990
1868 1991
double precision, allocatable :: work(:) 1869 1992 double precision, allocatable :: work(:)
integer, allocatable :: iwork(:) 1870 1993 integer, allocatable :: iwork(:)
integer :: lenw,neval,last 1871 1994 integer :: lenw,neval,last
1872 1995
! imsl value for limit is 500 1873 1996 ! imsl value for limit is 500
lenw=limit*4 1874 1997 lenw=limit*4
allocate(work(lenw)) 1875 1998 allocate(work(lenw))
allocate(iwork(limit)) 1876 1999 allocate(iwork(limit))
1877 2000
call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) 1878 2001 call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work)
1879 2002
deallocate(iwork) 1880 2003 deallocate(iwork)
deallocate(work) 1881 2004 deallocate(work)
end subroutine 1882 2005 end subroutine
1883 2006
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1884 2007 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Include the modified quadpack files 1885 2008 ! Include the modified quadpack files
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1886 2009 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
include "fvn_quadpack/dqag_2d_inner.f" 1887 2010 include "fvn_quadpack/dqag_2d_inner.f"
include "fvn_quadpack/dqk15_2d_inner.f" 1888 2011 include "fvn_quadpack/dqk15_2d_inner.f"
include "fvn_quadpack/dqk31_2d_outer.f" 1889 2012 include "fvn_quadpack/dqk31_2d_outer.f"
include "fvn_quadpack/d1mach.f" 1890 2013 include "fvn_quadpack/d1mach.f"
include "fvn_quadpack/dqk31_2d_inner.f" 1891 2014 include "fvn_quadpack/dqk31_2d_inner.f"
include "fvn_quadpack/dqage.f" 1892 2015 include "fvn_quadpack/dqage.f"
include "fvn_quadpack/dqk15.f" 1893 2016 include "fvn_quadpack/dqk15.f"
include "fvn_quadpack/dqk21.f" 1894 2017 include "fvn_quadpack/dqk21.f"
include "fvn_quadpack/dqk31.f" 1895 2018 include "fvn_quadpack/dqk31.f"
include "fvn_quadpack/dqk41.f" 1896 2019 include "fvn_quadpack/dqk41.f"
include "fvn_quadpack/dqk51.f" 1897 2020 include "fvn_quadpack/dqk51.f"
include "fvn_quadpack/dqk61.f" 1898 2021 include "fvn_quadpack/dqk61.f"
include "fvn_quadpack/dqk41_2d_outer.f" 1899 2022 include "fvn_quadpack/dqk41_2d_outer.f"
include "fvn_quadpack/dqk41_2d_inner.f" 1900 2023 include "fvn_quadpack/dqk41_2d_inner.f"
include "fvn_quadpack/dqag_2d_outer.f" 1901 2024 include "fvn_quadpack/dqag_2d_outer.f"
include "fvn_quadpack/dqpsrt.f" 1902 2025 include "fvn_quadpack/dqpsrt.f"
include "fvn_quadpack/dqag.f" 1903 2026 include "fvn_quadpack/dqag.f"
include "fvn_quadpack/dqage_2d_outer.f" 1904 2027 include "fvn_quadpack/dqage_2d_outer.f"
include "fvn_quadpack/dqage_2d_inner.f" 1905 2028 include "fvn_quadpack/dqage_2d_inner.f"
include "fvn_quadpack/dqk51_2d_outer.f" 1906 2029 include "fvn_quadpack/dqk51_2d_outer.f"
include "fvn_quadpack/dqk51_2d_inner.f" 1907 2030 include "fvn_quadpack/dqk51_2d_inner.f"
include "fvn_quadpack/dqk61_2d_outer.f" 1908 2031 include "fvn_quadpack/dqk61_2d_outer.f"
include "fvn_quadpack/dqk21_2d_outer.f" 1909 2032 include "fvn_quadpack/dqk21_2d_outer.f"
include "fvn_quadpack/dqk61_2d_inner.f" 1910 2033 include "fvn_quadpack/dqk61_2d_inner.f"
include "fvn_quadpack/dqk21_2d_inner.f" 1911 2034 include "fvn_quadpack/dqk21_2d_inner.f"
include "fvn_quadpack/dqk15_2d_outer.f" 1912 2035 include "fvn_quadpack/dqk15_2d_outer.f"
1913 2036
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1914 2037 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 1915 2038 !
! Trigonometric functions 1916 2039 ! Trigonometric functions
! 1917 2040 !
! fvn_z_acos, fvn_z_asin : complex arc cosine and sine 1918 2041 ! fvn_z_acos, fvn_z_asin : complex arc cosine and sine
! fvn_d_acosh : arc cosinus hyperbolic 1919 2042 ! fvn_d_acosh : arc cosinus hyperbolic
! fvn_d_asinh : arc sinus hyperbolic 1920 2043 ! fvn_d_asinh : arc sinus hyperbolic
! 1921 2044 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1922 2045 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
function fvn_z_acos(z) 1923 2046 function fvn_z_acos(z)
! double complex arccos function adapted from 1924 2047 ! double complex arccos function adapted from
! the c gsl library 1925 2048 ! the c gsl library
! http://www.gnu.org/software/gsl/ 1926 2049 ! http://www.gnu.org/software/gsl/
implicit none 1927 2050 implicit none
complex(kind=8) :: fvn_z_acos 1928 2051 complex(kind=8) :: fvn_z_acos
complex(kind=8) :: z 1929 2052 complex(kind=8) :: z
real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1 1930 2053 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 1931 2054 real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
complex(kind=8),parameter :: i=(0._8,1._8) 1932 2055 complex(kind=8),parameter :: i=(0._8,1._8)
real(kind=8) :: r_res,i_res 1933 2056 real(kind=8) :: r_res,i_res
1934 2057
rz=dreal(z) 1935 2058 rz=dreal(z)
iz=dimag(z) 1936 2059 iz=dimag(z)
if ( iz == 0._8 ) then 1937 2060 if ( iz == 0._8 ) then
fvn_z_acos=fvn_z_acos_real(rz) 1938 2061 fvn_z_acos=fvn_z_acos_real(rz)
return 1939 2062 return
end if 1940 2063 end if
1941 2064
x=dabs(rz) 1942 2065 x=dabs(rz)
y=dabs(iz) 1943 2066 y=dabs(iz)
r=fvn_d_hypot(x+1.,y) 1944 2067 r=fvn_d_hypot(x+1.,y)
s=fvn_d_hypot(x-1.,y) 1945 2068 s=fvn_d_hypot(x-1.,y)
a=0.5*(r + s) 1946 2069 a=0.5*(r + s)
b=x/a 1947 2070 b=x/a
y2=y*y 1948 2071 y2=y*y
1949 2072
if (b <= b_crossover) then 1950 2073 if (b <= b_crossover) then
r_res=dacos(b) 1951 2074 r_res=dacos(b)
else 1952 2075 else
if (x <= 1.) then 1953 2076 if (x <= 1.) then
d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x))) 1954 2077 d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
r_res=datan(dsqrt(d)/x) 1955 2078 r_res=datan(dsqrt(d)/x)
else 1956 2079 else
apx=a+x 1957 2080 apx=a+x
d=0.5*(apx/(r+x+1)+apx/(s + (x - 1))) 1958 2081 d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
r_res=datan((y*dsqrt(d))/x); 1959 2082 r_res=datan((y*dsqrt(d))/x);
end if 1960 2083 end if
end if 1961 2084 end if
1962 2085
if (a <= a_crossover) then 1963 2086 if (a <= a_crossover) then
if (x < 1.) then 1964 2087 if (x < 1.) then
am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x))) 1965 2088 am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
else 1966 2089 else
am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1))) 1967 2090 am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
end if 1968 2091 end if
i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1)))); 1969 2092 i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
else 1970 2093 else
i_res = dlog (a + dsqrt (a*a - 1.)); 1971 2094 i_res = dlog (a + dsqrt (a*a - 1.));
end if 1972 2095 end if
if (rz <0.) then 1973 2096 if (rz <0.) then
r_res=fvn_pi-r_res 1974 2097 r_res=fvn_pi-r_res
end if 1975 2098 end if
i_res=-sign(1._8,iz)*i_res 1976 2099 i_res=-sign(1._8,iz)*i_res
fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res) 1977 2100 fvn_z_acos=dcmplx(r_res)+fvn_i*dcmplx(i_res)
1978 2101
end function fvn_z_acos 1979 2102 end function fvn_z_acos
1980 2103
function fvn_z_acos_real(r) 1981 2104 function fvn_z_acos_real(r)
! return the double complex arc cosinus for a 1982 2105 ! return the double complex arc cosinus for a
! double precision argument 1983 2106 ! double precision argument
implicit none 1984 2107 implicit none
real(kind=8) :: r 1985 2108 real(kind=8) :: r
complex(kind=8) :: fvn_z_acos_real 1986 2109 complex(kind=8) :: fvn_z_acos_real
1987 2110
if (dabs(r)<=1._8) then 1988 2111 if (dabs(r)<=1._8) then
fvn_z_acos_real=dcmplx(dacos(r)) 1989 2112 fvn_z_acos_real=dcmplx(dacos(r))
return 1990 2113 return
end if 1991 2114 end if
if (r < 0._8) then 1992 2115 if (r < 0._8) then
fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r)) 1993 2116 fvn_z_acos_real=dcmplx(fvn_pi)-fvn_i*dcmplx(fvn_d_acosh(-r))
else 1994 2117 else
fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r)) 1995 2118 fvn_z_acos_real=fvn_i*dcmplx(fvn_d_acosh(r))
end if 1996 2119 end if
end function 1997 2120 end function
1998 2121
1999 2122
function fvn_z_asin(z) 2000 2123 function fvn_z_asin(z)
! double complex arcsin function derived from 2001 2124 ! double complex arcsin function derived from
! the c gsl library 2002 2125 ! the c gsl library
! http://www.gnu.org/software/gsl/ 2003 2126 ! http://www.gnu.org/software/gsl/
implicit none 2004 2127 implicit none
complex(kind=8) :: fvn_z_asin 2005 2128 complex(kind=8) :: fvn_z_asin
complex(kind=8) :: z 2006 2129 complex(kind=8) :: z
real(kind=8) :: rz,iz,x,y,a,b,y2,r,s,d,apx,am1 2007 2130 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 2008 2131 real(kind=8),parameter :: a_crossover=1.5_8,b_crossover = 0.6417_8
real(kind=8) :: r_res,i_res 2009 2132 real(kind=8) :: r_res,i_res
2010 2133
rz=dreal(z) 2011 2134 rz=dreal(z)
iz=dimag(z) 2012 2135 iz=dimag(z)
if ( iz == 0._8 ) then 2013 2136 if ( iz == 0._8 ) then
! z is real 2014 2137 ! z is real
fvn_z_asin=fvn_z_asin_real(rz) 2015 2138 fvn_z_asin=fvn_z_asin_real(rz)
return 2016 2139 return
end if 2017 2140 end if
2018 2141
x=dabs(rz) 2019 2142 x=dabs(rz)
y=dabs(iz) 2020 2143 y=dabs(iz)
r=fvn_d_hypot(x+1.,y) 2021 2144 r=fvn_d_hypot(x+1.,y)
s=fvn_d_hypot(x-1.,y) 2022 2145 s=fvn_d_hypot(x-1.,y)
a=0.5*(r + s) 2023 2146 a=0.5*(r + s)
b=x/a 2024 2147 b=x/a
y2=y*y 2025 2148 y2=y*y
2026 2149
if (b <= b_crossover) then 2027 2150 if (b <= b_crossover) then
r_res=dasin(b) 2028 2151 r_res=dasin(b)
else 2029 2152 else
if (x <= 1.) then 2030 2153 if (x <= 1.) then
d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x))) 2031 2154 d=0.5*(a+x)*(y2/(r+x+1)+(s + (1 - x)))
r_res=datan(x/dsqrt(d)) 2032 2155 r_res=datan(x/dsqrt(d))
else 2033 2156 else
apx=a+x 2034 2157 apx=a+x
d=0.5*(apx/(r+x+1)+apx/(s + (x - 1))) 2035 2158 d=0.5*(apx/(r+x+1)+apx/(s + (x - 1)))
r_res=datan(x/(y*dsqrt(d))); 2036 2159 r_res=datan(x/(y*dsqrt(d)));
end if 2037 2160 end if
end if 2038 2161 end if
2039 2162
if (a <= a_crossover) then 2040 2163 if (a <= a_crossover) then
if (x < 1.) then 2041 2164 if (x < 1.) then
am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x))) 2042 2165 am1=0.5*(y2 / (r + (x + 1)) + y2 / (s + (1 - x)))
else 2043 2166 else
am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1))) 2044 2167 am1=0.5*(y2 / (r + (x + 1)) + (s + (x - 1)))
end if 2045 2168 end if
i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1)))); 2046 2169 i_res = dlog(1.+(am1 + sqrt (am1 * (a + 1))));
else 2047 2170 else
i_res = dlog (a + dsqrt (a*a - 1.)); 2048 2171 i_res = dlog (a + dsqrt (a*a - 1.));
end if 2049 2172 end if
r_res=sign(1._8,rz)*r_res 2050 2173 r_res=sign(1._8,rz)*r_res
i_res=sign(1._8,iz)*i_res 2051 2174 i_res=sign(1._8,iz)*i_res
fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res) 2052 2175 fvn_z_asin=dcmplx(r_res)+fvn_i*dcmplx(i_res)
2053 2176
end function fvn_z_asin 2054 2177 end function fvn_z_asin
2055 2178
function fvn_z_asin_real(r) 2056 2179 function fvn_z_asin_real(r)
! return the double complex arc sinus for a 2057 2180 ! return the double complex arc sinus for a
! double precision argument 2058 2181 ! double precision argument
implicit none 2059 2182 implicit none
real(kind=8) :: r 2060 2183 real(kind=8) :: r
complex(kind=8) :: fvn_z_asin_real 2061 2184 complex(kind=8) :: fvn_z_asin_real
2062 2185
if (dabs(r)<=1._8) then 2063 2186 if (dabs(r)<=1._8) then
fvn_z_asin_real=dcmplx(dasin(r)) 2064 2187 fvn_z_asin_real=dcmplx(dasin(r))
return 2065 2188 return
end if 2066 2189 end if
if (r < 0._8) then 2067 2190 if (r < 0._8) then
fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r)) 2068 2191 fvn_z_asin_real=dcmplx(-fvn_pi/2._8)+fvn_i*dcmplx(fvn_d_acosh(-r))
else 2069 2192 else
fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r)) 2070 2193 fvn_z_asin_real=dcmplx(fvn_pi/2._8)-fvn_i*dcmplx(fvn_d_acosh(r))
end if 2071 2194 end if
end function fvn_z_asin_real 2072 2195 end function fvn_z_asin_real
2073 2196
function fvn_d_acosh(r) 2074 2197 function fvn_d_acosh(r)
! return the arc hyperbolic cosine 2075 2198 ! return the arc hyperbolic cosine
implicit none 2076 2199 implicit none
real(kind=8) :: r 2077 2200 real(kind=8) :: r
real(kind=8) :: fvn_d_acosh 2078 2201 real(kind=8) :: fvn_d_acosh
if (r >=1) then 2079 2202 if (r >=1) then
fvn_d_acosh=dlog(r+dsqrt(r*r-1)) 2080 2203 fvn_d_acosh=dlog(r+dsqrt(r*r-1))
else 2081 2204 else
!! TODO : Better error handling!!!!!! 2082 2205 !! TODO : Better error handling!!!!!!
stop "Argument to fvn_d_acosh lesser than 1" 2083 2206 stop "Argument to fvn_d_acosh lesser than 1"
end if 2084 2207 end if
end function fvn_d_acosh 2085 2208 end function fvn_d_acosh
2086 2209
function fvn_d_asinh(r) 2087 2210 function fvn_d_asinh(r)
! return the arc hyperbolic sine 2088 2211 ! return the arc hyperbolic sine
implicit none 2089 2212 implicit none
real(kind=8) :: r 2090 2213 real(kind=8) :: r
real(kind=8) :: fvn_d_asinh 2091 2214 real(kind=8) :: fvn_d_asinh
fvn_d_asinh=dlog(r+dsqrt(r*r+1)) 2092 2215 fvn_d_asinh=dlog(r+dsqrt(r*r+1))
end function fvn_d_asinh 2093 2216 end function fvn_d_asinh
2094 2217
function fvn_d_hypot(a,b) 2095 2218 function fvn_d_hypot(a,b)
implicit none 2096 2219 implicit none
! return the euclidian norm of vector(a,b) 2097 2220 ! return the euclidian norm of vector(a,b)
real(kind=8) :: a,b 2098 2221 real(kind=8) :: a,b
real(kind=8) :: fvn_d_hypot 2099 2222 real(kind=8) :: fvn_d_hypot
fvn_d_hypot=dsqrt(a*a+b*b) 2100 2223 fvn_d_hypot=dsqrt(a*a+b*b)
end function 2101 2224 end function
2102 2225
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2103 2226 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 2104 2227 !
! SPARSE RESOLUTION 2105 2228 ! SPARSE RESOLUTION
! 2106 2229 !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 2107 2230 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! 2108 2231 !
! Sparse resolution is done by interfaçing Tim Davi's UMFPACK 2109 2232 ! Sparse resolution is done by interfaçing Tim Davi's UMFPACK
! http://www.cise.ufl.edu/research/sparse/SuiteSparse/ 2110 2233 ! http://www.cise.ufl.edu/research/sparse/SuiteSparse/
! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig 2111 2234 ! Used packages from SuiteSparse : AMD,UMFPACK,UFconfig
! 2112 2235 !
! Solve Ax=B using UMFPACK 2113 2236 ! Solve Ax=B using UMFPACK
! 2114 2237 !
! Where A is a sparse matrix given in its triplet form 2115 2238 ! Where A is a sparse matrix given in its triplet form
! T -> non zero elements 2116 2239 ! T -> non zero elements
! Ti,Tj -> row and column index (1-based) of the given elt 2117 2240 ! Ti,Tj -> row and column index (1-based) of the given elt
! n : rank of matrix A 2118 2241 ! n : rank of matrix A
! nz : number of non zero elts 2119 2242 ! nz : number of non zero elts
! 2120 2243 !
! fvn_*_sparse_solve 2121 2244 ! fvn_*_sparse_solve
! * = zl : double complex + integer(8) 2122 2245 ! * = zl : double complex + integer(8)
! * = zi : double complex + integer(4) 2123 2246 ! * = zi : double complex + integer(4)
! 2124 2247 !
subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status) 2125 2248 subroutine fvn_zl_sparse_solve(n,nz,T,Ti,Tj,B,x,status)
implicit none 2126 2249 implicit none
integer(8), intent(in) :: n,nz 2127 2250 integer(8), intent(in) :: n,nz
complex(8),dimension(nz),intent(in) :: T 2128 2251 complex(8),dimension(nz),intent(in) :: T
integer(8),dimension(nz),intent(in) :: Ti,Tj 2129 2252 integer(8),dimension(nz),intent(in) :: Ti,Tj
complex(8),dimension(n),intent(in) :: B 2130 2253 complex(8),dimension(n),intent(in) :: B
complex(8),dimension(n),intent(out) :: x 2131 2254 complex(8),dimension(n),intent(out) :: x
integer(8), intent(out) :: status 2132 2255 integer(8), intent(out) :: status
2133 2256
integer(8),dimension(:),allocatable :: wTi,wTj 2134 2257 integer(8),dimension(:),allocatable :: wTi,wTj
real(8),dimension(:),allocatable :: Tx,Tz 2135 2258 real(8),dimension(:),allocatable :: Tx,Tz
real(8),dimension(:),allocatable :: Ax,Az 2136 2259 real(8),dimension(:),allocatable :: Ax,Az
integer(8),dimension(:),allocatable :: Ap,Ai 2137 2260 integer(8),dimension(:),allocatable :: Ap,Ai
integer(8) :: symbolic,numeric 2138 2261 integer(8) :: symbolic,numeric
real(8),dimension(:),allocatable :: xx,xz,bx,bz 2139 2262 real(8),dimension(:),allocatable :: xx,xz,bx,bz
real(8),dimension(90) :: info 2140 2263 real(8),dimension(90) :: info
real(8),dimension(20) :: control 2141 2264 real(8),dimension(20) :: control
integer(8) :: sys 2142 2265 integer(8) :: sys
2143 2266
2144 2267
status=0 2145 2268 status=0
2146 2269
! we use a working copy of Ti and Tj to perform 1-based to 0-based translation 2147 2270 ! we use a working copy of Ti and Tj to perform 1-based to 0-based translation
! Tx and Tz are the real and imaginary parts of T 2148 2271 ! Tx and Tz are the real and imaginary parts of T
allocate(wTi(nz),wTj(nz)) 2149 2272 allocate(wTi(nz),wTj(nz))
allocate(Tx(nz),Tz(nz)) 2150 2273 allocate(Tx(nz),Tz(nz))
Tx=dble(T) 2151 2274 Tx=dble(T)
Tz=aimag(T) 2152 2275 Tz=aimag(T)
wTi=Ti-1 2153 2276 wTi=Ti-1