lxsd / fvn

Download as
Browse Code ยป

Commit 6ac82e990ee0ab81340d03178f3bd7aec3d7de43

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

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

Showing 1 changed file with 0 additions and 1779 deletions Inline Diff

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