module fvn !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! fvn : a f95 module replacement for some imsl routines ! it uses lapack for linear algebra ! it uses modified quadpack for integration ! ! William Daniau 2007 ! william.daniau@femto-st.fr ! ! Routines naming scheme : ! ! fvn_x_name ! where x can be s : real ! d : real double precision ! c : complex ! z : double complex ! ! ! 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 ;-) ! ! Version 1.1 ! ! TO DO LIST : ! + Order eigenvalues and vectors in decreasing eigenvalue's modulus order -> atm ! eigenvalues are given with no particular order. ! + Generic interface for fvn_x_name family -> fvn_name ! + Make some parameters optional, status for example ! + use f95 kinds "double complex" -> complex(kind=8) ! + unify quadpack routines ! + ... ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! implicit none ! All quadpack routines are private to the module 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, & 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 contains !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Matrix inversion subroutines ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine fvn_s_matinv(d,a,inva,status) ! ! Matrix inversion of a real matrix using BLAS and LAPACK ! ! d (in) : matrix rank ! a (in) : input matrix ! inva (out) : inversed matrix ! status (ou) : =0 if something failed ! integer, intent(in) :: d real, intent(in) :: a(d,d) real, intent(out) :: inva(d,d) integer, intent(out) :: status integer, allocatable :: ipiv(:) real, allocatable :: work(:) real twork(1) integer :: info integer :: lwork status=1 allocate(ipiv(d)) ! copy a into inva using BLAS !call scopy(d*d,a,1,inva,1) inva(:,:)=a(:,:) ! LU factorization using LAPACK call sgetrf(d,d,inva,d,ipiv,info) ! if info is not equal to 0, something went wrong we exit setting status to 0 if (info /= 0) then status=0 deallocate(ipiv) return end if ! we use the query fonction of xxxtri to obtain the optimal workspace size call sgetri(d,inva,d,ipiv,twork,-1,info) lwork=int(twork(1)) allocate(work(lwork)) ! Matrix inversion using LAPACK call sgetri(d,inva,d,ipiv,work,lwork,info) ! again if info is not equal to 0, we exit setting status to 0 if (info /= 0) then status=0 end if deallocate(work) deallocate(ipiv) end subroutine subroutine fvn_d_matinv(d,a,inva,status) ! ! Matrix inversion of a double precision matrix using BLAS and LAPACK ! ! d (in) : matrix rank ! a (in) : input matrix ! inva (out) : inversed matrix ! status (ou) : =0 if something failed ! integer, intent(in) :: d double precision, intent(in) :: a(d,d) double precision, intent(out) :: inva(d,d) integer, intent(out) :: status integer, allocatable :: ipiv(:) double precision, allocatable :: work(:) double precision :: twork(1) integer :: info integer :: lwork status=1 allocate(ipiv(d)) ! copy a into inva using BLAS !call dcopy(d*d,a,1,inva,1) inva(:,:)=a(:,:) ! LU factorization using LAPACK call dgetrf(d,d,inva,d,ipiv,info) ! if info is not equal to 0, something went wrong we exit setting status to 0 if (info /= 0) then status=0 deallocate(ipiv) return end if ! we use the query fonction of xxxtri to obtain the optimal workspace size call dgetri(d,inva,d,ipiv,twork,-1,info) lwork=int(twork(1)) allocate(work(lwork)) ! Matrix inversion using LAPACK call dgetri(d,inva,d,ipiv,work,lwork,info) ! again if info is not equal to 0, we exit setting status to 0 if (info /= 0) then status=0 end if deallocate(work) deallocate(ipiv) end subroutine subroutine fvn_c_matinv(d,a,inva,status) ! ! Matrix inversion of a complex matrix using BLAS and LAPACK ! ! d (in) : matrix rank ! a (in) : input matrix ! inva (out) : inversed matrix ! status (ou) : =0 if something failed ! integer, intent(in) :: d complex, intent(in) :: a(d,d) complex, intent(out) :: inva(d,d) integer, intent(out) :: status integer, allocatable :: ipiv(:) complex, allocatable :: work(:) complex :: twork(1) integer :: info integer :: lwork status=1 allocate(ipiv(d)) ! copy a into inva using BLAS !call ccopy(d*d,a,1,inva,1) inva(:,:)=a(:,:) ! LU factorization using LAPACK call cgetrf(d,d,inva,d,ipiv,info) ! if info is not equal to 0, something went wrong we exit setting status to 0 if (info /= 0) then status=0 deallocate(ipiv) return end if ! we use the query fonction of xxxtri to obtain the optimal workspace size call cgetri(d,inva,d,ipiv,twork,-1,info) lwork=int(twork(1)) allocate(work(lwork)) ! Matrix inversion using LAPACK call cgetri(d,inva,d,ipiv,work,lwork,info) ! again if info is not equal to 0, we exit setting status to 0 if (info /= 0) then status=0 end if deallocate(work) deallocate(ipiv) end subroutine subroutine fvn_z_matinv(d,a,inva,status) ! ! Matrix inversion of a double complex matrix using BLAS and LAPACK ! ! d (in) : matrix rank ! a (in) : input matrix ! inva (out) : inversed matrix ! status (ou) : =0 if something failed ! integer, intent(in) :: d double complex, intent(in) :: a(d,d) double complex, intent(out) :: inva(d,d) integer, intent(out) :: status integer, allocatable :: ipiv(:) double complex, allocatable :: work(:) double complex :: twork(1) integer :: info integer :: lwork status=1 allocate(ipiv(d)) ! copy a into inva using BLAS !call zcopy(d*d,a,1,inva,1) inva(:,:)=a(:,:) ! LU factorization using LAPACK call zgetrf(d,d,inva,d,ipiv,info) ! if info is not equal to 0, something went wrong we exit setting status to 0 if (info /= 0) then status=0 deallocate(ipiv) return end if ! we use the query fonction of xxxtri to obtain the optimal workspace size call zgetri(d,inva,d,ipiv,twork,-1,info) lwork=int(twork(1)) allocate(work(lwork)) ! Matrix inversion using LAPACK call zgetri(d,inva,d,ipiv,work,lwork,info) ! again if info is not equal to 0, we exit setting status to 0 if (info /= 0) then status=0 end if deallocate(work) deallocate(ipiv) end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Determinants ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! function fvn_s_det(d,a,status) ! ! Evaluate the determinant of a square matrix using lapack LU factorization ! ! d (in) : matrix rank ! a (in) : The Matrix ! status (out) : =0 if LU factorization failed ! integer, intent(in) :: d real, intent(in) :: a(d,d) integer, intent(out) :: status real :: fvn_s_det real, allocatable :: wc_a(:,:) integer, allocatable :: ipiv(:) integer :: info,i status=1 allocate(wc_a(d,d)) allocate(ipiv(d)) wc_a(:,:)=a(:,:) call sgetrf(d,d,wc_a,d,ipiv,info) if (info/= 0) then status=0 fvn_s_det=0.e0 deallocate(ipiv) deallocate(wc_a) return end if fvn_s_det=1.e0 do i=1,d if (ipiv(i)==i) then fvn_s_det=fvn_s_det*wc_a(i,i) else fvn_s_det=-fvn_s_det*wc_a(i,i) end if end do deallocate(ipiv) deallocate(wc_a) end function function fvn_d_det(d,a,status) ! ! Evaluate the determinant of a square matrix using lapack LU factorization ! ! d (in) : matrix rank ! a (in) : The Matrix ! status (out) : =0 if LU factorization failed ! integer, intent(in) :: d double precision, intent(in) :: a(d,d) integer, intent(out) :: status double precision :: fvn_d_det double precision, allocatable :: wc_a(:,:) integer, allocatable :: ipiv(:) integer :: info,i status=1 allocate(wc_a(d,d)) allocate(ipiv(d)) wc_a(:,:)=a(:,:) call dgetrf(d,d,wc_a,d,ipiv,info) if (info/= 0) then status=0 fvn_d_det=0.d0 deallocate(ipiv) deallocate(wc_a) return end if fvn_d_det=1.d0 do i=1,d if (ipiv(i)==i) then fvn_d_det=fvn_d_det*wc_a(i,i) else fvn_d_det=-fvn_d_det*wc_a(i,i) end if end do deallocate(ipiv) deallocate(wc_a) end function function fvn_c_det(d,a,status) ! ! Evaluate the determinant of a square matrix using lapack LU factorization ! ! d (in) : matrix rank ! a (in) : The Matrix ! status (out) : =0 if LU factorization failed ! integer, intent(in) :: d complex, intent(in) :: a(d,d) integer, intent(out) :: status complex :: fvn_c_det complex, allocatable :: wc_a(:,:) integer, allocatable :: ipiv(:) integer :: info,i status=1 allocate(wc_a(d,d)) allocate(ipiv(d)) wc_a(:,:)=a(:,:) call cgetrf(d,d,wc_a,d,ipiv,info) if (info/= 0) then status=0 fvn_c_det=(0.e0,0.e0) deallocate(ipiv) deallocate(wc_a) return end if fvn_c_det=(1.e0,0.e0) do i=1,d if (ipiv(i)==i) then fvn_c_det=fvn_c_det*wc_a(i,i) else fvn_c_det=-fvn_c_det*wc_a(i,i) end if end do deallocate(ipiv) deallocate(wc_a) end function function fvn_z_det(d,a,status) ! ! Evaluate the determinant of a square matrix using lapack LU factorization ! ! d (in) : matrix rank ! a (in) : The Matrix ! det (out) : determinant ! status (out) : =0 if LU factorization failed ! integer, intent(in) :: d double complex, intent(in) :: a(d,d) integer, intent(out) :: status double complex :: fvn_z_det double complex, allocatable :: wc_a(:,:) integer, allocatable :: ipiv(:) integer :: info,i status=1 allocate(wc_a(d,d)) allocate(ipiv(d)) wc_a(:,:)=a(:,:) call zgetrf(d,d,wc_a,d,ipiv,info) if (info/= 0) then status=0 fvn_z_det=(0.d0,0.d0) deallocate(ipiv) deallocate(wc_a) return end if fvn_z_det=(1.d0,0.d0) do i=1,d if (ipiv(i)==i) then fvn_z_det=fvn_z_det*wc_a(i,i) else fvn_z_det=-fvn_z_det*wc_a(i,i) end if end do deallocate(ipiv) deallocate(wc_a) end function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Condition test ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! 1-norm ! fonction lapack slange,dlange,clange,zlange pour obtenir la 1-norm ! fonction lapack sgecon,dgecon,cgecon,zgecon pour calculer la rcond ! subroutine fvn_s_matcon(d,a,rcond,status) ! Matrix condition (reciprocal of condition number) ! ! d (in) : matrix rank ! a (in) : The Matrix ! rcond (out) : guess what ! status (out) : =0 if something went wrong ! integer, intent(in) :: d real, intent(in) :: a(d,d) real, intent(out) :: rcond integer, intent(out) :: status real, allocatable :: work(:) integer, allocatable :: iwork(:) real :: anorm real, allocatable :: wc_a(:,:) ! working copy of a integer :: info integer, allocatable :: ipiv(:) real, external :: slange status=1 anorm=slange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm allocate(wc_a(d,d)) !call scopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) allocate(ipiv(d)) call sgetrf(d,d,wc_a,d,ipiv,info) if (info /= 0) then status=0 deallocate(ipiv) deallocate(wc_a) return end if allocate(work(4*d)) allocate(iwork(d)) call sgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) if (info /= 0) then status=0 end if deallocate(iwork) deallocate(work) deallocate(ipiv) deallocate(wc_a) end subroutine subroutine fvn_d_matcon(d,a,rcond,status) ! Matrix condition (reciprocal of condition number) ! ! d (in) : matrix rank ! a (in) : The Matrix ! rcond (out) : guess what ! status (out) : =0 if something went wrong ! integer, intent(in) :: d double precision, intent(in) :: a(d,d) double precision, intent(out) :: rcond integer, intent(out) :: status double precision, allocatable :: work(:) integer, allocatable :: iwork(:) double precision :: anorm double precision, allocatable :: wc_a(:,:) ! working copy of a integer :: info integer, allocatable :: ipiv(:) double precision, external :: dlange status=1 anorm=dlange('1',d,d,a,d,work) ! work is unallocated as it is only used when computing infinity norm allocate(wc_a(d,d)) !call dcopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) allocate(ipiv(d)) call dgetrf(d,d,wc_a,d,ipiv,info) if (info /= 0) then status=0 deallocate(ipiv) deallocate(wc_a) return end if allocate(work(4*d)) allocate(iwork(d)) call dgecon('1',d,wc_a,d,anorm,rcond,work,iwork,info) if (info /= 0) then status=0 end if deallocate(iwork) deallocate(work) deallocate(ipiv) deallocate(wc_a) end subroutine subroutine fvn_c_matcon(d,a,rcond,status) ! Matrix condition (reciprocal of condition number) ! ! d (in) : matrix rank ! a (in) : The Matrix ! rcond (out) : guess what ! status (out) : =0 if something went wrong ! integer, intent(in) :: d complex, intent(in) :: a(d,d) real, intent(out) :: rcond integer, intent(out) :: status real, allocatable :: rwork(:) complex, allocatable :: work(:) integer, allocatable :: iwork(:) real :: anorm complex, allocatable :: wc_a(:,:) ! working copy of a integer :: info integer, allocatable :: ipiv(:) real, external :: clange status=1 anorm=clange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm allocate(wc_a(d,d)) !call ccopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) allocate(ipiv(d)) call cgetrf(d,d,wc_a,d,ipiv,info) if (info /= 0) then status=0 deallocate(ipiv) deallocate(wc_a) return end if allocate(work(2*d)) allocate(rwork(2*d)) call cgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) if (info /= 0) then status=0 end if deallocate(rwork) deallocate(work) deallocate(ipiv) deallocate(wc_a) end subroutine subroutine fvn_z_matcon(d,a,rcond,status) ! Matrix condition (reciprocal of condition number) ! ! d (in) : matrix rank ! a (in) : The Matrix ! rcond (out) : guess what ! status (out) : =0 if something went wrong ! integer, intent(in) :: d double complex, intent(in) :: a(d,d) double precision, intent(out) :: rcond integer, intent(out) :: status double complex, allocatable :: work(:) double precision, allocatable :: rwork(:) double precision :: anorm double complex, allocatable :: wc_a(:,:) ! working copy of a integer :: info integer, allocatable :: ipiv(:) double precision, external :: zlange status=1 anorm=zlange('1',d,d,a,d,rwork) ! rwork is unallocated as it is only used when computing infinity norm allocate(wc_a(d,d)) !call zcopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) allocate(ipiv(d)) call zgetrf(d,d,wc_a,d,ipiv,info) if (info /= 0) then status=0 deallocate(ipiv) deallocate(wc_a) return end if allocate(work(2*d)) allocate(rwork(2*d)) call zgecon('1',d,wc_a,d,anorm,rcond,work,rwork,info) if (info /= 0) then status=0 end if deallocate(rwork) deallocate(work) deallocate(ipiv) deallocate(wc_a) end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Valeurs propres/ Vecteurs propre ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine fvn_s_matev(d,a,evala,eveca,status) ! ! integer d (in) : matrice rank ! real a(d,d) (in) : The Matrix ! complex evala(d) (out) : eigenvalues ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector ! integer (out) : status =0 if something went wrong ! ! interfacing Lapack routine SGEEV integer, intent(in) :: d real, intent(in) :: a(d,d) complex, intent(out) :: evala(d) complex, intent(out) :: eveca(d,d) integer, intent(out) :: status real, allocatable :: wc_a(:,:) ! a working copy of a integer :: info integer :: lwork real, allocatable :: wr(:),wi(:) real :: vl ! unused but necessary for the call real, allocatable :: vr(:,:) real, allocatable :: work(:) real :: twork(1) integer i integer j ! making a working copy of a allocate(wc_a(d,d)) !call scopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) allocate(wr(d)) allocate(wi(d)) allocate(vr(d,d)) ! query optimal work size call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) lwork=int(twork(1)) allocate(work(lwork)) call sgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) if (info /= 0) then status=0 deallocate(work) deallocate(vr) deallocate(wi) deallocate(wr) deallocate(wc_a) return end if ! now fill in the results i=1 do while(i<=d) evala(i)=cmplx(wr(i),wi(i)) if (wi(i) == 0.) then ! eigenvalue is real eveca(:,i)=cmplx(vr(:,i),0.) else ! eigenvalue is complex evala(i+1)=cmplx(wr(i+1),wi(i+1)) eveca(:,i)=cmplx(vr(:,i),vr(:,i+1)) eveca(:,i+1)=cmplx(vr(:,i),-vr(:,i+1)) i=i+1 end if i=i+1 enddo deallocate(work) deallocate(vr) deallocate(wi) deallocate(wr) deallocate(wc_a) end subroutine subroutine fvn_d_matev(d,a,evala,eveca,status) ! ! integer d (in) : matrice rank ! double precision a(d,d) (in) : The Matrix ! double complex evala(d) (out) : eigenvalues ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector ! integer (out) : status =0 if something went wrong ! ! interfacing Lapack routine DGEEV integer, intent(in) :: d double precision, intent(in) :: a(d,d) double complex, intent(out) :: evala(d) double complex, intent(out) :: eveca(d,d) integer, intent(out) :: status double precision, allocatable :: wc_a(:,:) ! a working copy of a integer :: info integer :: lwork double precision, allocatable :: wr(:),wi(:) double precision :: vl ! unused but necessary for the call double precision, allocatable :: vr(:,:) double precision, allocatable :: work(:) double precision :: twork(1) integer i integer j ! making a working copy of a allocate(wc_a(d,d)) !call dcopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) allocate(wr(d)) allocate(wi(d)) allocate(vr(d,d)) ! query optimal work size call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,twork,-1,info) lwork=int(twork(1)) allocate(work(lwork)) call dgeev('N','V',d,wc_a,d,wr,wi,vl,1,vr,d,work,lwork,info) if (info /= 0) then status=0 deallocate(work) deallocate(vr) deallocate(wi) deallocate(wr) deallocate(wc_a) return end if ! now fill in the results i=1 do while(i<=d) evala(i)=dcmplx(wr(i),wi(i)) if (wi(i) == 0.) then ! eigenvalue is real eveca(:,i)=dcmplx(vr(:,i),0.) else ! eigenvalue is complex evala(i+1)=dcmplx(wr(i+1),wi(i+1)) eveca(:,i)=dcmplx(vr(:,i),vr(:,i+1)) eveca(:,i+1)=dcmplx(vr(:,i),-vr(:,i+1)) i=i+1 end if i=i+1 enddo deallocate(work) deallocate(vr) deallocate(wi) deallocate(wr) deallocate(wc_a) end subroutine subroutine fvn_c_matev(d,a,evala,eveca,status) ! ! integer d (in) : matrice rank ! complex a(d,d) (in) : The Matrix ! complex evala(d) (out) : eigenvalues ! complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector ! integer (out) : status =0 if something went wrong ! ! interfacing Lapack routine CGEEV integer, intent(in) :: d complex, intent(in) :: a(d,d) complex, intent(out) :: evala(d) complex, intent(out) :: eveca(d,d) integer, intent(out) :: status complex, allocatable :: wc_a(:,:) ! a working copy of a integer :: info integer :: lwork complex, allocatable :: work(:) complex :: twork(1) real, allocatable :: rwork(:) complex :: vl ! unused but necessary for the call status=1 ! making a working copy of a allocate(wc_a(d,d)) !call ccopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) ! query optimal work size call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) lwork=int(twork(1)) allocate(work(lwork)) allocate(rwork(2*d)) call cgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) if (info /= 0) then status=0 end if deallocate(rwork) deallocate(work) deallocate(wc_a) end subroutine subroutine fvn_z_matev(d,a,evala,eveca,status) ! ! integer d (in) : matrice rank ! double complex a(d,d) (in) : The Matrix ! double complex evala(d) (out) : eigenvalues ! double complex eveca(d,d) (out) : eveca(:,j) = jth eigenvector ! integer (out) : status =0 if something went wrong ! ! interfacing Lapack routine ZGEEV integer, intent(in) :: d double complex, intent(in) :: a(d,d) double complex, intent(out) :: evala(d) double complex, intent(out) :: eveca(d,d) integer, intent(out) :: status double complex, allocatable :: wc_a(:,:) ! a working copy of a integer :: info integer :: lwork double complex, allocatable :: work(:) double complex :: twork(1) double precision, allocatable :: rwork(:) double complex :: vl ! unused but necessary for the call status=1 ! making a working copy of a allocate(wc_a(d,d)) !call zcopy(d*d,a,1,wc_a,1) wc_a(:,:)=a(:,:) ! query optimal work size call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,twork,-1,rwork,info) lwork=int(twork(1)) allocate(work(lwork)) allocate(rwork(2*d)) call zgeev('N','V',d,wc_a,d,evala,vl,1,eveca,d,work,lwork,rwork,info) if (info /= 0) then status=0 end if deallocate(rwork) deallocate(work) deallocate(wc_a) end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Akima spline interpolation and spline evaluation ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Single precision subroutine fvn_s_akima(n,x,y,br,co) implicit none integer, intent(in) :: n real, intent(in) :: x(n) real, intent(in) :: y(n) real, intent(out) :: br(n) real, intent(out) :: co(4,n) real, allocatable :: var(:),z(:) real :: wi_1,wi integer :: i real :: dx,a,b ! br is just a copy of x br(:)=x(:) allocate(var(n)) allocate(z(n)) ! evaluate the variations do i=1, n-1 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) end do var(n+2)=2.e0*var(n+1)-var(n) var(n+3)=2.e0*var(n+2)-var(n+1) var(2)=2.e0*var(3)-var(4) var(1)=2.e0*var(2)-var(3) do i = 1, n wi_1=abs(var(i+3)-var(i+2)) wi=abs(var(i+1)-var(i)) if ((wi_1+wi).eq.0.e0) then z(i)=(var(i+2)+var(i+1))/2.e0 else z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) end if end do do i=1, n-1 dx=x(i+1)-x(i) 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 co(1,i)=y(i) co(2,i)=z(i) !co(3,i)=-(a-3.*b)/dx**2 ! méthode wd !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 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 ! etrangement la fonction csval corrige et donne la bonne valeur ... end do co(1,n)=y(n) co(2,n)=z(n) co(3,n)=0.e0 co(4,n)=0.e0 deallocate(z) deallocate(var) end subroutine ! Double precision subroutine fvn_d_akima(n,x,y,br,co) implicit none integer, intent(in) :: n double precision, intent(in) :: x(n) double precision, intent(in) :: y(n) double precision, intent(out) :: br(n) double precision, intent(out) :: co(4,n) double precision, allocatable :: var(:),z(:) double precision :: wi_1,wi integer :: i double precision :: dx,a,b ! br is just a copy of x br(:)=x(:) allocate(var(n)) allocate(z(n)) ! evaluate the variations do i=1, n-1 var(i+2)=(y(i+1)-y(i))/(x(i+1)-x(i)) end do var(n+2)=2.d0*var(n+1)-var(n) var(n+3)=2.d0*var(n+2)-var(n+1) var(2)=2.d0*var(3)-var(4) var(1)=2.d0*var(2)-var(3) do i = 1, n wi_1=dabs(var(i+3)-var(i+2)) wi=dabs(var(i+1)-var(i)) if ((wi_1+wi).eq.0.d0) then z(i)=(var(i+2)+var(i+1))/2.d0 else z(i)=(wi_1*var(i+1)+wi*var(i+2))/(wi_1+wi) end if end do do i=1, n-1 dx=x(i+1)-x(i) 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 co(1,i)=y(i) co(2,i)=z(i) !co(3,i)=-(a-3.*b)/dx**2 ! méthode wd !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 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 ! etrangement la fonction csval corrige et donne la bonne valeur ... end do co(1,n)=y(n) co(2,n)=z(n) co(3,n)=0.d0 co(4,n)=0.d0 deallocate(z) deallocate(var) end subroutine ! ! Single precision spline evaluation ! function fvn_s_spline_eval(x,n,br,co) implicit none real, intent(in) :: x ! x must be br(1)<= x <= br(n+1) otherwise value is extrapolated integer, intent(in) :: n ! number of intervals real, intent(in) :: br(n+1) ! breakpoints real, intent(in) :: co(4,n+1) ! spline coeeficients real :: fvn_s_spline_eval integer :: i real :: dx if (x<=br(1)) then i=1 else if (x>=br(n+1)) then i=n else i=1 do while(x>=br(i)) i=i+1 end do i=i-1 end if 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 end function ! Double precision spline evaluation function fvn_d_spline_eval(x,n,br,co) implicit none 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 double precision, intent(in) :: br(n+1) ! breakpoints double precision, intent(in) :: co(4,n+1) ! spline coeeficients double precision :: fvn_d_spline_eval integer :: i double precision :: dx if (x<=br(1)) then i=1 else if (x>=br(n+1)) then i=n else i=1 do while(x>=br(i)) i=i+1 end do i=i-1 end if 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 end function ! ! Muller ! ! ! ! William Daniau 2007 ! ! This routine is a fortran 90 port of Hans D. Mittelmann's routine muller.f ! http://plato.asu.edu/ftp/other_software/muller.f ! ! it can be used as a replacement for imsl routine dzanly with minor changes ! !----------------------------------------------------------------------- ! ! purpose - zeros of an analytic complex function ! using the muller method with deflation ! ! usage - call fvn_z_muller (f,eps,eps1,kn,n,nguess,x,itmax, ! infer,ier) ! ! arguments f - a complex function subprogram, f(z), written ! by the user specifying the equation whose ! roots are to be found. f must appear in ! an external statement in the calling pro- ! gram. ! eps - 1st stopping criterion. let fp(z)=f(z)/p ! where p = (z-z(1))*(z-z(2))*,,,*(z-z(k-1)) ! and z(1),...,z(k-1) are previously found ! roots. if ((cdabs(f(z)).le.eps) .and. ! (cdabs(fp(z)).le.eps)), then z is accepted ! as a root. (input) ! eps1 - 2nd stopping criterion. a root is accepted ! if two successive approximations to a given ! root agree within eps1. (input) ! note. if either or both of the stopping ! criteria are fulfilled, the root is ! accepted. ! kn - the number of known roots which must be stored ! in x(1),...,x(kn), prior to entry to muller ! nguess - the number of initial guesses provided. these ! guesses must be stored in x(kn+1),..., ! x(kn+nguess). nguess must be set equal ! to zero if no guesses are provided. (input) ! n - the number of new roots to be found by ! muller (input) ! x - a complex vector of length kn+n. x(1),..., ! x(kn) on input must contain any known ! roots. x(kn+1),..., x(kn+n) on input may, ! on user option, contain initial guesses for ! the n new roots which are to be computed. ! if the user does not provide an initial ! guess, zero is used. ! on output, x(kn+1),...,x(kn+n) contain the ! approximate roots found by muller. ! itmax - the maximum allowable number of iterations ! per root (input) ! infer - an integer vector of length kn+n. on ! output infer(j) contains the number of ! iterations used in finding the j-th root ! when convergence was achieved. if ! convergence was not obtained in itmax ! iterations, infer(j) will be greater than ! itmax (output). ! ier - error parameter (output) ! warning error ! ier = 33 indicates failure to converge with- ! in itmax iterations for at least one of ! the (n) new roots. ! ! ! remarks muller always returns the last approximation for root j ! in x(j). if the convergence criterion is satisfied, ! then infer(j) is less than or equal to itmax. if the ! convergence criterion is not satisified, then infer(j) ! is set to either itmax+1 or itmax+k, with k greater ! than 1. infer(j) = itmax+1 indicates that muller did ! not obtain convergence in the allowed number of iter- ! ations. in this case, the user may wish to set itmax ! to a larger value. infer(j) = itmax+k means that con- ! vergence was obtained (on iteration k) for the defla- ! ted function ! fp(z) = f(z)/((z-z(1)...(z-z(j-1))) ! ! but failed for f(z). in this case, better initial ! guesses might help or, it might be necessary to relax ! the convergence criterion. ! !----------------------------------------------------------------------- ! subroutine fvn_z_muller (f,eps,eps1,kn,nguess,n,x,itmax,infer,ier) implicit none double precision :: rzero,rten,rhun,rp01,ax,eps1,qz,eps,tpq double complex :: d,dd,den,fprt,frt,h,rt,t1,t2,t3, & tem,z0,z1,z2,bi,xx,xl,y0,y1,y2,x0, & zero,p1,one,four,p5 double complex, external :: f integer :: ickmax,kn,nguess,n,itmax,ier,knp1,knpn,i,l,ic, & knpng,jk,ick,nn,lm1,errcode double complex :: x(kn+n) integer :: infer(kn+n) data zero/(0.0,0.0)/,p1/(0.1,0.0)/, & one/(1.0,0.0)/,four/(4.0,0.0)/, & p5/(0.5,0.0)/, & rzero/0.0/,rten/10.0/,rhun/100.0/, & ax/0.1/,ickmax/3/,rp01/0.01/ ier = 0 if (n .lt. 1) then ! What the hell are doing here then ... return end if !eps1 = rten **(-nsig) eps1 = min(eps1,rp01) knp1 = kn+1 knpn = kn+n knpng = kn+nguess do i=1,knpn infer(i) = 0 if (i .gt. knpng) x(i) = zero end do l= knp1 ic=0 rloop: do while (l<=knpn) ! Main loop over new roots jk = 0 ick = 0 xl = x(l) icloop: do ic = 0 h = ax h = p1*h if (cdabs(xl) .gt. ax) h = p1*xl ! first three points are ! xl+h, xl-h, xl rt = xl+h call deflated_work(errcode) if (errcode == 1) then exit icloop end if z0 = fprt y0 = frt x0 = rt rt = xl-h call deflated_work(errcode) if (errcode == 1) then exit icloop end if z1 = fprt y1 = frt h = xl-rt d = h/(rt-x0) rt = xl call deflated_work(errcode) if (errcode == 1) then exit icloop end if z2 = fprt y2 = frt ! begin main algorithm iloop: do dd = one + d t1 = z0*d*d t2 = z1*dd*dd xx = z2*dd t3 = z2*d bi = t1-t2+xx+t3 den = bi*bi-four*(xx*t1-t3*(t2-xx)) ! use denominator of maximum amplitude t1 = cdsqrt(den) qz = rhun*max(cdabs(bi),cdabs(t1)) t2 = bi + t1 tpq = cdabs(t2)+qz if (tpq .eq. qz) t2 = zero t3 = bi - t1 tpq = cdabs(t3) + qz if (tpq .eq. qz) t3 = zero den = t2 qz = cdabs(t3)-cdabs(t2) if (qz .gt. rzero) den = t3 ! test for zero denominator if (cdabs(den) .eq. rzero) then call trans_rt() call deflated_work(errcode) if (errcode == 1) then exit icloop end if z2 = fprt y2 = frt cycle iloop end if d = -xx/den d = d+d h = d*h rt = rt + h ! check convergence of the first kind if (cdabs(h) .le. eps1*max(cdabs(rt),ax)) then if (ic .ne. 0) then exit icloop end if ic = 1 z0 = y1 z1 = y2 z2 = f(rt) xl = rt ick = ick+1 if (ick .le. ickmax) then cycle iloop end if ! warning error, itmax = maximum jk = itmax + jk ier = 33 end if if (ic .ne. 0) then cycle icloop end if call deflated_work(errcode) if (errcode == 1) then exit icloop end if do while ( (cdabs(fprt)-cdabs(z2)*rten) .ge. rzero) ! take remedial action to induce ! convergence d = d*p5 h = h*p5 rt = rt-h call deflated_work(errcode) if (errcode == 1) then exit icloop end if end do z0 = z1 z1 = z2 z2 = fprt y0 = y1 y1 = y2 y2 = frt end do iloop end do icloop x(l) = rt infer(l) = jk l = l+1 end do rloop contains subroutine trans_rt() tem = rten*eps1 if (cdabs(rt) .gt. ax) tem = tem*rt rt = rt+tem d = (h+tem)*d/h h = h+tem end subroutine trans_rt subroutine deflated_work(errcode) ! errcode=0 => no errors ! errcode=1 => jk>itmax or convergence of second kind achieved integer :: errcode,flag flag=1 loop1: do while(flag==1) errcode=0 jk = jk+1 if (jk .gt. itmax) then ier=33 errcode=1 return end if frt = f(rt) fprt = frt if (l /= 1) then lm1 = l-1 do i=1,lm1 tem = rt - x(i) if (cdabs(tem) .eq. rzero) then !if (ic .ne. 0) go to 15 !! ?? possible? call trans_rt() cycle loop1 end if fprt = fprt/tem end do end if flag=0 end do loop1 if (cdabs(fprt) .le. eps .and. cdabs(frt) .le. eps) then errcode=1 return end if end subroutine deflated_work end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! ! Integration ! ! Only double precision coded atm ! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine fvn_d_gauss_legendre(n,qx,qw) ! ! This routine compute the n Gauss Legendre abscissas and weights ! Adapted from Numerical Recipes routine gauleg ! ! n (in) : number of points ! qx(out) : abscissas ! qw(out) : weights ! implicit none double precision,parameter :: pi=3.141592653589793d0 integer, intent(in) :: n double precision, intent(out) :: qx(n),qw(n) integer :: m,i,j double precision :: z,z1,p1,p2,p3,pp m=(n+1)/2 do i=1,m z=cos(pi*(dble(i)-0.25d0)/(dble(n)+0.5d0)) iloop: do p1=1.d0 p2=0.d0 do j=1,n p3=p2 p2=p1 p1=((2.d0*dble(j)-1.d0)*z*p2-(dble(j)-1.d0)*p3)/dble(j) end do pp=dble(n)*(z*p1-p2)/(z*z-1.d0) z1=z z=z1-p1/pp if (dabs(z-z1)<=epsilon(z)) then exit iloop end if end do iloop qx(i)=-z qx(n+1-i)=z qw(i)=2.d0/((1.d0-z*z)*pp*pp) qw(n+1-i)=qw(i) end do end subroutine subroutine fvn_d_gl_integ(f,a,b,n,res) ! ! This is a simple non adaptative integration routine ! using n gauss legendre abscissas and weights ! ! f(in) : the function to integrate ! a(in) : lower bound ! b(in) : higher bound ! n(in) : number of gauss legendre pairs ! res(out): the evaluation of the integral ! double precision,external :: f double precision, intent(in) :: a,b integer, intent(in):: n double precision, intent(out) :: res double precision, allocatable :: qx(:),qw(:) double precision :: xm,xr integer :: i ! First compute n gauss legendre abs and weight allocate(qx(n)) allocate(qw(n)) call fvn_d_gauss_legendre(n,qx,qw) xm=0.5d0*(b+a) xr=0.5d0*(b-a) res=0.d0 do i=1,n res=res+qw(i)*f(xm+xr*qx(i)) end do res=xr*res deallocate(qw) deallocate(qx) end subroutine !!!!!!!!!!!!!!!!!!!!!!!! ! ! Simple and double adaptative Gauss Kronrod integration based on ! a modified version of quadpack ( http://www.netlib.org/quadpack ! ! Common parameters : ! ! key (in) ! epsabs ! epsrel ! ! !!!!!!!!!!!!!!!!!!!!!!!! subroutine fvn_d_integ_1_gk(f,a,b,epsabs,epsrel,key,res,abserr,ier,limit) ! ! Evaluate the integral of function f(x) between a and b ! ! f(in) : the function ! a(in) : lower bound ! b(in) : higher bound ! epsabs(in) : desired absolute error ! epsrel(in) : desired relative error ! key(in) : gauss kronrod rule ! 1: 7 - 15 points ! 2: 10 - 21 points ! 3: 15 - 31 points ! 4: 20 - 41 points ! 5: 25 - 51 points ! 6: 30 - 61 points ! ! limit(in) : maximum number of subintervals in the partition of the ! given integration interval (a,b). A value of 500 will give the same ! behaviour as the imsl routine dqdag ! ! res(out) : estimated integral value ! abserr(out) : estimated absolute error ! ier(out) : error flag from quadpack routines ! 0 : no error ! 1 : maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yield no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulaties. ! if the position of a local difficulty can ! be determined (i.e.singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used which is designed for ! handling the type of difficulty involved. ! 2 : the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! 3 : extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! 6 : the input is invalid, because ! (epsabs.le.0 and ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) ! or limit.lt.1 or lenw.lt.limit*4. ! result, abserr, neval, last are set ! to zero. ! except when lenw is invalid, iwork(1), ! work(limit*2+1) and work(limit*3+1) are ! set to zero, work(1) is set to a and ! work(limit+1) to b. implicit none double precision, external :: f double precision, intent(in) :: a,b,epsabs,epsrel integer, intent(in) :: key integer, intent(in) :: limit double precision, intent(out) :: res,abserr integer, intent(out) :: ier double precision, allocatable :: work(:) integer, allocatable :: iwork(:) integer :: lenw,neval,last ! imsl value for limit is 500 lenw=limit*4 allocate(iwork(limit)) allocate(work(lenw)) call dqag(f,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) deallocate(work) deallocate(iwork) end subroutine subroutine fvn_d_integ_2_gk(f,a,b,g,h,epsabs,epsrel,key,res,abserr,ier,limit) ! ! Evaluate the double integral of function f(x,y) for x between a and b and y between g(x) and h(x) ! ! f(in) : the function ! a(in) : lower bound ! b(in) : higher bound ! g(in) : external function describing lower bound for y ! h(in) : external function describing higher bound for y ! epsabs(in) : desired absolute error ! epsrel(in) : desired relative error ! key(in) : gauss kronrod rule ! 1: 7 - 15 points ! 2: 10 - 21 points ! 3: 15 - 31 points ! 4: 20 - 41 points ! 5: 25 - 51 points ! 6: 30 - 61 points ! ! limit(in) : maximum number of subintervals in the partition of the ! given integration interval (a,b). A value of 500 will give the same ! behaviour as the imsl routine dqdag ! ! res(out) : estimated integral value ! abserr(out) : estimated absolute error ! ier(out) : error flag from quadpack routines ! 0 : no error ! 1 : maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yield no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulaties. ! if the position of a local difficulty can ! be determined (i.e.singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used which is designed for ! handling the type of difficulty involved. ! 2 : the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! 3 : extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! 6 : the input is invalid, because ! (epsabs.le.0 and ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) ! or limit.lt.1 or lenw.lt.limit*4. ! result, abserr, neval, last are set ! to zero. ! except when lenw is invalid, iwork(1), ! work(limit*2+1) and work(limit*3+1) are ! set to zero, work(1) is set to a and ! work(limit+1) to b. implicit none double precision, external:: f,g,h double precision, intent(in) :: a,b,epsabs,epsrel integer, intent(in) :: key,limit integer, intent(out) :: ier double precision, intent(out) :: res,abserr double precision, allocatable :: work(:) integer, allocatable :: iwork(:) integer :: lenw,neval,last ! imsl value for limit is 500 lenw=limit*4 allocate(work(lenw)) allocate(iwork(limit)) call dqag_2d_outer(f,a,b,g,h,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) deallocate(iwork) deallocate(work) end subroutine subroutine fvn_d_integ_2_inner_gk(f,x,a,b,epsabs,epsrel,key,res,abserr,ier,limit) ! ! Evaluate the single integral of function f(x,y) for y between a and b with a ! given x value ! ! This function is used for the evaluation of the double integral fvn_d_integ_2_gk ! ! f(in) : the function ! x(in) : x ! a(in) : lower bound ! b(in) : higher bound ! epsabs(in) : desired absolute error ! epsrel(in) : desired relative error ! key(in) : gauss kronrod rule ! 1: 7 - 15 points ! 2: 10 - 21 points ! 3: 15 - 31 points ! 4: 20 - 41 points ! 5: 25 - 51 points ! 6: 30 - 61 points ! ! limit(in) : maximum number of subintervals in the partition of the ! given integration interval (a,b). A value of 500 will give the same ! behaviour as the imsl routine dqdag ! ! res(out) : estimated integral value ! abserr(out) : estimated absolute error ! ier(out) : error flag from quadpack routines ! 0 : no error ! 1 : maximum number of subdivisions allowed ! has been achieved. one can allow more ! subdivisions by increasing the value of ! limit (and taking the according dimension ! adjustments into account). however, if ! this yield no improvement it is advised ! to analyze the integrand in order to ! determine the integration difficulaties. ! if the position of a local difficulty can ! be determined (i.e.singularity, ! discontinuity within the interval) one ! will probably gain from splitting up the ! interval at this point and calling the ! integrator on the subranges. if possible, ! an appropriate special-purpose integrator ! should be used which is designed for ! handling the type of difficulty involved. ! 2 : the occurrence of roundoff error is ! detected, which prevents the requested ! tolerance from being achieved. ! 3 : extremely bad integrand behaviour occurs ! at some points of the integration ! interval. ! 6 : the input is invalid, because ! (epsabs.le.0 and ! epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) ! or limit.lt.1 or lenw.lt.limit*4. ! result, abserr, neval, last are set ! to zero. ! except when lenw is invalid, iwork(1), ! work(limit*2+1) and work(limit*3+1) are ! set to zero, work(1) is set to a and ! work(limit+1) to b. implicit none double precision, external:: f double precision, intent(in) :: x,a,b,epsabs,epsrel integer, intent(in) :: key,limit integer, intent(out) :: ier double precision, intent(out) :: res,abserr double precision, allocatable :: work(:) integer, allocatable :: iwork(:) integer :: lenw,neval,last ! imsl value for limit is 500 lenw=limit*4 allocate(work(lenw)) allocate(iwork(limit)) call dqag_2d_inner(f,x,a,b,epsabs,epsrel,key,res,abserr,neval,ier,limit,lenw,last,iwork,work) deallocate(iwork) deallocate(work) end subroutine !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! Include the modified quadpack files !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! include "fvn_quadpack/dqag_2d_inner.f" include "fvn_quadpack/dqk15_2d_inner.f" include "fvn_quadpack/dqk31_2d_outer.f" include "fvn_quadpack/d1mach.f" include "fvn_quadpack/dqk31_2d_inner.f" include "fvn_quadpack/dqage.f" include "fvn_quadpack/dqk15.f" include "fvn_quadpack/dqk21.f" include "fvn_quadpack/dqk31.f" include "fvn_quadpack/dqk41.f" include "fvn_quadpack/dqk51.f" include "fvn_quadpack/dqk61.f" include "fvn_quadpack/dqk41_2d_outer.f" include "fvn_quadpack/dqk41_2d_inner.f" include "fvn_quadpack/dqag_2d_outer.f" include "fvn_quadpack/dqpsrt.f" include "fvn_quadpack/dqag.f" include "fvn_quadpack/dqage_2d_outer.f" include "fvn_quadpack/dqage_2d_inner.f" include "fvn_quadpack/dqk51_2d_outer.f" include "fvn_quadpack/dqk51_2d_inner.f" include "fvn_quadpack/dqk61_2d_outer.f" include "fvn_quadpack/dqk21_2d_outer.f" include "fvn_quadpack/dqk61_2d_inner.f" include "fvn_quadpack/dqk21_2d_inner.f" include "fvn_quadpack/dqk15_2d_outer.f" end module fvn