lxsd / fvn

Download as
Browse Code ยป

Commit d32a4703383b961c2182b5ea68585cf33b41060b

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

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

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

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