diff --git a/stable/fvnlib.f90 b/stable/fvnlib.f90 deleted file mode 100644 index c4eceac..0000000 --- a/stable/fvnlib.f90 +++ /dev/null @@ -1,1779 +0,0 @@ - -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 ;-) Nevertheless, you -! may give credits to quadpack authors. -! -! 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 \ No newline at end of file