! Friday, September 12, 2003 at 9:28 ! ! calcule et estime les moments d'un ARMA(p,q)-GARCH(pp,qq) et les ecarts-types asymptotiques ! en utilisant la formule de Bartlett generalisee ! ! MODULE nrtype INTEGER, PARAMETER :: I4B = SELECTED_INT_KIND(9) INTEGER, PARAMETER :: I2B = SELECTED_INT_KIND(4) INTEGER, PARAMETER :: I1B = SELECTED_INT_KIND(2) INTEGER, PARAMETER :: SP = KIND(1.0) INTEGER, PARAMETER :: DP = KIND(1.0D0) INTEGER, PARAMETER :: SPC = KIND((1.0,1.0)) INTEGER, PARAMETER :: DPC = KIND((1.0D0,1.0D0)) INTEGER, PARAMETER :: LGT = KIND(.true.) REAL(SP), PARAMETER :: PI=3.141592653589793238462643383279502884197_sp REAL(SP), PARAMETER :: PIO2=1.57079632679489661923132169163975144209858_sp REAL(SP), PARAMETER :: TWOPI=6.283185307179586476925286766559005768394_sp REAL(SP), PARAMETER :: SQRT2=1.41421356237309504880168872420969807856967_sp REAL(SP), PARAMETER :: EULER=0.5772156649015328606065120900824024310422_sp REAL(DP), PARAMETER :: PI_D=3.141592653589793238462643383279502884197_dp REAL(DP), PARAMETER :: PIO2_D=1.57079632679489661923132169163975144209858_dp REAL(DP), PARAMETER :: TWOPI_D=6.283185307179586476925286766559005768394_dp TYPE sprs2_sp INTEGER(I4B) :: n,len REAL(SP), DIMENSION(:), POINTER :: val INTEGER(I4B), DIMENSION(:), POINTER :: irow INTEGER(I4B), DIMENSION(:), POINTER :: jcol END TYPE sprs2_sp TYPE sprs2_dp INTEGER(I4B) :: n,len REAL(DP), DIMENSION(:), POINTER :: val INTEGER(I4B), DIMENSION(:), POINTER :: irow INTEGER(I4B), DIMENSION(:), POINTER :: jcol END TYPE sprs2_dp END MODULE nrtype MODULE nrutil USE nrtype IMPLICIT NONE INTEGER(I4B), PARAMETER :: NPAR_ARTH=16,NPAR2_ARTH=8 INTEGER(I4B), PARAMETER :: NPAR_GEOP=4,NPAR2_GEOP=2 INTEGER(I4B), PARAMETER :: NPAR_CUMSUM=16 INTEGER(I4B), PARAMETER :: NPAR_CUMPROD=8 INTEGER(I4B), PARAMETER :: NPAR_POLY=8 INTEGER(I4B), PARAMETER :: NPAR_POLYTERM=8 INTERFACE array_copy MODULE PROCEDURE array_copy_r, array_copy_d, array_copy_i END INTERFACE INTERFACE swap MODULE PROCEDURE swap_i,swap_r,swap_rv,swap_c, & swap_cv,swap_cm,swap_z,swap_zv,swap_zm, & masked_swap_rs,masked_swap_rv,masked_swap_rm END INTERFACE INTERFACE reallocate MODULE PROCEDURE reallocate_rv,reallocate_rm,& reallocate_iv,reallocate_im,reallocate_hv END INTERFACE INTERFACE imaxloc MODULE PROCEDURE imaxloc_r,imaxloc_i END INTERFACE INTERFACE assert MODULE PROCEDURE assert1,assert2,assert3,assert4,assert_v END INTERFACE INTERFACE assert_eq MODULE PROCEDURE assert_eq2,assert_eq3,assert_eq4,assert_eqn END INTERFACE INTERFACE arth MODULE PROCEDURE arth_r, arth_d, arth_i END INTERFACE INTERFACE geop MODULE PROCEDURE geop_r, geop_d, geop_i, geop_c, geop_dv END INTERFACE INTERFACE cumsum MODULE PROCEDURE cumsum_r,cumsum_i END INTERFACE INTERFACE poly MODULE PROCEDURE poly_rr,poly_rrv,poly_dd,poly_ddv,& poly_rc,poly_cc,poly_msk_rrv,poly_msk_ddv END INTERFACE INTERFACE poly_term MODULE PROCEDURE poly_term_rr,poly_term_cc END INTERFACE INTERFACE outerprod MODULE PROCEDURE outerprod_r,outerprod_d END INTERFACE INTERFACE outerdiff MODULE PROCEDURE outerdiff_r,outerdiff_d,outerdiff_i END INTERFACE INTERFACE scatter_add MODULE PROCEDURE scatter_add_r,scatter_add_d END INTERFACE INTERFACE scatter_max MODULE PROCEDURE scatter_max_r,scatter_max_d END INTERFACE INTERFACE diagadd MODULE PROCEDURE diagadd_rv,diagadd_r END INTERFACE INTERFACE diagmult MODULE PROCEDURE diagmult_rv,diagmult_r END INTERFACE INTERFACE get_diag MODULE PROCEDURE get_diag_rv, get_diag_dv END INTERFACE INTERFACE put_diag MODULE PROCEDURE put_diag_rv, put_diag_r END INTERFACE CONTAINS !BL SUBROUTINE array_copy_r(src,dest,n_copied,n_not_copied) REAL(SP), DIMENSION(:), INTENT(IN) :: src REAL(SP), DIMENSION(:), INTENT(OUT) :: dest INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied n_copied=min(size(src),size(dest)) n_not_copied=size(src)-n_copied dest(1:n_copied)=src(1:n_copied) END SUBROUTINE array_copy_r !BL SUBROUTINE array_copy_d(src,dest,n_copied,n_not_copied) REAL(DP), DIMENSION(:), INTENT(IN) :: src REAL(DP), DIMENSION(:), INTENT(OUT) :: dest INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied n_copied=min(size(src),size(dest)) n_not_copied=size(src)-n_copied dest(1:n_copied)=src(1:n_copied) END SUBROUTINE array_copy_d !BL SUBROUTINE array_copy_i(src,dest,n_copied,n_not_copied) INTEGER(I4B), DIMENSION(:), INTENT(IN) :: src INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: dest INTEGER(I4B), INTENT(OUT) :: n_copied, n_not_copied n_copied=min(size(src),size(dest)) n_not_copied=size(src)-n_copied dest(1:n_copied)=src(1:n_copied) END SUBROUTINE array_copy_i !BL !BL SUBROUTINE swap_i(a,b) INTEGER(I4B), INTENT(INOUT) :: a,b INTEGER(I4B) :: dum dum=a a=b b=dum END SUBROUTINE swap_i !BL SUBROUTINE swap_r(a,b) REAL(SP), INTENT(INOUT) :: a,b REAL(SP) :: dum dum=a a=b b=dum END SUBROUTINE swap_r !BL SUBROUTINE swap_rv(a,b) REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b REAL(SP), DIMENSION(SIZE(a)) :: dum dum=a a=b b=dum END SUBROUTINE swap_rv !BL SUBROUTINE swap_c(a,b) COMPLEX(SPC), INTENT(INOUT) :: a,b COMPLEX(SPC) :: dum dum=a a=b b=dum END SUBROUTINE swap_c !BL SUBROUTINE swap_cv(a,b) COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: a,b COMPLEX(SPC), DIMENSION(SIZE(a)) :: dum dum=a a=b b=dum END SUBROUTINE swap_cv !BL SUBROUTINE swap_cm(a,b) COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: a,b COMPLEX(SPC), DIMENSION(size(a,1),size(a,2)) :: dum dum=a a=b b=dum END SUBROUTINE swap_cm !BL SUBROUTINE swap_z(a,b) COMPLEX(DPC), INTENT(INOUT) :: a,b COMPLEX(DPC) :: dum dum=a a=b b=dum END SUBROUTINE swap_z !BL SUBROUTINE swap_zv(a,b) COMPLEX(DPC), DIMENSION(:), INTENT(INOUT) :: a,b COMPLEX(DPC), DIMENSION(SIZE(a)) :: dum dum=a a=b b=dum END SUBROUTINE swap_zv !BL SUBROUTINE swap_zm(a,b) COMPLEX(DPC), DIMENSION(:,:), INTENT(INOUT) :: a,b COMPLEX(DPC), DIMENSION(size(a,1),size(a,2)) :: dum dum=a a=b b=dum END SUBROUTINE swap_zm !BL SUBROUTINE masked_swap_rs(a,b,mask) REAL(SP), INTENT(INOUT) :: a,b LOGICAL(LGT), INTENT(IN) :: mask REAL(SP) :: swp if (mask) then swp=a a=b b=swp end if END SUBROUTINE masked_swap_rs !BL SUBROUTINE masked_swap_rv(a,b,mask) REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask REAL(SP), DIMENSION(size(a)) :: swp where (mask) swp=a a=b b=swp end where END SUBROUTINE masked_swap_rv !BL SUBROUTINE masked_swap_rm(a,b,mask) REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b LOGICAL(LGT), DIMENSION(:,:), INTENT(IN) :: mask REAL(SP), DIMENSION(size(a,1),size(a,2)) :: swp where (mask) swp=a a=b b=swp end where END SUBROUTINE masked_swap_rm !BL !BL FUNCTION reallocate_rv(p,n) REAL(SP), DIMENSION(:), POINTER :: p, reallocate_rv INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B) :: nold,ierr allocate(reallocate_rv(n),stat=ierr) if (ierr /= 0) call & nrerror('reallocate_rv: problem in attempt to allocate memory') if (.not. associated(p)) RETURN nold=size(p) reallocate_rv(1:min(nold,n))=p(1:min(nold,n)) deallocate(p) END FUNCTION reallocate_rv !BL FUNCTION reallocate_iv(p,n) INTEGER(I4B), DIMENSION(:), POINTER :: p, reallocate_iv INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B) :: nold,ierr allocate(reallocate_iv(n),stat=ierr) if (ierr /= 0) call & nrerror('reallocate_iv: problem in attempt to allocate memory') if (.not. associated(p)) RETURN nold=size(p) reallocate_iv(1:min(nold,n))=p(1:min(nold,n)) deallocate(p) END FUNCTION reallocate_iv !BL FUNCTION reallocate_hv(p,n) CHARACTER(1), DIMENSION(:), POINTER :: p, reallocate_hv INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B) :: nold,ierr allocate(reallocate_hv(n),stat=ierr) if (ierr /= 0) call & nrerror('reallocate_hv: problem in attempt to allocate memory') if (.not. associated(p)) RETURN nold=size(p) reallocate_hv(1:min(nold,n))=p(1:min(nold,n)) deallocate(p) END FUNCTION reallocate_hv !BL FUNCTION reallocate_rm(p,n,m) REAL(SP), DIMENSION(:,:), POINTER :: p, reallocate_rm INTEGER(I4B), INTENT(IN) :: n,m INTEGER(I4B) :: nold,mold,ierr allocate(reallocate_rm(n,m),stat=ierr) if (ierr /= 0) call & nrerror('reallocate_rm: problem in attempt to allocate memory') if (.not. associated(p)) RETURN nold=size(p,1) mold=size(p,2) reallocate_rm(1:min(nold,n),1:min(mold,m))=& p(1:min(nold,n),1:min(mold,m)) deallocate(p) END FUNCTION reallocate_rm !BL FUNCTION reallocate_im(p,n,m) INTEGER(I4B), DIMENSION(:,:), POINTER :: p, reallocate_im INTEGER(I4B), INTENT(IN) :: n,m INTEGER(I4B) :: nold,mold,ierr allocate(reallocate_im(n,m),stat=ierr) if (ierr /= 0) call & nrerror('reallocate_im: problem in attempt to allocate memory') if (.not. associated(p)) RETURN nold=size(p,1) mold=size(p,2) reallocate_im(1:min(nold,n),1:min(mold,m))=& p(1:min(nold,n),1:min(mold,m)) deallocate(p) END FUNCTION reallocate_im !BL FUNCTION ifirstloc(mask) LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask INTEGER(I4B) :: ifirstloc INTEGER(I4B), DIMENSION(1) :: loc loc=maxloc(merge(1,0,mask)) ifirstloc=loc(1) if (.not. mask(ifirstloc)) ifirstloc=size(mask)+1 END FUNCTION ifirstloc !BL FUNCTION imaxloc_r(arr) REAL(SP), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B) :: imaxloc_r INTEGER(I4B), DIMENSION(1) :: imax imax=maxloc(arr(:)) imaxloc_r=imax(1) END FUNCTION imaxloc_r !BL FUNCTION imaxloc_i(iarr) INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr INTEGER(I4B), DIMENSION(1) :: imax INTEGER(I4B) :: imaxloc_i imax=maxloc(iarr(:)) imaxloc_i=imax(1) END FUNCTION imaxloc_i !BL FUNCTION iminloc(arr) REAL(SP), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), DIMENSION(1) :: imin INTEGER(I4B) :: iminloc imin=minloc(arr(:)) iminloc=imin(1) END FUNCTION iminloc !BL SUBROUTINE assert1(n1,string) CHARACTER(LEN=*), INTENT(IN) :: string LOGICAL, INTENT(IN) :: n1 if (.not. n1) then write (*,*) 'nrerror: an assertion failed with this tag:', & string STOP 'program terminated by assert1' end if END SUBROUTINE assert1 !BL SUBROUTINE assert2(n1,n2,string) CHARACTER(LEN=*), INTENT(IN) :: string LOGICAL, INTENT(IN) :: n1,n2 if (.not. (n1 .and. n2)) then write (*,*) 'nrerror: an assertion failed with this tag:', & string STOP 'program terminated by assert2' end if END SUBROUTINE assert2 !BL SUBROUTINE assert3(n1,n2,n3,string) CHARACTER(LEN=*), INTENT(IN) :: string LOGICAL, INTENT(IN) :: n1,n2,n3 if (.not. (n1 .and. n2 .and. n3)) then write (*,*) 'nrerror: an assertion failed with this tag:', & string STOP 'program terminated by assert3' end if END SUBROUTINE assert3 !BL SUBROUTINE assert4(n1,n2,n3,n4,string) CHARACTER(LEN=*), INTENT(IN) :: string LOGICAL, INTENT(IN) :: n1,n2,n3,n4 if (.not. (n1 .and. n2 .and. n3 .and. n4)) then write (*,*) 'nrerror: an assertion failed with this tag:', & string STOP 'program terminated by assert4' end if END SUBROUTINE assert4 !BL SUBROUTINE assert_v(n,string) CHARACTER(LEN=*), INTENT(IN) :: string LOGICAL, DIMENSION(:), INTENT(IN) :: n if (.not. all(n)) then write (*,*) 'nrerror: an assertion failed with this tag:', & string STOP 'program terminated by assert_v' end if END SUBROUTINE assert_v !BL FUNCTION assert_eq2(n1,n2,string) CHARACTER(LEN=*), INTENT(IN) :: string INTEGER, INTENT(IN) :: n1,n2 INTEGER :: assert_eq2 if (n1 == n2) then assert_eq2=n1 else write (*,*) 'nrerror: an assert_eq failed with this tag:', & string STOP 'program terminated by assert_eq2' end if END FUNCTION assert_eq2 !BL FUNCTION assert_eq3(n1,n2,n3,string) CHARACTER(LEN=*), INTENT(IN) :: string INTEGER, INTENT(IN) :: n1,n2,n3 INTEGER :: assert_eq3 if (n1 == n2 .and. n2 == n3) then assert_eq3=n1 else write (*,*) 'nrerror: an assert_eq failed with this tag:', & string STOP 'program terminated by assert_eq3' end if END FUNCTION assert_eq3 !BL FUNCTION assert_eq4(n1,n2,n3,n4,string) CHARACTER(LEN=*), INTENT(IN) :: string INTEGER, INTENT(IN) :: n1,n2,n3,n4 INTEGER :: assert_eq4 if (n1 == n2 .and. n2 == n3 .and. n3 == n4) then assert_eq4=n1 else write (*,*) 'nrerror: an assert_eq failed with this tag:', & string STOP 'program terminated by assert_eq4' end if END FUNCTION assert_eq4 !BL FUNCTION assert_eqn(nn,string) CHARACTER(LEN=*), INTENT(IN) :: string INTEGER, DIMENSION(:), INTENT(IN) :: nn INTEGER :: assert_eqn if (all(nn(2:) == nn(1))) then assert_eqn=nn(1) else write (*,*) 'nrerror: an assert_eq failed with this tag:', & string STOP 'program terminated by assert_eqn' end if END FUNCTION assert_eqn !BL SUBROUTINE nrerror(string) CHARACTER(LEN=*), INTENT(IN) :: string write (8,*) 'nrerror: ',string STOP 'program terminated by nrerror' END SUBROUTINE nrerror !BL FUNCTION arth_r(first,increment,n) REAL(SP), INTENT(IN) :: first,increment INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(n) :: arth_r INTEGER(I4B) :: k,k2 REAL(SP) :: temp if (n > 0) arth_r(1)=first if (n <= NPAR_ARTH) then do k=2,n arth_r(k)=arth_r(k-1)+increment end do else do k=2,NPAR2_ARTH arth_r(k)=arth_r(k-1)+increment end do temp=increment*NPAR2_ARTH k=NPAR2_ARTH do if (k >= n) exit k2=k+k arth_r(k+1:min(k2,n))=temp+arth_r(1:min(k,n-k)) temp=temp+temp k=k2 end do end if END FUNCTION arth_r !BL FUNCTION arth_d(first,increment,n) REAL(DP), INTENT(IN) :: first,increment INTEGER(I4B), INTENT(IN) :: n REAL(DP), DIMENSION(n) :: arth_d INTEGER(I4B) :: k,k2 REAL(DP) :: temp if (n > 0) arth_d(1)=first if (n <= NPAR_ARTH) then do k=2,n arth_d(k)=arth_d(k-1)+increment end do else do k=2,NPAR2_ARTH arth_d(k)=arth_d(k-1)+increment end do temp=increment*NPAR2_ARTH k=NPAR2_ARTH do if (k >= n) exit k2=k+k arth_d(k+1:min(k2,n))=temp+arth_d(1:min(k,n-k)) temp=temp+temp k=k2 end do end if END FUNCTION arth_d !BL FUNCTION arth_i(first,increment,n) INTEGER(I4B), INTENT(IN) :: first,increment,n INTEGER(I4B), DIMENSION(n) :: arth_i INTEGER(I4B) :: k,k2,temp if (n > 0) arth_i(1)=first if (n <= NPAR_ARTH) then do k=2,n arth_i(k)=arth_i(k-1)+increment end do else do k=2,NPAR2_ARTH arth_i(k)=arth_i(k-1)+increment end do temp=increment*NPAR2_ARTH k=NPAR2_ARTH do if (k >= n) exit k2=k+k arth_i(k+1:min(k2,n))=temp+arth_i(1:min(k,n-k)) temp=temp+temp k=k2 end do end if END FUNCTION arth_i !BL !BL FUNCTION geop_r(first,factor,n) REAL(SP), INTENT(IN) :: first,factor INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(n) :: geop_r INTEGER(I4B) :: k,k2 REAL(SP) :: temp if (n > 0) geop_r(1)=first if (n <= NPAR_GEOP) then do k=2,n geop_r(k)=geop_r(k-1)*factor end do else do k=2,NPAR2_GEOP geop_r(k)=geop_r(k-1)*factor end do temp=factor**NPAR2_GEOP k=NPAR2_GEOP do if (k >= n) exit k2=k+k geop_r(k+1:min(k2,n))=temp*geop_r(1:min(k,n-k)) temp=temp*temp k=k2 end do end if END FUNCTION geop_r !BL FUNCTION geop_d(first,factor,n) REAL(DP), INTENT(IN) :: first,factor INTEGER(I4B), INTENT(IN) :: n REAL(DP), DIMENSION(n) :: geop_d INTEGER(I4B) :: k,k2 REAL(DP) :: temp if (n > 0) geop_d(1)=first if (n <= NPAR_GEOP) then do k=2,n geop_d(k)=geop_d(k-1)*factor end do else do k=2,NPAR2_GEOP geop_d(k)=geop_d(k-1)*factor end do temp=factor**NPAR2_GEOP k=NPAR2_GEOP do if (k >= n) exit k2=k+k geop_d(k+1:min(k2,n))=temp*geop_d(1:min(k,n-k)) temp=temp*temp k=k2 end do end if END FUNCTION geop_d !BL FUNCTION geop_i(first,factor,n) INTEGER(I4B), INTENT(IN) :: first,factor,n INTEGER(I4B), DIMENSION(n) :: geop_i INTEGER(I4B) :: k,k2,temp if (n > 0) geop_i(1)=first if (n <= NPAR_GEOP) then do k=2,n geop_i(k)=geop_i(k-1)*factor end do else do k=2,NPAR2_GEOP geop_i(k)=geop_i(k-1)*factor end do temp=factor**NPAR2_GEOP k=NPAR2_GEOP do if (k >= n) exit k2=k+k geop_i(k+1:min(k2,n))=temp*geop_i(1:min(k,n-k)) temp=temp*temp k=k2 end do end if END FUNCTION geop_i !BL FUNCTION geop_c(first,factor,n) COMPLEX(SP), INTENT(IN) :: first,factor INTEGER(I4B), INTENT(IN) :: n COMPLEX(SP), DIMENSION(n) :: geop_c INTEGER(I4B) :: k,k2 COMPLEX(SP) :: temp if (n > 0) geop_c(1)=first if (n <= NPAR_GEOP) then do k=2,n geop_c(k)=geop_c(k-1)*factor end do else do k=2,NPAR2_GEOP geop_c(k)=geop_c(k-1)*factor end do temp=factor**NPAR2_GEOP k=NPAR2_GEOP do if (k >= n) exit k2=k+k geop_c(k+1:min(k2,n))=temp*geop_c(1:min(k,n-k)) temp=temp*temp k=k2 end do end if END FUNCTION geop_c !BL FUNCTION geop_dv(first,factor,n) REAL(DP), DIMENSION(:), INTENT(IN) :: first,factor INTEGER(I4B), INTENT(IN) :: n REAL(DP), DIMENSION(size(first),n) :: geop_dv INTEGER(I4B) :: k,k2 REAL(DP), DIMENSION(size(first)) :: temp if (n > 0) geop_dv(:,1)=first(:) if (n <= NPAR_GEOP) then do k=2,n geop_dv(:,k)=geop_dv(:,k-1)*factor(:) end do else do k=2,NPAR2_GEOP geop_dv(:,k)=geop_dv(:,k-1)*factor(:) end do temp=factor**NPAR2_GEOP k=NPAR2_GEOP do if (k >= n) exit k2=k+k geop_dv(:,k+1:min(k2,n))=geop_dv(:,1:min(k,n-k))*& spread(temp,2,size(geop_dv(:,1:min(k,n-k)),2)) temp=temp*temp k=k2 end do end if END FUNCTION geop_dv !BL !BL RECURSIVE FUNCTION cumsum_r(arr,seed) RESULT(ans) REAL(SP), DIMENSION(:), INTENT(IN) :: arr REAL(SP), OPTIONAL, INTENT(IN) :: seed REAL(SP), DIMENSION(size(arr)) :: ans INTEGER(I4B) :: n,j REAL(SP) :: sd n=size(arr) if (n == 0_i4b) RETURN sd=0.0_sp if (present(seed)) sd=seed ans(1)=arr(1)+sd if (n < NPAR_CUMSUM) then do j=2,n ans(j)=ans(j-1)+arr(j) end do else ans(2:n:2)=cumsum_r(arr(2:n:2)+arr(1:n-1:2),sd) ans(3:n:2)=ans(2:n-1:2)+arr(3:n:2) end if END FUNCTION cumsum_r !BL RECURSIVE FUNCTION cumsum_i(arr,seed) RESULT(ans) INTEGER(I4B), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), OPTIONAL, INTENT(IN) :: seed INTEGER(I4B), DIMENSION(size(arr)) :: ans INTEGER(I4B) :: n,j,sd n=size(arr) if (n == 0_i4b) RETURN sd=0_i4b if (present(seed)) sd=seed ans(1)=arr(1)+sd if (n < NPAR_CUMSUM) then do j=2,n ans(j)=ans(j-1)+arr(j) end do else ans(2:n:2)=cumsum_i(arr(2:n:2)+arr(1:n-1:2),sd) ans(3:n:2)=ans(2:n-1:2)+arr(3:n:2) end if END FUNCTION cumsum_i !BL !BL RECURSIVE FUNCTION cumprod(arr,seed) RESULT(ans) REAL(SP), DIMENSION(:), INTENT(IN) :: arr REAL(SP), OPTIONAL, INTENT(IN) :: seed REAL(SP), DIMENSION(size(arr)) :: ans INTEGER(I4B) :: n,j REAL(SP) :: sd n=size(arr) if (n == 0_i4b) RETURN sd=1.0_sp if (present(seed)) sd=seed ans(1)=arr(1)*sd if (n < NPAR_CUMPROD) then do j=2,n ans(j)=ans(j-1)*arr(j) end do else ans(2:n:2)=cumprod(arr(2:n:2)*arr(1:n-1:2),sd) ans(3:n:2)=ans(2:n-1:2)*arr(3:n:2) end if END FUNCTION cumprod !BL !BL FUNCTION poly_rr(x,coeffs) REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs REAL(SP) :: poly_rr REAL(SP) :: pow REAL(SP), DIMENSION(:), ALLOCATABLE :: vec INTEGER(I4B) :: i,n,nn n=size(coeffs) if (n <= 0) then poly_rr=0.0_sp else if (n < NPAR_POLY) then poly_rr=coeffs(n) do i=n-1,1,-1 poly_rr=x*poly_rr+coeffs(i) end do else allocate(vec(n+1)) pow=x vec(1:n)=coeffs do vec(n+1)=0.0_sp nn=ishft(n+1,-1) vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) if (nn == 1) exit pow=pow*pow n=nn end do poly_rr=vec(1) deallocate(vec) end if END FUNCTION poly_rr !BL FUNCTION poly_dd(x,coeffs) REAL(DP), INTENT(IN) :: x REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs REAL(DP) :: poly_dd REAL(DP) :: pow REAL(DP), DIMENSION(:), ALLOCATABLE :: vec INTEGER(I4B) :: i,n,nn n=size(coeffs) if (n <= 0) then poly_dd=0.0_dp else if (n < NPAR_POLY) then poly_dd=coeffs(n) do i=n-1,1,-1 poly_dd=x*poly_dd+coeffs(i) end do else allocate(vec(n+1)) pow=x vec(1:n)=coeffs do vec(n+1)=0.0_dp nn=ishft(n+1,-1) vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) if (nn == 1) exit pow=pow*pow n=nn end do poly_dd=vec(1) deallocate(vec) end if END FUNCTION poly_dd !BL FUNCTION poly_rc(x,coeffs) COMPLEX(SPC), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs COMPLEX(SPC) :: poly_rc COMPLEX(SPC) :: pow COMPLEX(SPC), DIMENSION(:), ALLOCATABLE :: vec INTEGER(I4B) :: i,n,nn n=size(coeffs) if (n <= 0) then poly_rc=0.0_sp else if (n < NPAR_POLY) then poly_rc=coeffs(n) do i=n-1,1,-1 poly_rc=x*poly_rc+coeffs(i) end do else allocate(vec(n+1)) pow=x vec(1:n)=coeffs do vec(n+1)=0.0_sp nn=ishft(n+1,-1) vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) if (nn == 1) exit pow=pow*pow n=nn end do poly_rc=vec(1) deallocate(vec) end if END FUNCTION poly_rc !BL FUNCTION poly_cc(x,coeffs) COMPLEX(SPC), INTENT(IN) :: x COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: coeffs COMPLEX(SPC) :: poly_cc COMPLEX(SPC) :: pow COMPLEX(SPC), DIMENSION(:), ALLOCATABLE :: vec INTEGER(I4B) :: i,n,nn n=size(coeffs) if (n <= 0) then poly_cc=0.0_sp else if (n < NPAR_POLY) then poly_cc=coeffs(n) do i=n-1,1,-1 poly_cc=x*poly_cc+coeffs(i) end do else allocate(vec(n+1)) pow=x vec(1:n)=coeffs do vec(n+1)=0.0_sp nn=ishft(n+1,-1) vec(1:nn)=vec(1:n:2)+pow*vec(2:n+1:2) if (nn == 1) exit pow=pow*pow n=nn end do poly_cc=vec(1) deallocate(vec) end if END FUNCTION poly_cc !BL FUNCTION poly_rrv(x,coeffs) REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs,x REAL(SP), DIMENSION(size(x)) :: poly_rrv INTEGER(I4B) :: i,n,m m=size(coeffs) n=size(x) if (m <= 0) then poly_rrv=0.0_sp else if (m < n .or. m < NPAR_POLY) then poly_rrv=coeffs(m) do i=m-1,1,-1 poly_rrv=x*poly_rrv+coeffs(i) end do else do i=1,n poly_rrv(i)=poly_rr(x(i),coeffs) end do end if END FUNCTION poly_rrv !BL FUNCTION poly_ddv(x,coeffs) REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs,x REAL(DP), DIMENSION(size(x)) :: poly_ddv INTEGER(I4B) :: i,n,m m=size(coeffs) n=size(x) if (m <= 0) then poly_ddv=0.0_dp else if (m < n .or. m < NPAR_POLY) then poly_ddv=coeffs(m) do i=m-1,1,-1 poly_ddv=x*poly_ddv+coeffs(i) end do else do i=1,n poly_ddv(i)=poly_dd(x(i),coeffs) end do end if END FUNCTION poly_ddv !BL FUNCTION poly_msk_rrv(x,coeffs,mask) REAL(SP), DIMENSION(:), INTENT(IN) :: coeffs,x LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask REAL(SP), DIMENSION(size(x)) :: poly_msk_rrv poly_msk_rrv=unpack(poly_rrv(pack(x,mask),coeffs),mask,0.0_sp) END FUNCTION poly_msk_rrv !BL FUNCTION poly_msk_ddv(x,coeffs,mask) REAL(DP), DIMENSION(:), INTENT(IN) :: coeffs,x LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: mask REAL(DP), DIMENSION(size(x)) :: poly_msk_ddv poly_msk_ddv=unpack(poly_ddv(pack(x,mask),coeffs),mask,0.0_dp) END FUNCTION poly_msk_ddv !BL !BL RECURSIVE FUNCTION poly_term_rr(a,b) RESULT(u) REAL(SP), DIMENSION(:), INTENT(IN) :: a REAL(SP), INTENT(IN) :: b REAL(SP), DIMENSION(size(a)) :: u INTEGER(I4B) :: n,j n=size(a) if (n <= 0) RETURN u(1)=a(1) if (n < NPAR_POLYTERM) then do j=2,n u(j)=a(j)+b*u(j-1) end do else u(2:n:2)=poly_term_rr(a(2:n:2)+a(1:n-1:2)*b,b*b) u(3:n:2)=a(3:n:2)+b*u(2:n-1:2) end if END FUNCTION poly_term_rr !BL RECURSIVE FUNCTION poly_term_cc(a,b) RESULT(u) COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a COMPLEX(SPC), INTENT(IN) :: b COMPLEX(SPC), DIMENSION(size(a)) :: u INTEGER(I4B) :: n,j n=size(a) if (n <= 0) RETURN u(1)=a(1) if (n < NPAR_POLYTERM) then do j=2,n u(j)=a(j)+b*u(j-1) end do else u(2:n:2)=poly_term_cc(a(2:n:2)+a(1:n-1:2)*b,b*b) u(3:n:2)=a(3:n:2)+b*u(2:n-1:2) end if END FUNCTION poly_term_cc !BL !BL FUNCTION zroots_unity(n,nn) INTEGER(I4B), INTENT(IN) :: n,nn COMPLEX(SPC), DIMENSION(nn) :: zroots_unity INTEGER(I4B) :: k REAL(SP) :: theta zroots_unity(1)=1.0 theta=TWOPI/n k=1 do if (k >= nn) exit zroots_unity(k+1)=cmplx(cos(k*theta),sin(k*theta),SPC) zroots_unity(k+2:min(2*k,nn))=zroots_unity(k+1)*& zroots_unity(2:min(k,nn-k)) k=2*k end do END FUNCTION zroots_unity !BL FUNCTION outerprod_r(a,b) REAL(SP), DIMENSION(:), INTENT(IN) :: a,b REAL(SP), DIMENSION(size(a),size(b)) :: outerprod_r outerprod_r = spread(a,dim=2,ncopies=size(b)) * & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerprod_r !BL FUNCTION outerprod_d(a,b) REAL(DP), DIMENSION(:), INTENT(IN) :: a,b REAL(DP), DIMENSION(size(a),size(b)) :: outerprod_d outerprod_d = spread(a,dim=2,ncopies=size(b)) * & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerprod_d !BL FUNCTION outerdiv(a,b) REAL(SP), DIMENSION(:), INTENT(IN) :: a,b REAL(SP), DIMENSION(size(a),size(b)) :: outerdiv outerdiv = spread(a,dim=2,ncopies=size(b)) / & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerdiv !BL FUNCTION outersum(a,b) REAL(SP), DIMENSION(:), INTENT(IN) :: a,b REAL(SP), DIMENSION(size(a),size(b)) :: outersum outersum = spread(a,dim=2,ncopies=size(b)) + & spread(b,dim=1,ncopies=size(a)) END FUNCTION outersum !BL FUNCTION outerdiff_r(a,b) REAL(SP), DIMENSION(:), INTENT(IN) :: a,b REAL(SP), DIMENSION(size(a),size(b)) :: outerdiff_r outerdiff_r = spread(a,dim=2,ncopies=size(b)) - & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerdiff_r !BL FUNCTION outerdiff_d(a,b) REAL(DP), DIMENSION(:), INTENT(IN) :: a,b REAL(DP), DIMENSION(size(a),size(b)) :: outerdiff_d outerdiff_d = spread(a,dim=2,ncopies=size(b)) - & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerdiff_d !BL FUNCTION outerdiff_i(a,b) INTEGER(I4B), DIMENSION(:), INTENT(IN) :: a,b INTEGER(I4B), DIMENSION(size(a),size(b)) :: outerdiff_i outerdiff_i = spread(a,dim=2,ncopies=size(b)) - & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerdiff_i !BL FUNCTION outerand(a,b) LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: a,b LOGICAL(LGT), DIMENSION(size(a),size(b)) :: outerand outerand = spread(a,dim=2,ncopies=size(b)) .and. & spread(b,dim=1,ncopies=size(a)) END FUNCTION outerand !BL SUBROUTINE scatter_add_r(dest,source,dest_index) REAL(SP), DIMENSION(:), INTENT(OUT) :: dest REAL(SP), DIMENSION(:), INTENT(IN) :: source INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index INTEGER(I4B) :: m,n,j,i n=assert_eq2(size(source),size(dest_index),'scatter_add_r') m=size(dest) do j=1,n i=dest_index(j) if (i > 0 .and. i <= m) dest(i)=dest(i)+source(j) end do END SUBROUTINE scatter_add_r SUBROUTINE scatter_add_d(dest,source,dest_index) REAL(DP), DIMENSION(:), INTENT(OUT) :: dest REAL(DP), DIMENSION(:), INTENT(IN) :: source INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index INTEGER(I4B) :: m,n,j,i n=assert_eq2(size(source),size(dest_index),'scatter_add_d') m=size(dest) do j=1,n i=dest_index(j) if (i > 0 .and. i <= m) dest(i)=dest(i)+source(j) end do END SUBROUTINE scatter_add_d SUBROUTINE scatter_max_r(dest,source,dest_index) REAL(SP), DIMENSION(:), INTENT(OUT) :: dest REAL(SP), DIMENSION(:), INTENT(IN) :: source INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index INTEGER(I4B) :: m,n,j,i n=assert_eq2(size(source),size(dest_index),'scatter_max_r') m=size(dest) do j=1,n i=dest_index(j) if (i > 0 .and. i <= m) dest(i)=max(dest(i),source(j)) end do END SUBROUTINE scatter_max_r SUBROUTINE scatter_max_d(dest,source,dest_index) REAL(DP), DIMENSION(:), INTENT(OUT) :: dest REAL(DP), DIMENSION(:), INTENT(IN) :: source INTEGER(I4B), DIMENSION(:), INTENT(IN) :: dest_index INTEGER(I4B) :: m,n,j,i n=assert_eq2(size(source),size(dest_index),'scatter_max_d') m=size(dest) do j=1,n i=dest_index(j) if (i > 0 .and. i <= m) dest(i)=max(dest(i),source(j)) end do END SUBROUTINE scatter_max_d !BL SUBROUTINE diagadd_rv(mat,diag) REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat REAL(SP), DIMENSION(:), INTENT(IN) :: diag INTEGER(I4B) :: j,n n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagadd_rv') do j=1,n mat(j,j)=mat(j,j)+diag(j) end do END SUBROUTINE diagadd_rv !BL SUBROUTINE diagadd_r(mat,diag) REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat REAL(SP), INTENT(IN) :: diag INTEGER(I4B) :: j,n n = min(size(mat,1),size(mat,2)) do j=1,n mat(j,j)=mat(j,j)+diag end do END SUBROUTINE diagadd_r !BL SUBROUTINE diagmult_rv(mat,diag) REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat REAL(SP), DIMENSION(:), INTENT(IN) :: diag INTEGER(I4B) :: j,n n = assert_eq2(size(diag),min(size(mat,1),size(mat,2)),'diagmult_rv') do j=1,n mat(j,j)=mat(j,j)*diag(j) end do END SUBROUTINE diagmult_rv !BL SUBROUTINE diagmult_r(mat,diag) REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat REAL(SP), INTENT(IN) :: diag INTEGER(I4B) :: j,n n = min(size(mat,1),size(mat,2)) do j=1,n mat(j,j)=mat(j,j)*diag end do END SUBROUTINE diagmult_r !BL FUNCTION get_diag_rv(mat) REAL(SP), DIMENSION(:,:), INTENT(IN) :: mat REAL(SP), DIMENSION(size(mat,1)) :: get_diag_rv INTEGER(I4B) :: j j=assert_eq2(size(mat,1),size(mat,2),'get_diag_rv') do j=1,size(mat,1) get_diag_rv(j)=mat(j,j) end do END FUNCTION get_diag_rv !BL FUNCTION get_diag_dv(mat) REAL(DP), DIMENSION(:,:), INTENT(IN) :: mat REAL(DP), DIMENSION(size(mat,1)) :: get_diag_dv INTEGER(I4B) :: j j=assert_eq2(size(mat,1),size(mat,2),'get_diag_dv') do j=1,size(mat,1) get_diag_dv(j)=mat(j,j) end do END FUNCTION get_diag_dv !BL SUBROUTINE put_diag_rv(diagv,mat) REAL(SP), DIMENSION(:), INTENT(IN) :: diagv REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat INTEGER(I4B) :: j,n n=assert_eq2(size(diagv),min(size(mat,1),size(mat,2)),'put_diag_rv') do j=1,n mat(j,j)=diagv(j) end do END SUBROUTINE put_diag_rv !BL SUBROUTINE put_diag_r(scal,mat) REAL(SP), INTENT(IN) :: scal REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: mat INTEGER(I4B) :: j,n n = min(size(mat,1),size(mat,2)) do j=1,n mat(j,j)=scal end do END SUBROUTINE put_diag_r !BL SUBROUTINE unit_matrix(mat) REAL(SP), DIMENSION(:,:), INTENT(OUT) :: mat INTEGER(I4B) :: i,n n=min(size(mat,1),size(mat,2)) mat(:,:)=0.0_sp do i=1,n mat(i,i)=1.0_sp end do END SUBROUTINE unit_matrix !BL FUNCTION upper_triangle(j,k,extra) INTEGER(I4B), INTENT(IN) :: j,k INTEGER(I4B), OPTIONAL, INTENT(IN) :: extra LOGICAL(LGT), DIMENSION(j,k) :: upper_triangle INTEGER(I4B) :: n n=0 if (present(extra)) n=extra upper_triangle=(outerdiff(arth_i(1,1,j),arth_i(1,1,k)) < n) END FUNCTION upper_triangle !BL FUNCTION lower_triangle(j,k,extra) INTEGER(I4B), INTENT(IN) :: j,k INTEGER(I4B), OPTIONAL, INTENT(IN) :: extra LOGICAL(LGT), DIMENSION(j,k) :: lower_triangle INTEGER(I4B) :: n n=0 if (present(extra)) n=extra lower_triangle=(outerdiff(arth_i(1,1,j),arth_i(1,1,k)) > -n) END FUNCTION lower_triangle !BL FUNCTION vabs(v) REAL(SP), DIMENSION(:), INTENT(IN) :: v REAL(SP) :: vabs vabs=sqrt(dot_product(v,v)) END FUNCTION vabs !BL END MODULE nrutil module constants_module !definition de quelques constantes utiles use nrtype implicit none integer, public, parameter :: large_int_kind = selected_real_kind(4), & ! large_real_kind = selected_real_kind(8) large_real_kind = kind(1.0D0) real(SP), public, parameter :: & trois=3.0_sp, demi=.5_sp, & two_pi = 2*pi, un=1.0_sp, deux=2.0_sp, zero=0.0_sp end module constants_module module mod_mat ! spectral_radius(a) donne le rayon spectral de a ! spectre(a) donne les valeurs propres de a ! adamar(x,y) est le produit d'adamar de 2 vecteurs ! prod(x,y) est la matrice produit x*y' de 2 vecteurs ! norm(a) est la norme somme des valeurs absolue de la matrice a ! normv(a) est la norme somme des valeurs absolue du vecteur a ! identity(ia) est la matrice identite de taille ia ! diag(vec) est la matrice diagonale dont la premiere diagonale est vec ! trace(a) est la trace de la matrice carree a ! kro(a,b) est le produit de kro de a et b ! gaussj(a,b) resoud ax=b ou a de taille n*n et b de taille n*m ! en sortie a contient l'inverse de a et b la solution ! ERREUR: gauss(a,b,n,m) resoud ax=b ou a de taille n*n et b de taille n*m ! en sortie a contient l'inverse de a et b la solution ! ERREUR: inv(a) contient l'inverse d'une matrice a de taille n ! inverse(a) contient l'inverse d'une matrice a de taille n ! hh matrice h de taille n(n+1)/2 x n telle vech(a)=hvec(a) ! kk matrice k de taille n(n+1)/2 x n**2 telle vec(a)=k'vech(a) ! vech d'une matrice a symetrique de taille n x n ! vec d'une matrice a de taille n x m ! jordan(A,P,Lambda) tel que A=P*Lambda*P' ou PP'=I, Lambda diagonale et A symetrique ! sqrtmat(A) est A**(1/2) ou A symetrique ! sqrtmat2(A) est A**(1/2)=L ou A=LL' symetrique definie positive ! sqrtinv2(A) est L^-1 ou A=LL' symetrique definie positive ! sqrtinv(A) est A**(-1/2) ou A symetrique ! invsim(A) est A**(-1) ou A symetrique ! writemat(A) ecrit la matrice A use constants_module; use nrtype USE nrtype; USE nrutil, ONLY : assert_eq implicit none private private :: elmhes, hqr, vech, hh, kk, normv,& balanc, adamar, tred2,tqli, choldc public :: identity, inverse, vec, spectral_radius,kro, kromatvec,krovecmat,krovecvec,spectre,prod,resgauss,& diag,trace,jordan,sqrtmat,sqrtinv,invsim,sqrtmat2,sqrtinv2,writemat,determinant,norm,gaussj,det contains SUBROUTINE gaussj(a,b) USE nrtype; USE nrutil, ONLY : assert_eq,nrerror,outerand,outerprod,swap IMPLICIT NONE REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b INTEGER(I4B), DIMENSION(size(a,1)) :: ipiv,indxr,indxc LOGICAL(LGT), DIMENSION(size(a,1)) :: lpiv REAL(SP) :: pivinv REAL(SP), DIMENSION(size(a,1)) :: dumc INTEGER(I4B), TARGET :: irc(2) INTEGER(I4B) :: i,l,n INTEGER(I4B), POINTER :: irow,icol n=assert_eq(size(a,1),size(a,2),size(b,1),'gaussj') irow => irc(1) icol => irc(2) ipiv=0 do i=1,n lpiv = (ipiv == 0) irc=maxloc(abs(a),outerand(lpiv,lpiv)) ipiv(icol)=ipiv(icol)+1 if (ipiv(icol) > 1) call nrerror('gaussj: singular matrix (1)') if (irow /= icol) then call swap(a(irow,:),a(icol,:)) call swap(b(irow,:),b(icol,:)) end if indxr(i)=irow indxc(i)=icol if (a(icol,icol) == 0.0) & call nrerror('gaussj: singular matrix (2)') pivinv=1.0_sp/a(icol,icol) a(icol,icol)=1.0 a(icol,:)=a(icol,:)*pivinv b(icol,:)=b(icol,:)*pivinv dumc=a(:,icol) a(:,icol)=0.0 a(icol,icol)=pivinv a(1:icol-1,:)=a(1:icol-1,:)-outerprod(dumc(1:icol-1),a(icol,:)) b(1:icol-1,:)=b(1:icol-1,:)-outerprod(dumc(1:icol-1),b(icol,:)) a(icol+1:,:)=a(icol+1:,:)-outerprod(dumc(icol+1:),a(icol,:)) b(icol+1:,:)=b(icol+1:,:)-outerprod(dumc(icol+1:),b(icol,:)) end do do l=n,1,-1 call swap(a(:,indxr(l)),a(:,indxc(l))) end do END SUBROUTINE gaussj ! writemat(A) ecrirt la matrice A subroutine writemat(A) real (SP), intent(in), dimension(:,:) :: A integer(I4B) :: i, n n=size(A,1) do i=1,n write(8,fmt='(1x,40(f7.3,1x))') A(i,:) end do return end subroutine writemat function inverse(a) real(SP), intent(in) :: a(:,:) real(SP) :: inverse(size(a,1),size(a,1)) real(sp) :: aa(size(a,1),size(a,1)),y(size(a,1),size(a,1)),d integer(I4B) :: j,indx(size(a,1)) aa=a y=identity(size(a,1)) call ludcmp(aa,indx,d) do j=1,size(a,1) call lubksb(aa,indx,y(:,j)) end do inverse=y return contains SUBROUTINE lubksb(a,indx,b) USE nrtype; USE nrutil, ONLY : assert_eq IMPLICIT NONE REAL(SP), DIMENSION(:,:), INTENT(IN) :: a INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx REAL(SP), DIMENSION(:), INTENT(INOUT) :: b INTEGER(I4B) :: i,n,ii,ll REAL(SP) :: summ n=assert_eq(size(a,1),size(a,2),size(indx),'lubksb') ii=0 do i=1,n ll=indx(i) summ=b(ll) b(ll)=b(i) if (ii /= 0) then summ=summ-dot_product(a(i,ii:i-1),b(ii:i-1)) else if (summ /= 0.0) then ii=i end if b(i)=summ end do do i=n,1,-1 b(i) = (b(i)-dot_product(a(i,i+1:n),b(i+1:n)))/a(i,i) end do END SUBROUTINE lubksb SUBROUTINE ludcmp(a,indx,d) USE nrtype; USE nrutil, ONLY : assert_eq,imaxloc,nrerror,outerprod,swap IMPLICIT NONE REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx REAL(SP), INTENT(OUT) :: d REAL(SP), DIMENSION(size(a,1)) :: vv REAL(SP), PARAMETER :: TINY=1.0e-20_sp INTEGER(I4B) :: j,n,imax n=assert_eq(size(a,1),size(a,2),size(indx),'ludcmp') d=1.0 vv=maxval(abs(a),dim=2) if (any(vv == 0.0)) call nrerror('singular matrix in ludcmp') vv=1.0_sp/vv do j=1,n imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j))) if (j /= imax) then call swap(a(imax,:),a(j,:)) d=-d vv(imax)=vv(j) end if indx(j)=imax if (a(j,j) == 0.0) a(j,j)=TINY a(j+1:n,j)=a(j+1:n,j)/a(j,j) a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n)) end do END SUBROUTINE ludcmp end function inverse function det(a) real(SP), intent(in) :: A(:,:) real(SP) :: det real(SP) :: aout(size(a,1),size(a,1)),d integer(I4B) :: i integer(I4B) :: indx(size(a,1)) aout=a call ludcmp(aout,indx,d) det=d do i=1,size(a,1) det=det*aout(i,i) end do return contains SUBROUTINE ludcmp(a,indx,d) USE nrtype; USE nrutil, ONLY : assert_eq,imaxloc,nrerror,outerprod,swap IMPLICIT NONE REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx REAL(SP), INTENT(OUT) :: d REAL(SP), DIMENSION(size(a,1)) :: vv REAL(SP), PARAMETER :: TINY=1.0e-20_sp INTEGER(I4B) :: j,n,imax n=assert_eq(size(a,1),size(a,2),size(indx),'ludcmp') d=1.0 vv=maxval(abs(a),dim=2) if (any(vv == 0.0)) call nrerror('singular matrix in ludcmp') vv=1.0_sp/vv do j=1,n imax=(j-1)+imaxloc(vv(j:n)*abs(a(j:n,j))) if (j /= imax) then call swap(a(imax,:),a(j,:)) d=-d vv(imax)=vv(j) end if indx(j)=imax if (a(j,j) == 0.0) a(j,j)=TINY a(j+1:n,j)=a(j+1:n,j)/a(j,j) a(j+1:n,j+1:n)=a(j+1:n,j+1:n)-outerprod(a(j+1:n,j),a(j,j+1:n)) end do END SUBROUTINE ludcmp end function det function determinant(a) real(SP), intent(in) :: A(:,:) real(SP) :: determinant real(SP) :: aout(size(a,1),size(a,1)),d integer(I4B) :: i call ludcmp(a,aout,d) determinant=d do i=1,size(a,1) determinant=determinant*aout(i,i) end do return contains subroutine ludcmp(ain,a,d) real (SP), intent(in) :: ain(:,:) real (SP), intent(out) :: a(:,:),d real (SP), parameter :: tiny=1.0e-20 integer(I4B), parameter :: nmax=100 real (SP) ::indx(size(a,1)),vv(nmax) real (SP) ::aamax,sum,dum integer(I4B) :: n,i,j,k,imax a=ain n=size(a,1) d=un do i=1,n aamax=zero do j=1,n if (abs(a(i,j)).gt.aamax) aamax=abs(a(i,j)) end do if (aamax.eq.0.) write(8,*) "singular matrix." vv(i)=un/aamax end do do j=1,n if (j.gt.1) then do i=1,j-1 sum=a(i,j) if (i.gt.1)then do k=1,i-1 sum=sum-a(i,k)*a(k,j) end do a(i,j)=sum end if end do end if aamax=zero do i=j,n sum=a(i,j) if (j.gt.1)then do k=1,j-1 sum=sum-a(i,k)*a(k,j) end do a(i,j)=sum end if dum=vv(i)*abs(sum) if (dum.ge.aamax) then imax=i aamax=dum end if end do if (j.ne.imax)then do k=1,n dum=a(imax,k) a(imax,k)=a(j,k) a(j,k)=dum end do d=-d vv(imax)=vv(j) end if indx(j)=imax if(j.ne.n)then if(a(j,j).eq.zero)a(j,j)=tiny dum=un/a(j,j) do i=j+1,n a(i,j)=a(i,j)*dum end do end if end do if(a(n,n).eq.zero)a(n,n)=tiny return end subroutine ludcmp end function determinant function sqrtinv2(A) real(SP), intent(in) :: A(:,:) real(SP) :: sqrtinv2(size(a,1),size(a,2)) real(SP) :: p(size(a,1)), b(size(a,1),size(a,2)),sum integer(I4B) :: n,i,j,k n=size(a,1) b=A call choldc(b,n,p) do i=1,n b(i,i)=un/p(i) do j=i+1,n sum=zero do k=i,j-1 sum=sum-b(j,k)*b(k,i) end do b(j,i)=sum/p(j) end do end do do i=1,n sqrtinv2(i,1:i)=b(i,1:i) sqrtinv2(i,i+1:n)=zero end do return end function sqrtinv2 function sqrtmat2(A) real(SP), intent(in) :: A(:,:) real(SP) :: sqrtmat2(size(a,1),size(a,2)) real(SP) :: p(size(a,1)), b(size(a,1),size(a,2)) integer(I4B) :: n,i n=size(a,1) b=A call choldc(b,n,p) do i=1,n sqrtmat2(i,1:i-1)=b(i,1:i-1) sqrtmat2(i,i)=p(i) sqrtmat2(i,i+1:n)=zero end do return end function sqrtmat2 subroutine choldc(a,n,p) real(SP), intent(in out) :: a(:,:) real(SP), intent(out) :: p(:) integer(I4B), intent(in) :: n integer(I4B) :: i,j,k real(SP) :: sum !given a positive-de .nite symmetric matrix a(1:n,1:n),with physical dimension np this !routine constructs its cholesky decomposition,a = l · l t .on input,only the upper triangle !of a need be given;it is not modi ed.the cholesky factor l is returned in the lower triangle !of a except for its diagonal elements which are returned in p(1:n). do i=1,n do j=i,n sum=a(i,j) do k=i-1,1,-1 sum=sum-a(i,k)*a(j,k) end do if(i.eq.j)then if(sum.le.zero)write(8,*)"choldc failed" p(i)=sqrt(sum) else a(j,i)=sum/p(i) end if end do end do return end subroutine choldc function invsim(A) real(SP), intent(in) :: A(:,:) real(SP) :: invsim(size(a,1),size(a,2)) real(SP) :: P(size(a,1),size(a,2)), Lambda(size(a,1),size(a,2)) integer(I4B) :: i,n call jordan(A,P,Lambda) n=size(a,1) do i=1,n Lambda(i,i)=un/Lambda(i,i) end do invsim=matmul(P,matmul(Lambda,transpose(P))) return end function invsim function sqrtinv(A) real(SP), intent(in) :: A(:,:) real(SP) :: sqrtinv(size(a,1),size(a,2)) real(SP) :: P(size(a,1),size(a,2)), Lambda(size(a,1),size(a,2)) integer(I4B) :: i,n call jordan(A,P,Lambda) n=size(a,1) do i=1,n Lambda(i,i)=un/Lambda(i,i) end do sqrtinv=matmul(P,matmul(sqrt(Lambda),transpose(P))) return end function sqrtinv function sqrtmat(A) real(SP), intent(in) :: A(:,:) real(SP) :: sqrtmat(size(a,1),size(a,2)) real(SP) :: P(size(a,1),size(a,2)), Lambda(size(a,1),size(a,2)) call jordan(A,P,Lambda) sqrtmat=matmul(P,matmul(sqrt(Lambda),transpose(P))) return end function sqrtmat subroutine jordan(A,P,Lambda) real(SP), intent(in) :: A(:,:) real(SP), intent(out) :: P(:,:), Lambda(:,:) real(SP) :: d(size(a,1)),e(size(a,1)) P=a call tred2(P,d,e) call tqli(d,e,P) Lambda=diag(d) return end subroutine jordan subroutine tred2(a,d,e) real(SP), intent(in out) :: a(:,:) real(SP), intent(out) :: d(:), e(:) integer(I4B) :: n,i,l,k,j real(SP) :: h,scale,f,g,hh n=size(a,1) if(n.gt.1)then do i=n,2,-1 l=i-1 h=zero scale=zero if(l.gt.1)then do k=1,l scale=scale+abs(a(i,k)) end do if(scale.eq.zero)then e(i)=a(i,l) else do k=1,l a(i,k)=a(i,k)/scale h=h+a(i,k)**2 end do f=a(i,l) g=-sign(sqrt(h),f) e(i)=scale*g h=h-f*g a(i,l)=f-g f=zero do j=1,l a(j,i)=a(i,j)/h g=zero do k=1,j g=g+a(j,k)*a(i,k) end do if(l.gt.j)then do k=j+1,l g=g+a(k,j)*a(i,k) end do end if e(j)=g/h f=f+e(j)*a(i,j) end do hh=f/(h+h) do j=1,l f=a(i,j) g=e(j)-hh*f e(j)=g do k=1,j a(j,k)=a(j,k)-f*e(k)-g*a(i,k) end do end do end if else e(i)=a(i,l) end if d(i)=h end do end if d(1)=zero e(1)=zero do i=1,n l=i-1 if(d(i).ne.zero)then do j=1,l g=zero do k=1,l g=g+a(i,k)*a(k,j) end do do k=1,l a(k,j)=a(k,j)-g*a(k,i) end do end do end if d(i)=a(i,i) a(i,i)=un if(l.ge.1)then do j=1,l a(i,j)=zero a(j,i)=zero end do end if end do return end subroutine tred2 subroutine tqli(d,e,z) real (SP), intent(in out) :: d(:),e(:) real (SP), intent(out) :: z(:,:) integer(I4B) :: n,i,l,iter,m,k real(SP) :: dd,g,r,s,c,p,f,b n=size(d) ! z=identity(n) if (n.gt.1) then do i=2,n e(i-1)=e(i) end do e(n)=zero do l=1,n iter=0 1 do m=l,n-1 dd=abs(d(m))+abs(d(m+1)) if (abs(e(m))+dd.eq.dd) go to 2 end do m=n 2 if(m.ne.l)then if(iter.eq.30)write(8,*) 'too many iterations' iter=iter+1 g=(d(l+1)-d(l))/(deux*e(l)) r=sqrt(g**2+un) g=d(m)-d(l)+e(l)/(g+sign(r,g)) s=un c=un p=zero do i=m-1,l,-1 f=s*e(i) b=c*e(i) if(abs(f).ge.abs(g))then c=g/f r=sqrt(c**2+un) e(i+1)=f*r s=un/r c=c*s else s=f/g r=sqrt(s**2+un) e(i+1)=g*r c=un/r s=s*c end if g=d(i+1)-p r=(d(i)-g)*s+deux*c*b p=s*r d(i+1)=g+p g=c*r-b do k=1,n f=z(k,i+1) z(k,i+1)=s*z(k,i)+c*f z(k,i)=c*z(k,i)-s*f end do end do d(l)=d(l)-p e(l)=g e(m)=zero go to 1 end if end do end if return end subroutine tqli ! prod(x,y) est la matrice produit x*y' de 2 vecteurs function prod(x,y) result (prod_result) real (SP), intent(in), dimension(:) :: x, y real (SP), dimension(size(x),size(y)) :: prod_result integer(I4B) :: i, j, n, m n=size(x) m=size(y) do i=1,n do j=1,m prod_result(i,j)=x(i)*y(j) end do end do return end function prod function vech(a) result (vech_result) ! vech d'une matrice a symetrique de taille n x n real(SP), intent(in), dimension(:,:) :: a real(SP), dimension((size(a,1)*(size(a,1)+1))/2) :: vech_result integer(I4B) :: i, j, id1, n n=size(a,1) id1=0 do j=1,n do i=j,n id1=id1+1 vech_result(id1)=a(i,j) end do end do return end function vech function vec(a) result (vec_result) ! vec d'une matrice a de taille n x m real(SP), intent(in), dimension(:,:) :: a real(SP), dimension(size(a,1)*size(a,2)) :: vec_result integer(I4B) :: i, j, n, m n=size(a,1) m=size(a,2) do j=1, m do i=1,n vec_result((j-1)*n+i)=a(i,j) end do end do return end function vec function kk(n) result (kk_result) ! matrice k de taille n(n+1)/2 x n**2 telle vec(a)=k'vech(a) integer(I4B), intent(in) :: n real(SP), dimension((n*(n+1))/2,n**2) :: kk_result real(SP), dimension(n**2,(n*(n+1))/2) :: k integer(I4B) :: i, j k=0.0_sp do j=1,n do i=j,n k((j-1)*n+i,((j-1)*(2*n-j))/2+i)=1.0_sp k((i-1)*n+j,((j-1)*(2*n-j))/2+i)=1.0_sp end do end do kk_result=transpose(k) return end function kk function hh(n) result (hh_result) ! matrice h de taille n(n+1)/2 x n**2 telle vech(a)=hvec(a) integer(I4B), intent(in) :: n real(SP), dimension((n*(n+1))/2,n**2) :: hh_result integer(I4B) :: i, j, id1, id2 hh_result=0.0_sp id1=0 id2=0 do i=n,1,-1 id2=id2+n id1=id1+i do j=1,i hh_result(id1+1-j,id2+1-j)=1.0_sp end do end do return end function hh function invgauss(aa,bb,n,m) result (invgauss_result) ! resoud aax=bb ou aa (n,n) x (n,m) et bb (n,m) ! en sortie invgauss_result est remplacee par l'inverse de a integer(I4B), intent(in) :: n, m real (SP), intent(in), dimension(:,:) :: aa, bb real (SP), dimension(size(aa,1),size(aa,2)) :: invgauss_result ! local variables real (SP), dimension(size(aa,1),size(aa,2)) :: a real (SP), dimension(size(bb,1),size(bb,2)) :: b real (SP) :: big, dum, pivinv integer(I4B), dimension(n) :: ipiv, indxr, indxc integer(I4B) :: i, j, k, l, ll, irow, icol a=aa b=bb ipiv=0 do i=1,n big=0.0_sp do j=1,n if (ipiv(j) /= 1) then do k=1,n if (ipiv(k) == 0) then if (abs(a(j,k)) >= big) then big=abs(a(j,k)) irow=j icol=k end if else if (ipiv(k) > 1) then write(unit=4,fmt=*)"matrice singuliere" return end if end if end do end if end do ipiv(icol)=ipiv(icol)+1 if (irow /= icol) then do l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum end do do l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum end do end if indxr(i)=irow indxc(i)=icol if (a(icol,icol)== 0.0_sp) then write(unit=4,fmt=*)"matrice singuliere" return end if pivinv=1.0_sp/a(icol,icol) a(icol,icol)=1.0_sp do l=1,n a(icol,l)=a(icol,l)*pivinv end do do l=1,m b(icol,l)=b(icol,l)*pivinv end do do ll=1,n if (ll /= icol) then dum=a(ll,icol) a(ll,icol)=0.0_sp do l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum end do do l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum end do end if end do end do do l=n,1,-1 if (indxr(l) /= indxr(l)) then do k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum end do end if end do invgauss_result=a return end function invgauss function inv(a) result (inv_result) real (SP), intent(in), dimension(:,:) :: a real (SP), dimension(size(a,1),size(a,1)) :: inv_result real (SP), dimension(size(a,1),1) :: zero zero=0.0_sp inv_result=invgauss(a,zero,size(a,1),1) return end function inv function balanc(aa,n) result (balanc_result) real(SP), intent(in), dimension(:,:) :: aa integer(I4B), intent(in) :: n real(SP), dimension(size(aa,1),size(aa,2)) :: balanc_result real(SP), parameter :: rradix=2.0_sp real(SP), parameter :: sqrdx=rradix**2 real(SP), dimension(size(aa,1),size(aa,2)) :: a integer(I4B) :: i, j, last real(SP) :: c, r, g, f, s a=aa boucle1 : do last=1 do i=1,n c=0.0_sp r=0.0_sp do j=1,n if (j /= i) then c=c+abs(a(j,i)) r=r+abs(a(i,j)) end if end do if (c /= 0.0_sp .and. r /= 0.0_sp)then g=r/rradix f=1.0_sp s=c+r boucle2 : do if (c < g) then f=f*rradix c=c*sqrdx else exit boucle2 end if end do boucle2 g=r*rradix boucle3 : do if (c > g) then f=f/rradix c=c/rradix else exit boucle3 end if end do boucle3 if ((c+r)/f < 0.950_sp*s) then last=0 g=1.0_sp/f do j=1,n a(i,j)=a(i,j)*g end do do j=1,n a(j,i)=a(j,i)*f end do end if end if end do if (last /= 0)then exit boucle1 end if end do boucle1 balanc_result=a return end function balanc function elmhes(aa,n) result(elmhes_result) real (SP), intent(in), dimension(:,:) :: aa integer(I4B), intent(in) :: n real (SP), dimension(size(aa,1),size(aa,2)) :: elmhes_result real (SP), dimension(size(aa,1),size(aa,2)) :: a integer(I4B) :: m, i, j real (SP) :: x, y a=aa a=balanc(a,n) do m=2,n-1 x=0.0_sp i=m do j=m,n if (abs(a(j,m-1)) > abs(x)) then x=a(j,m-1) i=j end if end do if (i /= m) then do j=m-1,n y=a(i,j) a(i,j)=a(m,j) a(m,j)=y end do do j=1,n y=a(j,i) a(j,i)=a(j,m) a(j,m)=y end do end if if (x /= 0.00_sp) then do i=m+1,n y=a(i,m-1) if (y /= 0.00_sp) then y=y/x a(i,m-1)=y do j=m,n a(i,j)=a(i,j)-y*a(m,j) end do do j=1,n a(j,m)=a(j,m)+y*a(j,i) end do end if end do end if end do elmhes_result=a return end function elmhes function hqr(aa,n) result(hqr_result) real(SP), intent(in), dimension(:,:) :: aa integer(I4B), intent(in) :: n complex(SP), dimension(size(aa,1)) :: hqr_result real(SP), dimension(size(aa,1),size(aa,2)) :: a real(SP), dimension(size(aa,1)) :: wr, wi integer(I4B) :: m, i, j, k, nn, l, its,ii real(SP) :: anorm, t, s, x, y, w, p, q, z, r, u, v a=aa anorm=abs(a(1,1)) do i=2,n do j=i-1,n anorm=anorm+abs(a(i,j)) end do end do nn=n t=0.0_sp boucle1 : do if (nn >= 1) then its=0 boucle2 : do boucle1x :do ii=1,1 do l=nn,2,-1 s=abs(a(l-1,l-1))+abs(a(l,l)) if (s == 0.0_sp)then s=anorm end if if (abs(a(l,l-1))+s == s) then exit boucle1x end if end do l=1 end do boucle1x x=a(nn,nn) if (l == nn) then wr(nn)=x+t wi(nn)=0.0_sp nn=nn-1 exit boucle2 else y=a(nn-1,nn-1) w=a(nn,nn-1)*a(nn-1,nn) if (l == nn-1) then p=0.50_sp*(y-x) q=p**2+w z=sqrt(abs(q)) x=x+t if (q >= 0.0_sp) then z=p+sign(z,p) wr(nn)=x+z wr(nn-1)=wr(nn) if (z /= 0.0_sp) then wr(nn)=x-w/z end if wi(nn)=0.0_sp wi(nn-1)=0.0_sp else wr(nn)=x+p wr(nn-1)=wr(nn) wi(nn)=z wi(nn-1)=-z end if nn=nn-2 exit boucle2 else if (its == 30) then hqr_result(:)=cmplx(wr(:),wi(:),sp) write(unit=4,fmt=*)"hqr fails: too many iterations" return end if if (its == 10 .or. its == 20) then t=t+x do i=1,nn a(i,i)=a(i,i)-x end do s=abs(a(nn,nn-1))+abs(a(nn-1,nn-2)) x=0.750_sp*s y=x w=-0.43750_sp*s**2 end if its=its+1 do m=nn-2,l,-1 z=a(m,m) r=x-z s=y-z p=(r*s-w)/a(m+1,m)+a(m,m+1) q=a(m+1,m+1)-z-r-s r=a(m+2,m+1) s=abs(p)+abs(q)+abs(r) p=p/s q=q/s r=r/s if (m == l) then exit end if u=abs(a(m,m-1))*(abs(q)+abs(r)) v=abs(p)*(abs(a(m-1,m-1))+abs(z)+abs(a(m+1,m+1))) if( u+v == v) then exit end if end do do i=m+2,nn a(i,i-2)=0.0_sp if (i /= m+2)then a(i,i-3)=0.0_sp end if end do do k=m,nn-1 if (k /= m) then p=a(k,k-1) q=a(k+1,k-1) r=0.0_sp if (k /= nn-1) then r=a(k+2,k-1) end if x=abs(p)+abs(q)+abs(r) if (x /= 0.0_sp) then p=p/x q=q/x r=r/x end if end if s=sign(sqrt(p**2+q**2+r**2),p) if (s /= 0.0_sp) then if(k==m)then if (l /= m) then a(k,k-1)=-a(k,k-1) end if else a(k,k-1)=-s*x end if p=p+s x=p/s y=q/s z=r/s q=q/p r=r/p do j=k,nn p=a(k,j)+q*a(k+1,j) if (k /= nn-1) then p=p+r*a(k+2,j) a(k+2,j)=a(k+2,j)-p*z end if a(k+1,j)=a(k+1,j)-p*y a(k,j)=a(k,j)-p*x end do do i=l,min(nn,k+3) p=x*a(i,k)+y*a(i,k+1) if (k /= nn-1) then p=p+z*a(i,k+2) a(i,k+2)=a(i,k+2)-p*r end if a(i,k+1)=a(i,k+1)-p*q a(i,k)=a(i,k)-p end do end if end do end if end if end do boucle2 else exit boucle1 end if end do boucle1 hqr_result(:)=cmplx(wr(:),wi(:),sp) return end function hqr function spectral_radius(a) result (spectral_result) ! spectral_radius(a,n) donne le rayon spectral de a carree de taille n real (SP), intent(in), dimension(:,:) :: a real (SP) :: spectral_result real (SP), dimension(size(a,1),size(a,2)) :: b integer(I4B) :: n n=size(a,1) b= elmhes(a,n) spectral_result=maxval(abs(hqr(b,n))) return end function spectral_radius function spectre(a) result(spectre_result) ! spectre(a) donne le rayon spectral de a carree de taille n real (SP), intent(in), dimension(:,:) :: a complex(SP), dimension(size(a,1)) :: spectre_result real (SP), dimension(size(a,1),size(a,1)) :: b integer(I4B) :: n n=size(a,1) b=elmhes(a,n) spectre_result=hqr(b,n) return end function spectre function adamar(x,y) result(adamar_result) real (SP), intent(in), dimension(:) :: x, y real (SP), dimension(size(x)) :: adamar_result integer(I4B) :: i, n n=size(x) do i=1,n adamar_result(i)=x(i)*y(i) end do return end function adamar function resgauss(aa,bb,n,m) result(resgauss_result) ! resoud aax=bb ou a (n,n) x (n,m) et bb (n,m) ! en sortie resgauss_result contient la solution integer(I4B), intent(in) :: n, m real (SP), intent(in), dimension(:,:) :: aa, bb real (SP),dimension(n,m) :: resgauss_result ! local variables real (SP) :: big, dum, pivinv integer(I4B), dimension(n) :: ipiv, indxr, indxc real (SP), dimension(n,n) :: a real (SP), dimension(n,m) :: b integer(I4B) :: i, j, k, l, ll, irow, icol a=aa b=bb ipiv=0 do i=1,n big=0.0_sp do j=1,n if (ipiv(j) /= 1) then do k=1,n if (ipiv(k) == 0) then if (abs(a(j,k)) >= big) then big=abs(a(j,k)) irow=j icol=k end if else if (ipiv(k) > 1) then resgauss_result=b write(unit=4,fmt=*)"matrice singuliere" return end if end if end do end if end do ipiv(icol)=ipiv(icol)+1 if (irow /= icol) then do l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum end do do l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum end do end if indxr(i)=irow indxc(i)=icol if (a(icol,icol)== 0.0_sp) then resgauss_result=b write(unit=4,fmt=*)"matrice singuliere" return end if pivinv=1.0_sp/a(icol,icol) a(icol,icol)=1.0_sp do l=1,n a(icol,l)=a(icol,l)*pivinv end do do l=1,m b(icol,l)=b(icol,l)*pivinv end do do ll=1,n if (ll /= icol) then dum=a(ll,icol) a(ll,icol)=0.0_sp do l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum end do do l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum end do end if end do end do do l=n,1,-1 if (indxr(l) /= indxr(l)) then do k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum end do end if end do resgauss_result=b return end function resgauss subroutine gauss(a,b,n,m) ! resoud ax=b ou a (n,n) x (n,m) et b (n,m) ! en sortie a est remplacee par l'inverse de a et b contient la solution integer(I4B), intent(in) :: n, m real (SP), intent(in out), dimension(:,:) :: a, b ! local variables real (SP) :: big, dum, pivinv integer(I4B), dimension(n) :: ipiv, indxr, indxc integer(I4B) :: i, j, k, l, ll, irow, icol ipiv=0 do i=1,n big=0.0_sp do j=1,n if (ipiv(j) /= 1) then do k=1,n if (ipiv(k) == 0) then if (abs(a(j,k)) >= big) then big=abs(a(j,k)) irow=j icol=k end if else if (ipiv(k) > 1) then write(unit=4,fmt=*)"matrice singuliere" return end if end if end do end if end do ipiv(icol)=ipiv(icol)+1 if (irow /= icol) then do l=1,n dum=a(irow,l) a(irow,l)=a(icol,l) a(icol,l)=dum end do do l=1,m dum=b(irow,l) b(irow,l)=b(icol,l) b(icol,l)=dum end do end if indxr(i)=irow indxc(i)=icol if (a(icol,icol)== 0.0_sp) then write(unit=4,fmt=*)"matrice singuliere" return end if pivinv=1.0_sp/a(icol,icol) a(icol,icol)=1.0_sp do l=1,n a(icol,l)=a(icol,l)*pivinv end do do l=1,m b(icol,l)=b(icol,l)*pivinv end do do ll=1,n if (ll /= icol) then dum=a(ll,icol) a(ll,icol)=0.0_sp do l=1,n a(ll,l)=a(ll,l)-a(icol,l)*dum end do do l=1,m b(ll,l)=b(ll,l)-b(icol,l)*dum end do end if end do end do do l=n,1,-1 if (indxr(l) /= indxr(l)) then do k=1,n dum=a(k,indxr(l)) a(k,indxr(l))=a(k,indxc(l)) a(k,indxc(l))=dum end do end if end do return end subroutine gauss function krovecvec(a,b) result(kro_result) !!!! ! kro=produit de kronecker du vecteur a et du vecteur b !!!! real (SP), intent(in), dimension(:) :: a real (SP), intent(in), dimension(:) :: b real (SP), dimension(size(a,1)*size(b)) :: kro_result integer(I4B) :: i, ia, ib ia=size(a) ib=size(b) do i=1,ia kro_result((i-1)*ia+1:(i-1)*ia+ib)=a(i)*b(1:ib) end do return end function krovecvec function krovecmat(b,a) result(kro_result) !!!! ! kro=produit de kronecker du vecteur b et de la matrice a !!!! real (SP), intent(in), dimension(:,:) :: a real (SP), intent(in), dimension(:) :: b real (SP), dimension(size(a,1)*size(b),size(a,2)) :: kro_result integer(I4B) :: i1, ia, ib, ja ia=size(a,1) ib=size(b) ja=size(a,2) do i1=1,ib kro_result((i1-1)*ib+1:(i1-1)*ib+ia,1:ja)=b(i1)*a(1:ia,1:ja) end do return end function krovecmat function kromatvec(a,b) result(kro_result) !!!! ! kro=produit de kronecker de la matrice a et du vecteur b !!!! real (SP), intent(in), dimension(:,:) :: a real (SP), intent(in), dimension(:) :: b real (SP), dimension(size(a,1)*size(b),size(a,2)) :: kro_result integer(I4B) :: i1, i2, j1, ia, ib, ja ia=size(a,1) ib=size(b) ja=size(a,2) do i1=1,ia do i2=1,ib do j1=1,ja kro_result((i1-1)*ib+i2,j1)=a(i1,j1)*b(i2) end do end do end do return end function kromatvec function kro(a,b) result(kro_result) !!!! ! kro=produit de kronecker de a et b !!!! real (SP), intent(in), dimension(:,:) :: a, b real (SP), dimension(size(a,1)*size(b,1),size(a,2)*size(b,2)) :: kro_result integer(I4B) :: i1, i2, j1, j2, ia, ib, ja, jb ia=size(a,1) ib=size(b,1) ja=size(a,2) jb=size(b,2) do i1=1,ia do i2=1,ib do j1=1,ja do j2=1,jb kro_result((i1-1)*ib+i2,(j1-1)*jb+j2)=a(i1,j1)*b(i2,j2) end do end do end do end do return end function kro function trace(a) result(trace_result) !!!! ! trace de a !!!! real (SP), intent(in), dimension(:,:) :: a real (SP) :: trace_result integer(I4B) :: i,n n=size(a,1) if(size(a,2)/=n)then write(unit=4,fmt=*)"erreur dans trace" return end if trace_result=0.0_sp do i=1,n trace_result=trace_result+a(i,i) end do return end function trace function diag(vec) result(diag_result) !!!! ! diag= matrice diagonale de taille nxn dont la premiere diagonale est vec !!!! real (SP), intent(in), dimension(:) :: vec real (SP), dimension(size(vec),size(vec)) :: diag_result integer(I4B) :: i,n n=size(vec) diag_result=0.0_sp do i=1,n diag_result(i,i)=vec(i) end do return end function diag function identity(ia) result(identity_result) !!!! ! identity= matrice identite de taille iaxia !!!! integer(I4B), intent(in) :: ia real (SP), dimension(ia,ia) :: identity_result integer(I4B) :: i identity_result=0.0_sp do i=1,ia identity_result(i,i)=1.0_sp end do return end function identity function norm(a) result(norm_result) !!!! ! norm=somme des valeurs absolues des elements de la matrice a !!!! real (SP), intent(in), dimension(:,:) :: a real (SP) :: norm_result norm_result=sum(abs(a)) return end function norm function normv(a) result(normv_result) !!!! ! norm=somme des valeurs absolues des elements du vecteur a !!!! real (SP), intent(in), dimension(:) :: a real (SP) :: normv_result normv_result=sum(abs(a)) return end function normv end module mod_mat MODULE nr INTERFACE SUBROUTINE airy(x,ai,bi,aip,bip) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: ai,bi,aip,bip END SUBROUTINE airy END INTERFACE INTERFACE SUBROUTINE amebsa(p,y,pb,yb,ftol,func,iter,temptr) USE nrtype INTEGER(I4B), INTENT(INOUT) :: iter REAL(SP), INTENT(INOUT) :: yb REAL(SP), INTENT(IN) :: ftol,temptr REAL(SP), DIMENSION(:), INTENT(INOUT) :: y,pb REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: p INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE amebsa END INTERFACE INTERFACE SUBROUTINE amoeba(p,y,ftol,func,iter) USE nrtype INTEGER(I4B), INTENT(OUT) :: iter REAL(SP), INTENT(IN) :: ftol REAL(SP), DIMENSION(:), INTENT(INOUT) :: y REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: p INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE amoeba END INTERFACE INTERFACE SUBROUTINE anneal(x,y,iorder) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(INOUT) :: iorder REAL(SP), DIMENSION(:), INTENT(IN) :: x,y END SUBROUTINE anneal END INTERFACE INTERFACE SUBROUTINE asolve(b,x,itrnsp) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: b REAL(DP), DIMENSION(:), INTENT(OUT) :: x INTEGER(I4B), INTENT(IN) :: itrnsp END SUBROUTINE asolve END INTERFACE INTERFACE SUBROUTINE atimes(x,r,itrnsp) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: x REAL(DP), DIMENSION(:), INTENT(OUT) :: r INTEGER(I4B), INTENT(IN) :: itrnsp END SUBROUTINE atimes END INTERFACE INTERFACE SUBROUTINE avevar(data,ave,var) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data REAL(SP), INTENT(OUT) :: ave,var END SUBROUTINE avevar END INTERFACE INTERFACE SUBROUTINE balanc(a) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a END SUBROUTINE balanc END INTERFACE INTERFACE SUBROUTINE banbks(a,m1,m2,al,indx,b) USE nrtype INTEGER(I4B), INTENT(IN) :: m1,m2 INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx REAL(SP), DIMENSION(:,:), INTENT(IN) :: a,al REAL(SP), DIMENSION(:), INTENT(INOUT) :: b END SUBROUTINE banbks END INTERFACE INTERFACE SUBROUTINE bandec(a,m1,m2,al,indx,d) USE nrtype INTEGER(I4B), INTENT(IN) :: m1,m2 INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx REAL(SP), INTENT(OUT) :: d REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:,:), INTENT(OUT) :: al END SUBROUTINE bandec END INTERFACE INTERFACE SUBROUTINE banmul(a,m1,m2,x,b) USE nrtype INTEGER(I4B), INTENT(IN) :: m1,m2 REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(OUT) :: b REAL(SP), DIMENSION(:,:), INTENT(IN) :: a END SUBROUTINE banmul END INTERFACE INTERFACE SUBROUTINE bcucof(y,y1,y2,y12,d1,d2,c) USE nrtype REAL(SP), INTENT(IN) :: d1,d2 REAL(SP), DIMENSION(4), INTENT(IN) :: y,y1,y2,y12 REAL(SP), DIMENSION(4,4), INTENT(OUT) :: c END SUBROUTINE bcucof END INTERFACE INTERFACE SUBROUTINE bcuint(y,y1,y2,y12,x1l,x1u,x2l,x2u,x1,x2,ansy,& ansy1,ansy2) USE nrtype REAL(SP), DIMENSION(4), INTENT(IN) :: y,y1,y2,y12 REAL(SP), INTENT(IN) :: x1l,x1u,x2l,x2u,x1,x2 REAL(SP), INTENT(OUT) :: ansy,ansy1,ansy2 END SUBROUTINE bcuint END INTERFACE INTERFACE beschb SUBROUTINE beschb_s(x,gam1,gam2,gampl,gammi) USE nrtype REAL(DP), INTENT(IN) :: x REAL(DP), INTENT(OUT) :: gam1,gam2,gampl,gammi END SUBROUTINE beschb_s !BL SUBROUTINE beschb_v(x,gam1,gam2,gampl,gammi) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: x REAL(DP), DIMENSION(:), INTENT(OUT) :: gam1,gam2,gampl,gammi END SUBROUTINE beschb_v END INTERFACE INTERFACE bessi FUNCTION bessi_s(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: x REAL(SP) :: bessi_s END FUNCTION bessi_s !BL FUNCTION bessi_v(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessi_v END FUNCTION bessi_v END INTERFACE INTERFACE bessi0 FUNCTION bessi0_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessi0_s END FUNCTION bessi0_s !BL FUNCTION bessi0_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessi0_v END FUNCTION bessi0_v END INTERFACE INTERFACE bessi1 FUNCTION bessi1_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessi1_s END FUNCTION bessi1_s !BL FUNCTION bessi1_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessi1_v END FUNCTION bessi1_v END INTERFACE INTERFACE SUBROUTINE bessik(x,xnu,ri,rk,rip,rkp) USE nrtype REAL(SP), INTENT(IN) :: x,xnu REAL(SP), INTENT(OUT) :: ri,rk,rip,rkp END SUBROUTINE bessik END INTERFACE INTERFACE bessj FUNCTION bessj_s(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: x REAL(SP) :: bessj_s END FUNCTION bessj_s !BL FUNCTION bessj_v(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessj_v END FUNCTION bessj_v END INTERFACE INTERFACE bessj0 FUNCTION bessj0_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessj0_s END FUNCTION bessj0_s !BL FUNCTION bessj0_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessj0_v END FUNCTION bessj0_v END INTERFACE INTERFACE bessj1 FUNCTION bessj1_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessj1_s END FUNCTION bessj1_s !BL FUNCTION bessj1_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessj1_v END FUNCTION bessj1_v END INTERFACE INTERFACE bessjy SUBROUTINE bessjy_s(x,xnu,rj,ry,rjp,ryp) USE nrtype REAL(SP), INTENT(IN) :: x,xnu REAL(SP), INTENT(OUT) :: rj,ry,rjp,ryp END SUBROUTINE bessjy_s !BL SUBROUTINE bessjy_v(x,xnu,rj,ry,rjp,ryp) USE nrtype REAL(SP), INTENT(IN) :: xnu REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(OUT) :: rj,rjp,ry,ryp END SUBROUTINE bessjy_v END INTERFACE INTERFACE bessk FUNCTION bessk_s(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: x REAL(SP) :: bessk_s END FUNCTION bessk_s !BL FUNCTION bessk_v(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessk_v END FUNCTION bessk_v END INTERFACE INTERFACE bessk0 FUNCTION bessk0_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessk0_s END FUNCTION bessk0_s !BL FUNCTION bessk0_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessk0_v END FUNCTION bessk0_v END INTERFACE INTERFACE bessk1 FUNCTION bessk1_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessk1_s END FUNCTION bessk1_s !BL FUNCTION bessk1_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessk1_v END FUNCTION bessk1_v END INTERFACE INTERFACE bessy FUNCTION bessy_s(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: x REAL(SP) :: bessy_s END FUNCTION bessy_s !BL FUNCTION bessy_v(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessy_v END FUNCTION bessy_v END INTERFACE INTERFACE bessy0 FUNCTION bessy0_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessy0_s END FUNCTION bessy0_s !BL FUNCTION bessy0_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessy0_v END FUNCTION bessy0_v END INTERFACE INTERFACE bessy1 FUNCTION bessy1_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: bessy1_s END FUNCTION bessy1_s !BL FUNCTION bessy1_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: bessy1_v END FUNCTION bessy1_v END INTERFACE INTERFACE beta FUNCTION beta_s(z,w) USE nrtype REAL(SP), INTENT(IN) :: z,w REAL(SP) :: beta_s END FUNCTION beta_s !BL FUNCTION beta_v(z,w) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: z,w REAL(SP), DIMENSION(size(z)) :: beta_v END FUNCTION beta_v END INTERFACE INTERFACE betacf FUNCTION betacf_s(a,b,x) USE nrtype REAL(SP), INTENT(IN) :: a,b,x REAL(SP) :: betacf_s END FUNCTION betacf_s !BL FUNCTION betacf_v(a,b,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,x REAL(SP), DIMENSION(size(x)) :: betacf_v END FUNCTION betacf_v END INTERFACE INTERFACE betai FUNCTION betai_s(a,b,x) USE nrtype REAL(SP), INTENT(IN) :: a,b,x REAL(SP) :: betai_s END FUNCTION betai_s !BL FUNCTION betai_v(a,b,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,x REAL(SP), DIMENSION(size(a)) :: betai_v END FUNCTION betai_v END INTERFACE INTERFACE bico FUNCTION bico_s(n,k) USE nrtype INTEGER(I4B), INTENT(IN) :: n,k REAL(SP) :: bico_s END FUNCTION bico_s !BL FUNCTION bico_v(n,k) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(IN) :: n,k REAL(SP), DIMENSION(size(n)) :: bico_v END FUNCTION bico_v END INTERFACE INTERFACE FUNCTION bnldev(pp,n) USE nrtype REAL(SP), INTENT(IN) :: pp INTEGER(I4B), INTENT(IN) :: n REAL(SP) :: bnldev END FUNCTION bnldev END INTERFACE INTERFACE FUNCTION brent(ax,bx,cx,func,tol,xmin) USE nrtype REAL(SP), INTENT(IN) :: ax,bx,cx,tol REAL(SP), INTENT(OUT) :: xmin REAL(SP) :: brent INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION brent END INTERFACE INTERFACE SUBROUTINE broydn(x,check) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: x LOGICAL(LGT), INTENT(OUT) :: check END SUBROUTINE broydn END INTERFACE INTERFACE SUBROUTINE bsstep(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal REAL(SP), INTENT(INOUT) :: x REAL(SP), INTENT(IN) :: htry,eps REAL(SP), INTENT(OUT) :: hdid,hnext INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE bsstep END INTERFACE INTERFACE SUBROUTINE caldat(julian,mm,id,iyyy) USE nrtype INTEGER(I4B), INTENT(IN) :: julian INTEGER(I4B), INTENT(OUT) :: mm,id,iyyy END SUBROUTINE caldat END INTERFACE INTERFACE FUNCTION chder(a,b,c) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), DIMENSION(:), INTENT(IN) :: c REAL(SP), DIMENSION(size(c)) :: chder END FUNCTION chder END INTERFACE INTERFACE chebev FUNCTION chebev_s(a,b,c,x) USE nrtype REAL(SP), INTENT(IN) :: a,b,x REAL(SP), DIMENSION(:), INTENT(IN) :: c REAL(SP) :: chebev_s END FUNCTION chebev_s !BL FUNCTION chebev_v(a,b,c,x) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), DIMENSION(:), INTENT(IN) :: c,x REAL(SP), DIMENSION(size(x)) :: chebev_v END FUNCTION chebev_v END INTERFACE INTERFACE FUNCTION chebft(a,b,n,func) USE nrtype REAL(SP), INTENT(IN) :: a,b INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(n) :: chebft INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END FUNCTION chebft END INTERFACE INTERFACE FUNCTION chebpc(c) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: c REAL(SP), DIMENSION(size(c)) :: chebpc END FUNCTION chebpc END INTERFACE INTERFACE FUNCTION chint(a,b,c) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), DIMENSION(:), INTENT(IN) :: c REAL(SP), DIMENSION(size(c)) :: chint END FUNCTION chint END INTERFACE INTERFACE SUBROUTINE choldc(a,p) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:), INTENT(OUT) :: p END SUBROUTINE choldc END INTERFACE INTERFACE SUBROUTINE cholsl(a,p,b,x) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: a REAL(SP), DIMENSION(:), INTENT(IN) :: p,b REAL(SP), DIMENSION(:), INTENT(INOUT) :: x END SUBROUTINE cholsl END INTERFACE INTERFACE SUBROUTINE chsone(bins,ebins,knstrn,df,chsq,prob) USE nrtype INTEGER(I4B), INTENT(IN) :: knstrn REAL(SP), INTENT(OUT) :: df,chsq,prob REAL(SP), DIMENSION(:), INTENT(IN) :: bins,ebins END SUBROUTINE chsone END INTERFACE INTERFACE SUBROUTINE chstwo(bins1,bins2,knstrn,df,chsq,prob) USE nrtype INTEGER(I4B), INTENT(IN) :: knstrn REAL(SP), INTENT(OUT) :: df,chsq,prob REAL(SP), DIMENSION(:), INTENT(IN) :: bins1,bins2 END SUBROUTINE chstwo END INTERFACE INTERFACE SUBROUTINE cisi(x,ci,si) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: ci,si END SUBROUTINE cisi END INTERFACE INTERFACE SUBROUTINE cntab1(nn,chisq,df,prob,cramrv,ccc) USE nrtype INTEGER(I4B), DIMENSION(:,:), INTENT(IN) :: nn REAL(SP), INTENT(OUT) :: chisq,df,prob,cramrv,ccc END SUBROUTINE cntab1 END INTERFACE INTERFACE SUBROUTINE cntab2(nn,h,hx,hy,hygx,hxgy,uygx,uxgy,uxy) USE nrtype INTEGER(I4B), DIMENSION(:,:), INTENT(IN) :: nn REAL(SP), INTENT(OUT) :: h,hx,hy,hygx,hxgy,uygx,uxgy,uxy END SUBROUTINE cntab2 END INTERFACE INTERFACE FUNCTION convlv(data,respns,isign) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data REAL(SP), DIMENSION(:), INTENT(IN) :: respns INTEGER(I4B), INTENT(IN) :: isign REAL(SP), DIMENSION(size(data)) :: convlv END FUNCTION convlv END INTERFACE INTERFACE FUNCTION correl(data1,data2) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 REAL(SP), DIMENSION(size(data1)) :: correl END FUNCTION correl END INTERFACE INTERFACE SUBROUTINE cosft1(y) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y END SUBROUTINE cosft1 END INTERFACE INTERFACE SUBROUTINE cosft2(y,isign) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE cosft2 END INTERFACE INTERFACE SUBROUTINE covsrt(covar,maska) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: covar LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska END SUBROUTINE covsrt END INTERFACE INTERFACE SUBROUTINE cyclic(a,b,c,alpha,beta,r,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN):: a,b,c,r REAL(SP), INTENT(IN) :: alpha,beta REAL(SP), DIMENSION(:), INTENT(OUT):: x END SUBROUTINE cyclic END INTERFACE INTERFACE SUBROUTINE daub4(a,isign) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: a INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE daub4 END INTERFACE INTERFACE dawson FUNCTION dawson_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: dawson_s END FUNCTION dawson_s !BL FUNCTION dawson_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: dawson_v END FUNCTION dawson_v END INTERFACE INTERFACE FUNCTION dbrent(ax,bx,cx,func,dbrent_dfunc,tol,xmin) USE nrtype REAL(SP), INTENT(IN) :: ax,bx,cx,tol REAL(SP), INTENT(OUT) :: xmin REAL(SP) :: dbrent INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func !BL FUNCTION dbrent_dfunc(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: dbrent_dfunc END FUNCTION dbrent_dfunc END INTERFACE END FUNCTION dbrent END INTERFACE INTERFACE SUBROUTINE ddpoly(c,x,pd) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: c REAL(SP), DIMENSION(:), INTENT(OUT) :: pd END SUBROUTINE ddpoly END INTERFACE INTERFACE FUNCTION decchk(string,ch) USE nrtype CHARACTER(1), DIMENSION(:), INTENT(IN) :: string CHARACTER(1), INTENT(OUT) :: ch LOGICAL(LGT) :: decchk END FUNCTION decchk END INTERFACE INTERFACE SUBROUTINE dfpmin(p,gtol,iter,fret,func,dfunc) USE nrtype INTEGER(I4B), INTENT(OUT) :: iter REAL(SP), INTENT(IN) :: gtol REAL(SP), INTENT(OUT) :: fret REAL(SP), DIMENSION(:), INTENT(INOUT) :: p INTERFACE FUNCTION func(p) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: p REAL(SP) :: func END FUNCTION func !BL FUNCTION dfunc(p) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: p REAL(SP), DIMENSION(size(p)) :: dfunc END FUNCTION dfunc END INTERFACE END SUBROUTINE dfpmin END INTERFACE INTERFACE FUNCTION dfridr(func,x,h,err) USE nrtype REAL(SP), INTENT(IN) :: x,h REAL(SP), INTENT(OUT) :: err REAL(SP) :: dfridr INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION dfridr END INTERFACE INTERFACE SUBROUTINE dftcor(w,delta,a,b,endpts,corre,corim,corfac) USE nrtype REAL(SP), INTENT(IN) :: w,delta,a,b REAL(SP), INTENT(OUT) :: corre,corim,corfac REAL(SP), DIMENSION(:), INTENT(IN) :: endpts END SUBROUTINE dftcor END INTERFACE INTERFACE SUBROUTINE dftint(func,a,b,w,cosint,sinint) USE nrtype REAL(SP), INTENT(IN) :: a,b,w REAL(SP), INTENT(OUT) :: cosint,sinint INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END SUBROUTINE dftint END INTERFACE INTERFACE SUBROUTINE difeq(k,k1,k2,jsf,is1,isf,indexv,s,y) USE nrtype INTEGER(I4B), INTENT(IN) :: is1,isf,jsf,k,k1,k2 INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indexv REAL(SP), DIMENSION(:,:), INTENT(OUT) :: s REAL(SP), DIMENSION(:,:), INTENT(IN) :: y END SUBROUTINE difeq END INTERFACE INTERFACE FUNCTION eclass(lista,listb,n) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(IN) :: lista,listb INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B), DIMENSION(n) :: eclass END FUNCTION eclass END INTERFACE INTERFACE FUNCTION eclazz(equiv,n) USE nrtype INTERFACE FUNCTION equiv(i,j) USE nrtype LOGICAL(LGT) :: equiv INTEGER(I4B), INTENT(IN) :: i,j END FUNCTION equiv END INTERFACE INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B), DIMENSION(n) :: eclazz END FUNCTION eclazz END INTERFACE INTERFACE FUNCTION ei(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: ei END FUNCTION ei END INTERFACE INTERFACE SUBROUTINE eigsrt(d,v) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: d REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: v END SUBROUTINE eigsrt END INTERFACE INTERFACE elle FUNCTION elle_s(phi,ak) USE nrtype REAL(SP), INTENT(IN) :: phi,ak REAL(SP) :: elle_s END FUNCTION elle_s !BL FUNCTION elle_v(phi,ak) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: phi,ak REAL(SP), DIMENSION(size(phi)) :: elle_v END FUNCTION elle_v END INTERFACE INTERFACE ellf FUNCTION ellf_s(phi,ak) USE nrtype REAL(SP), INTENT(IN) :: phi,ak REAL(SP) :: ellf_s END FUNCTION ellf_s !BL FUNCTION ellf_v(phi,ak) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: phi,ak REAL(SP), DIMENSION(size(phi)) :: ellf_v END FUNCTION ellf_v END INTERFACE INTERFACE ellpi FUNCTION ellpi_s(phi,en,ak) USE nrtype REAL(SP), INTENT(IN) :: phi,en,ak REAL(SP) :: ellpi_s END FUNCTION ellpi_s !BL FUNCTION ellpi_v(phi,en,ak) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: phi,en,ak REAL(SP), DIMENSION(size(phi)) :: ellpi_v END FUNCTION ellpi_v END INTERFACE INTERFACE SUBROUTINE elmhes(a) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a END SUBROUTINE elmhes END INTERFACE INTERFACE erf FUNCTION erf_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: erf_s END FUNCTION erf_s !BL FUNCTION erf_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: erf_v END FUNCTION erf_v END INTERFACE INTERFACE erfc FUNCTION erfc_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: erfc_s END FUNCTION erfc_s !BL FUNCTION erfc_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: erfc_v END FUNCTION erfc_v END INTERFACE INTERFACE erfcc FUNCTION erfcc_s(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: erfcc_s END FUNCTION erfcc_s !BL FUNCTION erfcc_v(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: erfcc_v END FUNCTION erfcc_v END INTERFACE INTERFACE SUBROUTINE eulsum(sum,term,jterm) USE nrtype REAL(SP), INTENT(INOUT) :: sum REAL(SP), INTENT(IN) :: term INTEGER(I4B), INTENT(IN) :: jterm END SUBROUTINE eulsum END INTERFACE INTERFACE FUNCTION evlmem(fdt,d,xms) USE nrtype REAL(SP), INTENT(IN) :: fdt,xms REAL(SP), DIMENSION(:), INTENT(IN) :: d REAL(SP) :: evlmem END FUNCTION evlmem END INTERFACE INTERFACE expdev SUBROUTINE expdev_s(harvest) USE nrtype REAL(SP), INTENT(OUT) :: harvest END SUBROUTINE expdev_s !BL SUBROUTINE expdev_v(harvest) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest END SUBROUTINE expdev_v END INTERFACE INTERFACE FUNCTION expint(n,x) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: x REAL(SP) :: expint END FUNCTION expint END INTERFACE INTERFACE factln FUNCTION factln_s(n) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP) :: factln_s END FUNCTION factln_s !BL FUNCTION factln_v(n) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(IN) :: n REAL(SP), DIMENSION(size(n)) :: factln_v END FUNCTION factln_v END INTERFACE INTERFACE factrl FUNCTION factrl_s(n) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP) :: factrl_s END FUNCTION factrl_s !BL FUNCTION factrl_v(n) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(IN) :: n REAL(SP), DIMENSION(size(n)) :: factrl_v END FUNCTION factrl_v END INTERFACE INTERFACE SUBROUTINE fasper(x,y,ofac,hifac,px,py,jmax,prob) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), INTENT(IN) :: ofac,hifac INTEGER(I4B), INTENT(OUT) :: jmax REAL(SP), INTENT(OUT) :: prob REAL(SP), DIMENSION(:), POINTER :: px,py END SUBROUTINE fasper END INTERFACE INTERFACE SUBROUTINE fdjac(x,fvec,df) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: fvec REAL(SP), DIMENSION(:), INTENT(INOUT) :: x REAL(SP), DIMENSION(:,:), INTENT(OUT) :: df END SUBROUTINE fdjac END INTERFACE INTERFACE SUBROUTINE fgauss(x,a,y,dyda) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,a REAL(SP), DIMENSION(:), INTENT(OUT) :: y REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dyda END SUBROUTINE fgauss END INTERFACE INTERFACE SUBROUTINE fit(x,y,a,b,siga,sigb,chi2,q,sig) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), INTENT(OUT) :: a,b,siga,sigb,chi2,q REAL(SP), DIMENSION(:), OPTIONAL, INTENT(IN) :: sig END SUBROUTINE fit END INTERFACE INTERFACE SUBROUTINE fitexy(x,y,sigx,sigy,a,b,siga,sigb,chi2,q) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sigx,sigy REAL(SP), INTENT(OUT) :: a,b,siga,sigb,chi2,q END SUBROUTINE fitexy END INTERFACE INTERFACE SUBROUTINE fixrts(d) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: d END SUBROUTINE fixrts END INTERFACE INTERFACE FUNCTION fleg(x,n) USE nrtype REAL(SP), INTENT(IN) :: x INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(n) :: fleg END FUNCTION fleg END INTERFACE INTERFACE SUBROUTINE flmoon(n,nph,jd,frac) USE nrtype INTEGER(I4B), INTENT(IN) :: n,nph INTEGER(I4B), INTENT(OUT) :: jd REAL(SP), INTENT(OUT) :: frac END SUBROUTINE flmoon END INTERFACE INTERFACE four1 SUBROUTINE four1_dp(data,isign) USE nrtype COMPLEX(DPC), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE four1_dp !BL SUBROUTINE four1_sp(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE four1_sp END INTERFACE INTERFACE SUBROUTINE four1_alt(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE four1_alt END INTERFACE INTERFACE SUBROUTINE four1_gather(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE four1_gather END INTERFACE INTERFACE SUBROUTINE four2(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data INTEGER(I4B),INTENT(IN) :: isign END SUBROUTINE four2 END INTERFACE INTERFACE SUBROUTINE four2_alt(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE four2_alt END INTERFACE INTERFACE SUBROUTINE four3(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data INTEGER(I4B),INTENT(IN) :: isign END SUBROUTINE four3 END INTERFACE INTERFACE SUBROUTINE four3_alt(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE four3_alt END INTERFACE INTERFACE SUBROUTINE fourcol(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE fourcol END INTERFACE INTERFACE SUBROUTINE fourcol_3d(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE fourcol_3d END INTERFACE INTERFACE SUBROUTINE fourn_gather(data,nn,isign) USE nrtype COMPLEX(SPC), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), DIMENSION(:), INTENT(IN) :: nn INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE fourn_gather END INTERFACE INTERFACE fourrow SUBROUTINE fourrow_dp(data,isign) USE nrtype COMPLEX(DPC), DIMENSION(:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE fourrow_dp !BL SUBROUTINE fourrow_sp(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE fourrow_sp END INTERFACE INTERFACE SUBROUTINE fourrow_3d(data,isign) USE nrtype COMPLEX(SPC), DIMENSION(:,:,:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE fourrow_3d END INTERFACE INTERFACE FUNCTION fpoly(x,n) USE nrtype REAL(SP), INTENT(IN) :: x INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(n) :: fpoly END FUNCTION fpoly END INTERFACE INTERFACE SUBROUTINE fred2(a,b,t,f,w,g,ak) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), DIMENSION(:), INTENT(OUT) :: t,f,w INTERFACE FUNCTION g(t) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: t REAL(SP), DIMENSION(size(t)) :: g END FUNCTION g !BL FUNCTION ak(t,s) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: t,s REAL(SP), DIMENSION(size(t),size(s)) :: ak END FUNCTION ak END INTERFACE END SUBROUTINE fred2 END INTERFACE INTERFACE FUNCTION fredin(x,a,b,t,f,w,g,ak) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), DIMENSION(:), INTENT(IN) :: x,t,f,w REAL(SP), DIMENSION(size(x)) :: fredin INTERFACE FUNCTION g(t) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: t REAL(SP), DIMENSION(size(t)) :: g END FUNCTION g !BL FUNCTION ak(t,s) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: t,s REAL(SP), DIMENSION(size(t),size(s)) :: ak END FUNCTION ak END INTERFACE END FUNCTION fredin END INTERFACE INTERFACE SUBROUTINE frenel(x,s,c) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: s,c END SUBROUTINE frenel END INTERFACE INTERFACE SUBROUTINE frprmn(p,ftol,iter,fret) USE nrtype INTEGER(I4B), INTENT(OUT) :: iter REAL(SP), INTENT(IN) :: ftol REAL(SP), INTENT(OUT) :: fret REAL(SP), DIMENSION(:), INTENT(INOUT) :: p END SUBROUTINE frprmn END INTERFACE INTERFACE SUBROUTINE ftest(data1,data2,f,prob) USE nrtype REAL(SP), INTENT(OUT) :: f,prob REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 END SUBROUTINE ftest END INTERFACE INTERFACE FUNCTION gamdev(ia) USE nrtype INTEGER(I4B), INTENT(IN) :: ia REAL(SP) :: gamdev END FUNCTION gamdev END INTERFACE INTERFACE gammln FUNCTION gammln_s(xx) USE nrtype REAL(SP), INTENT(IN) :: xx REAL(SP) :: gammln_s END FUNCTION gammln_s !BL FUNCTION gammln_v(xx) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xx REAL(SP), DIMENSION(size(xx)) :: gammln_v END FUNCTION gammln_v END INTERFACE INTERFACE gammp FUNCTION gammp_s(a,x) USE nrtype REAL(SP), INTENT(IN) :: a,x REAL(SP) :: gammp_s END FUNCTION gammp_s !BL FUNCTION gammp_v(a,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,x REAL(SP), DIMENSION(size(a)) :: gammp_v END FUNCTION gammp_v END INTERFACE INTERFACE gammq FUNCTION gammq_s(a,x) USE nrtype REAL(SP), INTENT(IN) :: a,x REAL(SP) :: gammq_s END FUNCTION gammq_s !BL FUNCTION gammq_v(a,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,x REAL(SP), DIMENSION(size(a)) :: gammq_v END FUNCTION gammq_v END INTERFACE INTERFACE gasdev SUBROUTINE gasdev_s(harvest) USE nrtype REAL(SP), INTENT(OUT) :: harvest END SUBROUTINE gasdev_s !BL SUBROUTINE gasdev_v(harvest) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest END SUBROUTINE gasdev_v END INTERFACE INTERFACE SUBROUTINE gaucof(a,b,amu0,x,w) USE nrtype REAL(SP), INTENT(IN) :: amu0 REAL(SP), DIMENSION(:), INTENT(INOUT) :: a,b REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w END SUBROUTINE gaucof END INTERFACE INTERFACE SUBROUTINE gauher(x,w) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w END SUBROUTINE gauher END INTERFACE INTERFACE SUBROUTINE gaujac(x,w,alf,bet) USE nrtype REAL(SP), INTENT(IN) :: alf,bet REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w END SUBROUTINE gaujac END INTERFACE INTERFACE SUBROUTINE gaulag(x,w,alf) USE nrtype REAL(SP), INTENT(IN) :: alf REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w END SUBROUTINE gaulag END INTERFACE INTERFACE SUBROUTINE gauleg(x1,x2,x,w) USE nrtype REAL(SP), INTENT(IN) :: x1,x2 REAL(SP), DIMENSION(:), INTENT(OUT) :: x,w END SUBROUTINE gauleg END INTERFACE INTERFACE SUBROUTINE gaussj(a,b) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a,b END SUBROUTINE gaussj END INTERFACE INTERFACE gcf FUNCTION gcf_s(a,x,gln) USE nrtype REAL(SP), INTENT(IN) :: a,x REAL(SP), OPTIONAL, INTENT(OUT) :: gln REAL(SP) :: gcf_s END FUNCTION gcf_s !BL FUNCTION gcf_v(a,x,gln) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,x REAL(SP), DIMENSION(:), OPTIONAL, INTENT(OUT) :: gln REAL(SP), DIMENSION(size(a)) :: gcf_v END FUNCTION gcf_v END INTERFACE INTERFACE FUNCTION golden(ax,bx,cx,func,tol,xmin) USE nrtype REAL(SP), INTENT(IN) :: ax,bx,cx,tol REAL(SP), INTENT(OUT) :: xmin REAL(SP) :: golden INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION golden END INTERFACE INTERFACE gser FUNCTION gser_s(a,x,gln) USE nrtype REAL(SP), INTENT(IN) :: a,x REAL(SP), OPTIONAL, INTENT(OUT) :: gln REAL(SP) :: gser_s END FUNCTION gser_s !BL FUNCTION gser_v(a,x,gln) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,x REAL(SP), DIMENSION(:), OPTIONAL, INTENT(OUT) :: gln REAL(SP), DIMENSION(size(a)) :: gser_v END FUNCTION gser_v END INTERFACE INTERFACE SUBROUTINE hqr(a,wr,wi) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: wr,wi REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a END SUBROUTINE hqr END INTERFACE INTERFACE SUBROUTINE hunt(xx,x,jlo) USE nrtype INTEGER(I4B), INTENT(INOUT) :: jlo REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: xx END SUBROUTINE hunt END INTERFACE INTERFACE SUBROUTINE hypdrv(s,ry,rdyds) USE nrtype REAL(SP), INTENT(IN) :: s REAL(SP), DIMENSION(:), INTENT(IN) :: ry REAL(SP), DIMENSION(:), INTENT(OUT) :: rdyds END SUBROUTINE hypdrv END INTERFACE INTERFACE FUNCTION hypgeo(a,b,c,z) USE nrtype COMPLEX(SPC), INTENT(IN) :: a,b,c,z COMPLEX(SPC) :: hypgeo END FUNCTION hypgeo END INTERFACE INTERFACE SUBROUTINE hypser(a,b,c,z,series,deriv) USE nrtype COMPLEX(SPC), INTENT(IN) :: a,b,c,z COMPLEX(SPC), INTENT(OUT) :: series,deriv END SUBROUTINE hypser END INTERFACE INTERFACE FUNCTION icrc(crc,buf,jinit,jrev) USE nrtype CHARACTER(1), DIMENSION(:), INTENT(IN) :: buf INTEGER(I2B), INTENT(IN) :: crc,jinit INTEGER(I4B), INTENT(IN) :: jrev INTEGER(I2B) :: icrc END FUNCTION icrc END INTERFACE INTERFACE FUNCTION igray(n,is) USE nrtype INTEGER(I4B), INTENT(IN) :: n,is INTEGER(I4B) :: igray END FUNCTION igray END INTERFACE INTERFACE RECURSIVE SUBROUTINE index_bypack(arr,index,partial) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), DIMENSION(:), INTENT(INOUT) :: index INTEGER, OPTIONAL, INTENT(IN) :: partial END SUBROUTINE index_bypack END INTERFACE INTERFACE indexx SUBROUTINE indexx_sp(arr,index) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index END SUBROUTINE indexx_sp SUBROUTINE indexx_i4b(iarr,index) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index END SUBROUTINE indexx_i4b END INTERFACE INTERFACE FUNCTION interp(uc) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: uc REAL(DP), DIMENSION(2*size(uc,1)-1,2*size(uc,1)-1) :: interp END FUNCTION interp END INTERFACE INTERFACE FUNCTION rank(indx) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx INTEGER(I4B), DIMENSION(size(indx)) :: rank END FUNCTION rank END INTERFACE INTERFACE FUNCTION irbit1(iseed) USE nrtype INTEGER(I4B), INTENT(INOUT) :: iseed INTEGER(I4B) :: irbit1 END FUNCTION irbit1 END INTERFACE INTERFACE FUNCTION irbit2(iseed) USE nrtype INTEGER(I4B), INTENT(INOUT) :: iseed INTEGER(I4B) :: irbit2 END FUNCTION irbit2 END INTERFACE INTERFACE SUBROUTINE jacobi(a,d,v,nrot) USE nrtype INTEGER(I4B), INTENT(OUT) :: nrot REAL(SP), DIMENSION(:), INTENT(OUT) :: d REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v END SUBROUTINE jacobi END INTERFACE INTERFACE SUBROUTINE jacobn(x,y,dfdx,dfdy) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dfdx REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dfdy END SUBROUTINE jacobn END INTERFACE INTERFACE FUNCTION julday(mm,id,iyyy) USE nrtype INTEGER(I4B), INTENT(IN) :: mm,id,iyyy INTEGER(I4B) :: julday END FUNCTION julday END INTERFACE INTERFACE SUBROUTINE kendl1(data1,data2,tau,z,prob) USE nrtype REAL(SP), INTENT(OUT) :: tau,z,prob REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 END SUBROUTINE kendl1 END INTERFACE INTERFACE SUBROUTINE kendl2(tab,tau,z,prob) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: tab REAL(SP), INTENT(OUT) :: tau,z,prob END SUBROUTINE kendl2 END INTERFACE INTERFACE FUNCTION kermom(y,m) USE nrtype REAL(DP), INTENT(IN) :: y INTEGER(I4B), INTENT(IN) :: m REAL(DP), DIMENSION(m) :: kermom END FUNCTION kermom END INTERFACE INTERFACE SUBROUTINE ks2d1s(x1,y1,quadvl,d1,prob) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x1,y1 REAL(SP), INTENT(OUT) :: d1,prob INTERFACE SUBROUTINE quadvl(x,y,fa,fb,fc,fd) USE nrtype REAL(SP), INTENT(IN) :: x,y REAL(SP), INTENT(OUT) :: fa,fb,fc,fd END SUBROUTINE quadvl END INTERFACE END SUBROUTINE ks2d1s END INTERFACE INTERFACE SUBROUTINE ks2d2s(x1,y1,x2,y2,d,prob) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x1,y1,x2,y2 REAL(SP), INTENT(OUT) :: d,prob END SUBROUTINE ks2d2s END INTERFACE INTERFACE SUBROUTINE ksone(data,func,d,prob) USE nrtype REAL(SP), INTENT(OUT) :: d,prob REAL(SP), DIMENSION(:), INTENT(INOUT) :: data INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END SUBROUTINE ksone END INTERFACE INTERFACE SUBROUTINE kstwo(data1,data2,d,prob) USE nrtype REAL(SP), INTENT(OUT) :: d,prob REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 END SUBROUTINE kstwo END INTERFACE INTERFACE SUBROUTINE laguer(a,x,its) USE nrtype INTEGER(I4B), INTENT(OUT) :: its COMPLEX(SPC), INTENT(INOUT) :: x COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a END SUBROUTINE laguer END INTERFACE INTERFACE SUBROUTINE lfit(x,y,sig,a,maska,covar,chisq,funcs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig REAL(SP), DIMENSION(:), INTENT(INOUT) :: a LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: covar REAL(SP), INTENT(OUT) :: chisq INTERFACE SUBROUTINE funcs(x,arr) USE nrtype REAL(SP),INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(OUT) :: arr END SUBROUTINE funcs END INTERFACE END SUBROUTINE lfit END INTERFACE INTERFACE SUBROUTINE linbcg(b,x,itol,tol,itmax,iter,err) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: b REAL(DP), DIMENSION(:), INTENT(INOUT) :: x INTEGER(I4B), INTENT(IN) :: itol,itmax REAL(DP), INTENT(IN) :: tol INTEGER(I4B), INTENT(OUT) :: iter REAL(DP), INTENT(OUT) :: err END SUBROUTINE linbcg END INTERFACE INTERFACE SUBROUTINE linmin(p,xi,fret) USE nrtype REAL(SP), INTENT(OUT) :: fret REAL(SP), DIMENSION(:), TARGET, INTENT(INOUT) :: p,xi END SUBROUTINE linmin END INTERFACE INTERFACE SUBROUTINE lnsrch(xold,fold,g,p,x,f,stpmax,check,func) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xold,g REAL(SP), DIMENSION(:), INTENT(INOUT) :: p REAL(SP), INTENT(IN) :: fold,stpmax REAL(SP), DIMENSION(:), INTENT(OUT) :: x REAL(SP), INTENT(OUT) :: f LOGICAL(LGT), INTENT(OUT) :: check INTERFACE FUNCTION func(x) USE nrtype REAL(SP) :: func REAL(SP), DIMENSION(:), INTENT(IN) :: x END FUNCTION func END INTERFACE END SUBROUTINE lnsrch END INTERFACE INTERFACE FUNCTION locate(xx,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xx REAL(SP), INTENT(IN) :: x INTEGER(I4B) :: locate END FUNCTION locate END INTERFACE INTERFACE FUNCTION lop(u) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: u REAL(DP), DIMENSION(size(u,1),size(u,1)) :: lop END FUNCTION lop END INTERFACE INTERFACE SUBROUTINE lubksb(a,indx,b) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: a INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx REAL(SP), DIMENSION(:), INTENT(INOUT) :: b END SUBROUTINE lubksb END INTERFACE INTERFACE SUBROUTINE ludcmp(a,indx,d) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: indx REAL(SP), INTENT(OUT) :: d END SUBROUTINE ludcmp END INTERFACE INTERFACE SUBROUTINE machar(ibeta,it,irnd,ngrd,machep,negep,iexp,minexp,& maxexp,eps,epsneg,xmin,xmax) USE nrtype INTEGER(I4B), INTENT(OUT) :: ibeta,iexp,irnd,it,machep,maxexp,& minexp,negep,ngrd REAL(SP), INTENT(OUT) :: eps,epsneg,xmax,xmin END SUBROUTINE machar END INTERFACE INTERFACE SUBROUTINE medfit(x,y,a,b,abdev) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), INTENT(OUT) :: a,b,abdev END SUBROUTINE medfit END INTERFACE INTERFACE SUBROUTINE memcof(data,xms,d) USE nrtype REAL(SP), INTENT(OUT) :: xms REAL(SP), DIMENSION(:), INTENT(IN) :: data REAL(SP), DIMENSION(:), INTENT(OUT) :: d END SUBROUTINE memcof END INTERFACE INTERFACE SUBROUTINE mgfas(u,maxcyc) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u INTEGER(I4B), INTENT(IN) :: maxcyc END SUBROUTINE mgfas END INTERFACE INTERFACE SUBROUTINE mglin(u,ncycle) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u INTEGER(I4B), INTENT(IN) :: ncycle END SUBROUTINE mglin END INTERFACE INTERFACE SUBROUTINE midexp(funk,aa,bb,s,n) USE nrtype REAL(SP), INTENT(IN) :: aa,bb REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION funk(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: funk END FUNCTION funk END INTERFACE END SUBROUTINE midexp END INTERFACE INTERFACE SUBROUTINE midinf(funk,aa,bb,s,n) USE nrtype REAL(SP), INTENT(IN) :: aa,bb REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION funk(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: funk END FUNCTION funk END INTERFACE END SUBROUTINE midinf END INTERFACE INTERFACE SUBROUTINE midpnt(func,a,b,s,n) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END SUBROUTINE midpnt END INTERFACE INTERFACE SUBROUTINE midsql(funk,aa,bb,s,n) USE nrtype REAL(SP), INTENT(IN) :: aa,bb REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION funk(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: funk END FUNCTION funk END INTERFACE END SUBROUTINE midsql END INTERFACE INTERFACE SUBROUTINE midsqu(funk,aa,bb,s,n) USE nrtype REAL(SP), INTENT(IN) :: aa,bb REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION funk(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: funk END FUNCTION funk END INTERFACE END SUBROUTINE midsqu END INTERFACE INTERFACE RECURSIVE SUBROUTINE miser(func,regn,ndim,npts,dith,ave,var) USE nrtype INTERFACE FUNCTION func(x) USE nrtype REAL(SP) :: func REAL(SP), DIMENSION(:), INTENT(IN) :: x END FUNCTION func END INTERFACE REAL(SP), DIMENSION(:), INTENT(IN) :: regn INTEGER(I4B), INTENT(IN) :: ndim,npts REAL(SP), INTENT(IN) :: dith REAL(SP), INTENT(OUT) :: ave,var END SUBROUTINE miser END INTERFACE INTERFACE SUBROUTINE mmid(y,dydx,xs,htot,nstep,yout,derivs) USE nrtype INTEGER(I4B), INTENT(IN) :: nstep REAL(SP), INTENT(IN) :: xs,htot REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx REAL(SP), DIMENSION(:), INTENT(OUT) :: yout INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE mmid END INTERFACE INTERFACE SUBROUTINE mnbrak(ax,bx,cx,fa,fb,fc,func) USE nrtype REAL(SP), INTENT(INOUT) :: ax,bx REAL(SP), INTENT(OUT) :: cx,fa,fb,fc INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE mnbrak END INTERFACE INTERFACE SUBROUTINE mnewt(ntrial,x,tolx,tolf,usrfun) USE nrtype INTEGER(I4B), INTENT(IN) :: ntrial REAL(SP), INTENT(IN) :: tolx,tolf REAL(SP), DIMENSION(:), INTENT(INOUT) :: x INTERFACE SUBROUTINE usrfun(x,fvec,fjac) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(OUT) :: fvec REAL(SP), DIMENSION(:,:), INTENT(OUT) :: fjac END SUBROUTINE usrfun END INTERFACE END SUBROUTINE mnewt END INTERFACE INTERFACE SUBROUTINE moment(data,ave,adev,sdev,var,skew,curt) USE nrtype REAL(SP), INTENT(OUT) :: ave,adev,sdev,var,skew,curt REAL(SP), DIMENSION(:), INTENT(IN) :: data END SUBROUTINE moment END INTERFACE INTERFACE SUBROUTINE mp2dfr(a,s,n,m) USE nrtype INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B), INTENT(OUT) :: m CHARACTER(1), DIMENSION(:), INTENT(INOUT) :: a CHARACTER(1), DIMENSION(:), INTENT(OUT) :: s END SUBROUTINE mp2dfr END INTERFACE INTERFACE SUBROUTINE mpdiv(q,r,u,v,n,m) USE nrtype CHARACTER(1), DIMENSION(:), INTENT(OUT) :: q,r CHARACTER(1), DIMENSION(:), INTENT(IN) :: u,v INTEGER(I4B), INTENT(IN) :: n,m END SUBROUTINE mpdiv END INTERFACE INTERFACE SUBROUTINE mpinv(u,v,n,m) USE nrtype CHARACTER(1), DIMENSION(:), INTENT(OUT) :: u CHARACTER(1), DIMENSION(:), INTENT(IN) :: v INTEGER(I4B), INTENT(IN) :: n,m END SUBROUTINE mpinv END INTERFACE INTERFACE SUBROUTINE mpmul(w,u,v,n,m) USE nrtype CHARACTER(1), DIMENSION(:), INTENT(IN) :: u,v CHARACTER(1), DIMENSION(:), INTENT(OUT) :: w INTEGER(I4B), INTENT(IN) :: n,m END SUBROUTINE mpmul END INTERFACE INTERFACE SUBROUTINE mppi(n) USE nrtype INTEGER(I4B), INTENT(IN) :: n END SUBROUTINE mppi END INTERFACE INTERFACE SUBROUTINE mprove(a,alud,indx,b,x) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: a,alud INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indx REAL(SP), DIMENSION(:), INTENT(IN) :: b REAL(SP), DIMENSION(:), INTENT(INOUT) :: x END SUBROUTINE mprove END INTERFACE INTERFACE SUBROUTINE mpsqrt(w,u,v,n,m) USE nrtype CHARACTER(1), DIMENSION(:), INTENT(OUT) :: w,u CHARACTER(1), DIMENSION(:), INTENT(IN) :: v INTEGER(I4B), INTENT(IN) :: n,m END SUBROUTINE mpsqrt END INTERFACE INTERFACE SUBROUTINE mrqcof(x,y,sig,a,maska,alpha,beta,chisq,funcs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,a,sig REAL(SP), DIMENSION(:), INTENT(OUT) :: beta REAL(SP), DIMENSION(:,:), INTENT(OUT) :: alpha REAL(SP), INTENT(OUT) :: chisq LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska INTERFACE SUBROUTINE funcs(x,a,yfit,dyda) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,a REAL(SP), DIMENSION(:), INTENT(OUT) :: yfit REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dyda END SUBROUTINE funcs END INTERFACE END SUBROUTINE mrqcof END INTERFACE INTERFACE SUBROUTINE mrqmin(x,y,sig,a,maska,covar,alpha,chisq,funcs,alamda) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig REAL(SP), DIMENSION(:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:,:), INTENT(OUT) :: covar,alpha REAL(SP), INTENT(OUT) :: chisq REAL(SP), INTENT(INOUT) :: alamda LOGICAL(LGT), DIMENSION(:), INTENT(IN) :: maska INTERFACE SUBROUTINE funcs(x,a,yfit,dyda) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,a REAL(SP), DIMENSION(:), INTENT(OUT) :: yfit REAL(SP), DIMENSION(:,:), INTENT(OUT) :: dyda END SUBROUTINE funcs END INTERFACE END SUBROUTINE mrqmin END INTERFACE INTERFACE SUBROUTINE newt(x,check) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: x LOGICAL(LGT), INTENT(OUT) :: check END SUBROUTINE newt END INTERFACE INTERFACE SUBROUTINE odeint(ystart,x1,x2,eps,h1,hmin,derivs,rkqs) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: ystart REAL(SP), INTENT(IN) :: x1,x2,eps,h1,hmin INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs !BL SUBROUTINE rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal REAL(SP), INTENT(INOUT) :: x REAL(SP), INTENT(IN) :: htry,eps REAL(SP), INTENT(OUT) :: hdid,hnext INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE rkqs END INTERFACE END SUBROUTINE odeint END INTERFACE INTERFACE SUBROUTINE orthog(anu,alpha,beta,a,b) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: anu,alpha,beta REAL(SP), DIMENSION(:), INTENT(OUT) :: a,b END SUBROUTINE orthog END INTERFACE INTERFACE SUBROUTINE pade(cof,resid) USE nrtype REAL(DP), DIMENSION(:), INTENT(INOUT) :: cof REAL(SP), INTENT(OUT) :: resid END SUBROUTINE pade END INTERFACE INTERFACE FUNCTION pccheb(d) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: d REAL(SP), DIMENSION(size(d)) :: pccheb END FUNCTION pccheb END INTERFACE INTERFACE SUBROUTINE pcshft(a,b,d) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), DIMENSION(:), INTENT(INOUT) :: d END SUBROUTINE pcshft END INTERFACE INTERFACE SUBROUTINE pearsn(x,y,r,prob,z) USE nrtype REAL(SP), INTENT(OUT) :: r,prob,z REAL(SP), DIMENSION(:), INTENT(IN) :: x,y END SUBROUTINE pearsn END INTERFACE INTERFACE SUBROUTINE period(x,y,ofac,hifac,px,py,jmax,prob) USE nrtype INTEGER(I4B), INTENT(OUT) :: jmax REAL(SP), INTENT(IN) :: ofac,hifac REAL(SP), INTENT(OUT) :: prob REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), DIMENSION(:), POINTER :: px,py END SUBROUTINE period END INTERFACE INTERFACE plgndr FUNCTION plgndr_s(l,m,x) USE nrtype INTEGER(I4B), INTENT(IN) :: l,m REAL(SP), INTENT(IN) :: x REAL(SP) :: plgndr_s END FUNCTION plgndr_s !BL FUNCTION plgndr_v(l,m,x) USE nrtype INTEGER(I4B), INTENT(IN) :: l,m REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: plgndr_v END FUNCTION plgndr_v END INTERFACE INTERFACE FUNCTION poidev(xm) USE nrtype REAL(SP), INTENT(IN) :: xm REAL(SP) :: poidev END FUNCTION poidev END INTERFACE INTERFACE FUNCTION polcoe(x,y) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), DIMENSION(size(x)) :: polcoe END FUNCTION polcoe END INTERFACE INTERFACE FUNCTION polcof(xa,ya) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya REAL(SP), DIMENSION(size(xa)) :: polcof END FUNCTION polcof END INTERFACE INTERFACE SUBROUTINE poldiv(u,v,q,r) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: u,v REAL(SP), DIMENSION(:), INTENT(OUT) :: q,r END SUBROUTINE poldiv END INTERFACE INTERFACE SUBROUTINE polin2(x1a,x2a,ya,x1,x2,y,dy) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x1a,x2a REAL(SP), DIMENSION(:,:), INTENT(IN) :: ya REAL(SP), INTENT(IN) :: x1,x2 REAL(SP), INTENT(OUT) :: y,dy END SUBROUTINE polin2 END INTERFACE INTERFACE SUBROUTINE polint(xa,ya,x,y,dy) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: y,dy END SUBROUTINE polint END INTERFACE INTERFACE SUBROUTINE powell(p,xi,ftol,iter,fret) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: p REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: xi INTEGER(I4B), INTENT(OUT) :: iter REAL(SP), INTENT(IN) :: ftol REAL(SP), INTENT(OUT) :: fret END SUBROUTINE powell END INTERFACE INTERFACE FUNCTION predic(data,d,nfut) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data,d INTEGER(I4B), INTENT(IN) :: nfut REAL(SP), DIMENSION(nfut) :: predic END FUNCTION predic END INTERFACE INTERFACE FUNCTION probks(alam) USE nrtype REAL(SP), INTENT(IN) :: alam REAL(SP) :: probks END FUNCTION probks END INTERFACE INTERFACE psdes SUBROUTINE psdes_s(lword,rword) USE nrtype INTEGER(I4B), INTENT(INOUT) :: lword,rword END SUBROUTINE psdes_s !BL SUBROUTINE psdes_v(lword,rword) USE nrtype INTEGER(I4B), DIMENSION(:), INTENT(INOUT) :: lword,rword END SUBROUTINE psdes_v END INTERFACE INTERFACE SUBROUTINE pwt(a,isign) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: a INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE pwt END INTERFACE INTERFACE SUBROUTINE pwtset(n) USE nrtype INTEGER(I4B), INTENT(IN) :: n END SUBROUTINE pwtset END INTERFACE INTERFACE pythag FUNCTION pythag_dp(a,b) USE nrtype REAL(DP), INTENT(IN) :: a,b REAL(DP) :: pythag_dp END FUNCTION pythag_dp !BL FUNCTION pythag_sp(a,b) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP) :: pythag_sp END FUNCTION pythag_sp END INTERFACE INTERFACE SUBROUTINE pzextr(iest,xest,yest,yz,dy) USE nrtype INTEGER(I4B), INTENT(IN) :: iest REAL(SP), INTENT(IN) :: xest REAL(SP), DIMENSION(:), INTENT(IN) :: yest REAL(SP), DIMENSION(:), INTENT(OUT) :: yz,dy END SUBROUTINE pzextr END INTERFACE INTERFACE SUBROUTINE qrdcmp(a,c,d,sing) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:), INTENT(OUT) :: c,d LOGICAL(LGT), INTENT(OUT) :: sing END SUBROUTINE qrdcmp END INTERFACE INTERFACE FUNCTION qromb(func,a,b) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP) :: qromb INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END FUNCTION qromb END INTERFACE INTERFACE FUNCTION qromo(func,a,b,choose) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP) :: qromo INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE INTERFACE SUBROUTINE choose(funk,aa,bb,s,n) USE nrtype REAL(SP), INTENT(IN) :: aa,bb REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION funk(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: funk END FUNCTION funk END INTERFACE END SUBROUTINE choose END INTERFACE END FUNCTION qromo END INTERFACE INTERFACE SUBROUTINE qroot(p,b,c,eps) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: p REAL(SP), INTENT(INOUT) :: b,c REAL(SP), INTENT(IN) :: eps END SUBROUTINE qroot END INTERFACE INTERFACE SUBROUTINE qrsolv(a,c,d,b) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: a REAL(SP), DIMENSION(:), INTENT(IN) :: c,d REAL(SP), DIMENSION(:), INTENT(INOUT) :: b END SUBROUTINE qrsolv END INTERFACE INTERFACE SUBROUTINE qrupdt(r,qt,u,v) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: r,qt REAL(SP), DIMENSION(:), INTENT(INOUT) :: u REAL(SP), DIMENSION(:), INTENT(IN) :: v END SUBROUTINE qrupdt END INTERFACE INTERFACE FUNCTION qsimp(func,a,b) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP) :: qsimp INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END FUNCTION qsimp END INTERFACE INTERFACE FUNCTION qtrap(func,a,b) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP) :: qtrap INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END FUNCTION qtrap END INTERFACE INTERFACE SUBROUTINE quadct(x,y,xx,yy,fa,fb,fc,fd) USE nrtype REAL(SP), INTENT(IN) :: x,y REAL(SP), DIMENSION(:), INTENT(IN) :: xx,yy REAL(SP), INTENT(OUT) :: fa,fb,fc,fd END SUBROUTINE quadct END INTERFACE INTERFACE SUBROUTINE quadmx(a) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(OUT) :: a END SUBROUTINE quadmx END INTERFACE INTERFACE SUBROUTINE quadvl(x,y,fa,fb,fc,fd) USE nrtype REAL(SP), INTENT(IN) :: x,y REAL(SP), INTENT(OUT) :: fa,fb,fc,fd END SUBROUTINE quadvl END INTERFACE INTERFACE FUNCTION ran(idum) INTEGER(selected_int_kind(9)), INTENT(INOUT) :: idum REAL :: ran END FUNCTION ran END INTERFACE INTERFACE ran0 SUBROUTINE ran0_s(harvest) USE nrtype REAL(SP), INTENT(OUT) :: harvest END SUBROUTINE ran0_s !BL SUBROUTINE ran0_v(harvest) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest END SUBROUTINE ran0_v END INTERFACE INTERFACE ran1 SUBROUTINE ran1_s(harvest) USE nrtype REAL(SP), INTENT(OUT) :: harvest END SUBROUTINE ran1_s !BL SUBROUTINE ran1_v(harvest) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest END SUBROUTINE ran1_v END INTERFACE INTERFACE ran2 SUBROUTINE ran2_s(harvest) USE nrtype REAL(SP), INTENT(OUT) :: harvest END SUBROUTINE ran2_s !BL SUBROUTINE ran2_v(harvest) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest END SUBROUTINE ran2_v END INTERFACE INTERFACE ran3 SUBROUTINE ran3_s(harvest) USE nrtype REAL(SP), INTENT(OUT) :: harvest END SUBROUTINE ran3_s !BL SUBROUTINE ran3_v(harvest) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: harvest END SUBROUTINE ran3_v END INTERFACE INTERFACE SUBROUTINE ratint(xa,ya,x,y,dy) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: y,dy END SUBROUTINE ratint END INTERFACE INTERFACE SUBROUTINE ratlsq(func,a,b,mm,kk,cof,dev) USE nrtype REAL(DP), INTENT(IN) :: a,b INTEGER(I4B), INTENT(IN) :: mm,kk REAL(DP), DIMENSION(:), INTENT(OUT) :: cof REAL(DP), INTENT(OUT) :: dev INTERFACE FUNCTION func(x) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: x REAL(DP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END SUBROUTINE ratlsq END INTERFACE INTERFACE ratval FUNCTION ratval_s(x,cof,mm,kk) USE nrtype REAL(DP), INTENT(IN) :: x INTEGER(I4B), INTENT(IN) :: mm,kk REAL(DP), DIMENSION(mm+kk+1), INTENT(IN) :: cof REAL(DP) :: ratval_s END FUNCTION ratval_s !BL FUNCTION ratval_v(x,cof,mm,kk) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: x INTEGER(I4B), INTENT(IN) :: mm,kk REAL(DP), DIMENSION(mm+kk+1), INTENT(IN) :: cof REAL(DP), DIMENSION(size(x)) :: ratval_v END FUNCTION ratval_v END INTERFACE INTERFACE rc FUNCTION rc_s(x,y) USE nrtype REAL(SP), INTENT(IN) :: x,y REAL(SP) :: rc_s END FUNCTION rc_s !BL FUNCTION rc_v(x,y) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), DIMENSION(size(x)) :: rc_v END FUNCTION rc_v END INTERFACE INTERFACE rd FUNCTION rd_s(x,y,z) USE nrtype REAL(SP), INTENT(IN) :: x,y,z REAL(SP) :: rd_s END FUNCTION rd_s !BL FUNCTION rd_v(x,y,z) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,z REAL(SP), DIMENSION(size(x)) :: rd_v END FUNCTION rd_v END INTERFACE INTERFACE realft SUBROUTINE realft_dp(data,isign,zdata) USE nrtype REAL(DP), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign COMPLEX(DPC), DIMENSION(:), OPTIONAL, TARGET :: zdata END SUBROUTINE realft_dp !BL SUBROUTINE realft_sp(data,isign,zdata) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: data INTEGER(I4B), INTENT(IN) :: isign COMPLEX(SPC), DIMENSION(:), OPTIONAL, TARGET :: zdata END SUBROUTINE realft_sp END INTERFACE INTERFACE RECURSIVE FUNCTION recur1(a,b) RESULT(u) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,b REAL(SP), DIMENSION(size(a)) :: u END FUNCTION recur1 END INTERFACE INTERFACE FUNCTION recur2(a,b,c) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,c REAL(SP), DIMENSION(size(a)) :: recur2 END FUNCTION recur2 END INTERFACE INTERFACE SUBROUTINE relax(u,rhs) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u REAL(DP), DIMENSION(:,:), INTENT(IN) :: rhs END SUBROUTINE relax END INTERFACE INTERFACE SUBROUTINE relax2(u,rhs) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u REAL(DP), DIMENSION(:,:), INTENT(IN) :: rhs END SUBROUTINE relax2 END INTERFACE INTERFACE FUNCTION resid(u,rhs) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: u,rhs REAL(DP), DIMENSION(size(u,1),size(u,1)) :: resid END FUNCTION resid END INTERFACE INTERFACE rf FUNCTION rf_s(x,y,z) USE nrtype REAL(SP), INTENT(IN) :: x,y,z REAL(SP) :: rf_s END FUNCTION rf_s !BL FUNCTION rf_v(x,y,z) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,z REAL(SP), DIMENSION(size(x)) :: rf_v END FUNCTION rf_v END INTERFACE INTERFACE rj FUNCTION rj_s(x,y,z,p) USE nrtype REAL(SP), INTENT(IN) :: x,y,z,p REAL(SP) :: rj_s END FUNCTION rj_s !BL FUNCTION rj_v(x,y,z,p) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,z,p REAL(SP), DIMENSION(size(x)) :: rj_v END FUNCTION rj_v END INTERFACE INTERFACE SUBROUTINE rk4(y,dydx,x,h,yout,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx REAL(SP), INTENT(IN) :: x,h REAL(SP), DIMENSION(:), INTENT(OUT) :: yout INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE rk4 END INTERFACE INTERFACE SUBROUTINE rkck(y,dydx,x,h,yout,yerr,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx REAL(SP), INTENT(IN) :: x,h REAL(SP), DIMENSION(:), INTENT(OUT) :: yout,yerr INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE rkck END INTERFACE INTERFACE SUBROUTINE rkdumb(vstart,x1,x2,nstep,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: vstart REAL(SP), INTENT(IN) :: x1,x2 INTEGER(I4B), INTENT(IN) :: nstep INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE rkdumb END INTERFACE INTERFACE SUBROUTINE rkqs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal REAL(SP), INTENT(INOUT) :: x REAL(SP), INTENT(IN) :: htry,eps REAL(SP), INTENT(OUT) :: hdid,hnext INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE rkqs END INTERFACE INTERFACE SUBROUTINE rlft2(data,spec,speq,isign) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: data COMPLEX(SPC), DIMENSION(:,:), INTENT(OUT) :: spec COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: speq INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE rlft2 END INTERFACE INTERFACE SUBROUTINE rlft3(data,spec,speq,isign) USE nrtype REAL(SP), DIMENSION(:,:,:), INTENT(INOUT) :: data COMPLEX(SPC), DIMENSION(:,:,:), INTENT(OUT) :: spec COMPLEX(SPC), DIMENSION(:,:), INTENT(OUT) :: speq INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE rlft3 END INTERFACE INTERFACE SUBROUTINE rotate(r,qt,i,a,b) USE nrtype REAL(SP), DIMENSION(:,:), TARGET, INTENT(INOUT) :: r,qt INTEGER(I4B), INTENT(IN) :: i REAL(SP), INTENT(IN) :: a,b END SUBROUTINE rotate END INTERFACE INTERFACE SUBROUTINE rsolv(a,d,b) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: a REAL(SP), DIMENSION(:), INTENT(IN) :: d REAL(SP), DIMENSION(:), INTENT(INOUT) :: b END SUBROUTINE rsolv END INTERFACE INTERFACE FUNCTION rstrct(uf) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: uf REAL(DP), DIMENSION((size(uf,1)+1)/2,(size(uf,1)+1)/2) :: rstrct END FUNCTION rstrct END INTERFACE INTERFACE FUNCTION rtbis(func,x1,x2,xacc) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,xacc REAL(SP) :: rtbis INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION rtbis END INTERFACE INTERFACE FUNCTION rtflsp(func,x1,x2,xacc) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,xacc REAL(SP) :: rtflsp INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION rtflsp END INTERFACE INTERFACE FUNCTION rtnewt(funcd,x1,x2,xacc) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,xacc REAL(SP) :: rtnewt INTERFACE SUBROUTINE funcd(x,fval,fderiv) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: fval,fderiv END SUBROUTINE funcd END INTERFACE END FUNCTION rtnewt END INTERFACE INTERFACE FUNCTION rtsafe(funcd,x1,x2,xacc) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,xacc REAL(SP) :: rtsafe INTERFACE SUBROUTINE funcd(x,fval,fderiv) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: fval,fderiv END SUBROUTINE funcd END INTERFACE END FUNCTION rtsafe END INTERFACE INTERFACE FUNCTION rtsec(func,x1,x2,xacc) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,xacc REAL(SP) :: rtsec INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION rtsec END INTERFACE INTERFACE SUBROUTINE rzextr(iest,xest,yest,yz,dy) USE nrtype INTEGER(I4B), INTENT(IN) :: iest REAL(SP), INTENT(IN) :: xest REAL(SP), DIMENSION(:), INTENT(IN) :: yest REAL(SP), DIMENSION(:), INTENT(OUT) :: yz,dy END SUBROUTINE rzextr END INTERFACE INTERFACE FUNCTION savgol(nl,nrr,ld,m) USE nrtype INTEGER(I4B), INTENT(IN) :: nl,nrr,ld,m REAL(SP), DIMENSION(nl+nrr+1) :: savgol END FUNCTION savgol END INTERFACE INTERFACE SUBROUTINE scrsho(func) USE nrtype INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE scrsho END INTERFACE INTERFACE FUNCTION select(k,arr) USE nrtype INTEGER(I4B), INTENT(IN) :: k REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr REAL(SP) :: select END FUNCTION select END INTERFACE INTERFACE FUNCTION select_bypack(k,arr) USE nrtype INTEGER(I4B), INTENT(IN) :: k REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr REAL(SP) :: select_bypack END FUNCTION select_bypack END INTERFACE INTERFACE SUBROUTINE select_heap(arr,heap) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: arr REAL(SP), DIMENSION(:), INTENT(OUT) :: heap END SUBROUTINE select_heap END INTERFACE INTERFACE FUNCTION select_inplace(k,arr) USE nrtype INTEGER(I4B), INTENT(IN) :: k REAL(SP), DIMENSION(:), INTENT(IN) :: arr REAL(SP) :: select_inplace END FUNCTION select_inplace END INTERFACE INTERFACE SUBROUTINE simplx(a,m1,m2,m3,icase,izrov,iposv) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a INTEGER(I4B), INTENT(IN) :: m1,m2,m3 INTEGER(I4B), INTENT(OUT) :: icase INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: izrov,iposv END SUBROUTINE simplx END INTERFACE INTERFACE SUBROUTINE simpr(y,dydx,dfdx,dfdy,xs,htot,nstep,yout,derivs) USE nrtype REAL(SP), INTENT(IN) :: xs,htot REAL(SP), DIMENSION(:), INTENT(IN) :: y,dydx,dfdx REAL(SP), DIMENSION(:,:), INTENT(IN) :: dfdy INTEGER(I4B), INTENT(IN) :: nstep REAL(SP), DIMENSION(:), INTENT(OUT) :: yout INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE simpr END INTERFACE INTERFACE SUBROUTINE sinft(y) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y END SUBROUTINE sinft END INTERFACE INTERFACE SUBROUTINE slvsm2(u,rhs) USE nrtype REAL(DP), DIMENSION(3,3), INTENT(OUT) :: u REAL(DP), DIMENSION(3,3), INTENT(IN) :: rhs END SUBROUTINE slvsm2 END INTERFACE INTERFACE SUBROUTINE slvsml(u,rhs) USE nrtype REAL(DP), DIMENSION(3,3), INTENT(OUT) :: u REAL(DP), DIMENSION(3,3), INTENT(IN) :: rhs END SUBROUTINE slvsml END INTERFACE INTERFACE SUBROUTINE sncndn(uu,emmc,sn,cn,dn) USE nrtype REAL(SP), INTENT(IN) :: uu,emmc REAL(SP), INTENT(OUT) :: sn,cn,dn END SUBROUTINE sncndn END INTERFACE INTERFACE FUNCTION snrm(sx,itol) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: sx INTEGER(I4B), INTENT(IN) :: itol REAL(DP) :: snrm END FUNCTION snrm END INTERFACE INTERFACE SUBROUTINE sobseq(x,init) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: x INTEGER(I4B), OPTIONAL, INTENT(IN) :: init END SUBROUTINE sobseq END INTERFACE INTERFACE SUBROUTINE solvde(itmax,conv,slowc,scalv,indexv,nb,y) USE nrtype INTEGER(I4B), INTENT(IN) :: itmax,nb REAL(SP), INTENT(IN) :: conv,slowc REAL(SP), DIMENSION(:), INTENT(IN) :: scalv INTEGER(I4B), DIMENSION(:), INTENT(IN) :: indexv REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: y END SUBROUTINE solvde END INTERFACE INTERFACE SUBROUTINE sor(a,b,c,d,e,f,u,rjac) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: a,b,c,d,e,f REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: u REAL(DP), INTENT(IN) :: rjac END SUBROUTINE sor END INTERFACE INTERFACE SUBROUTINE sort(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort END INTERFACE INTERFACE SUBROUTINE sort2(arr,slave) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr,slave END SUBROUTINE sort2 END INTERFACE INTERFACE SUBROUTINE sort3(arr,slave1,slave2) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr,slave1,slave2 END SUBROUTINE sort3 END INTERFACE INTERFACE SUBROUTINE sort_bypack(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort_bypack END INTERFACE INTERFACE SUBROUTINE sort_byreshape(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort_byreshape END INTERFACE INTERFACE SUBROUTINE sort_heap(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort_heap END INTERFACE INTERFACE SUBROUTINE sort_pick(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort_pick END INTERFACE INTERFACE SUBROUTINE sort_radix(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort_radix END INTERFACE INTERFACE SUBROUTINE sort_shell(arr) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: arr END SUBROUTINE sort_shell END INTERFACE INTERFACE SUBROUTINE spctrm(p,k,ovrlap,unit,n_window) USE nrtype REAL(SP), DIMENSION(:), INTENT(OUT) :: p INTEGER(I4B), INTENT(IN) :: k LOGICAL(LGT), INTENT(IN) :: ovrlap INTEGER(I4B), OPTIONAL, INTENT(IN) :: n_window,unit END SUBROUTINE spctrm END INTERFACE INTERFACE SUBROUTINE spear(data1,data2,d,zd,probd,rs,probrs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 REAL(SP), INTENT(OUT) :: d,zd,probd,rs,probrs END SUBROUTINE spear END INTERFACE INTERFACE sphbes SUBROUTINE sphbes_s(n,x,sj,sy,sjp,syp) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: x REAL(SP), INTENT(OUT) :: sj,sy,sjp,syp END SUBROUTINE sphbes_s !BL SUBROUTINE sphbes_v(n,x,sj,sy,sjp,syp) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(OUT) :: sj,sy,sjp,syp END SUBROUTINE sphbes_v END INTERFACE INTERFACE SUBROUTINE splie2(x1a,x2a,ya,y2a) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x1a,x2a REAL(SP), DIMENSION(:,:), INTENT(IN) :: ya REAL(SP), DIMENSION(:,:), INTENT(OUT) :: y2a END SUBROUTINE splie2 END INTERFACE INTERFACE FUNCTION splin2(x1a,x2a,ya,y2a,x1,x2) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x1a,x2a REAL(SP), DIMENSION(:,:), INTENT(IN) :: ya,y2a REAL(SP), INTENT(IN) :: x1,x2 REAL(SP) :: splin2 END FUNCTION splin2 END INTERFACE INTERFACE SUBROUTINE spline(x,y,yp1,ypn,y2) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y REAL(SP), INTENT(IN) :: yp1,ypn REAL(SP), DIMENSION(:), INTENT(OUT) :: y2 END SUBROUTINE spline END INTERFACE INTERFACE FUNCTION splint(xa,ya,y2a,x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: xa,ya,y2a REAL(SP), INTENT(IN) :: x REAL(SP) :: splint END FUNCTION splint END INTERFACE INTERFACE sprsax SUBROUTINE sprsax_dp(sa,x,b) USE nrtype TYPE(sprs2_dp), INTENT(IN) :: sa REAL(DP), DIMENSION (:), INTENT(IN) :: x REAL(DP), DIMENSION (:), INTENT(OUT) :: b END SUBROUTINE sprsax_dp !BL SUBROUTINE sprsax_sp(sa,x,b) USE nrtype TYPE(sprs2_sp), INTENT(IN) :: sa REAL(SP), DIMENSION (:), INTENT(IN) :: x REAL(SP), DIMENSION (:), INTENT(OUT) :: b END SUBROUTINE sprsax_sp END INTERFACE INTERFACE sprsdiag SUBROUTINE sprsdiag_dp(sa,b) USE nrtype TYPE(sprs2_dp), INTENT(IN) :: sa REAL(DP), DIMENSION(:), INTENT(OUT) :: b END SUBROUTINE sprsdiag_dp !BL SUBROUTINE sprsdiag_sp(sa,b) USE nrtype TYPE(sprs2_sp), INTENT(IN) :: sa REAL(SP), DIMENSION(:), INTENT(OUT) :: b END SUBROUTINE sprsdiag_sp END INTERFACE INTERFACE sprsin SUBROUTINE sprsin_sp(a,thresh,sa) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: a REAL(SP), INTENT(IN) :: thresh TYPE(sprs2_sp), INTENT(OUT) :: sa END SUBROUTINE sprsin_sp !BL SUBROUTINE sprsin_dp(a,thresh,sa) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: a REAL(DP), INTENT(IN) :: thresh TYPE(sprs2_dp), INTENT(OUT) :: sa END SUBROUTINE sprsin_dp END INTERFACE INTERFACE SUBROUTINE sprstp(sa) USE nrtype TYPE(sprs2_sp), INTENT(INOUT) :: sa END SUBROUTINE sprstp END INTERFACE INTERFACE sprstx SUBROUTINE sprstx_dp(sa,x,b) USE nrtype TYPE(sprs2_dp), INTENT(IN) :: sa REAL(DP), DIMENSION (:), INTENT(IN) :: x REAL(DP), DIMENSION (:), INTENT(OUT) :: b END SUBROUTINE sprstx_dp !BL SUBROUTINE sprstx_sp(sa,x,b) USE nrtype TYPE(sprs2_sp), INTENT(IN) :: sa REAL(SP), DIMENSION (:), INTENT(IN) :: x REAL(SP), DIMENSION (:), INTENT(OUT) :: b END SUBROUTINE sprstx_sp END INTERFACE INTERFACE SUBROUTINE stifbs(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal REAL(SP), INTENT(IN) :: htry,eps REAL(SP), INTENT(INOUT) :: x REAL(SP), INTENT(OUT) :: hdid,hnext INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE stifbs END INTERFACE INTERFACE SUBROUTINE stiff(y,dydx,x,htry,eps,yscal,hdid,hnext,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: y REAL(SP), DIMENSION(:), INTENT(IN) :: dydx,yscal REAL(SP), INTENT(INOUT) :: x REAL(SP), INTENT(IN) :: htry,eps REAL(SP), INTENT(OUT) :: hdid,hnext INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE stiff END INTERFACE INTERFACE SUBROUTINE stoerm(y,d2y,xs,htot,nstep,yout,derivs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: y,d2y REAL(SP), INTENT(IN) :: xs,htot INTEGER(I4B), INTENT(IN) :: nstep REAL(SP), DIMENSION(:), INTENT(OUT) :: yout INTERFACE SUBROUTINE derivs(x,y,dydx) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP), DIMENSION(:), INTENT(IN) :: y REAL(SP), DIMENSION(:), INTENT(OUT) :: dydx END SUBROUTINE derivs END INTERFACE END SUBROUTINE stoerm END INTERFACE INTERFACE svbksb SUBROUTINE svbksb_dp(u,w,v,b,x) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(IN) :: u,v REAL(DP), DIMENSION(:), INTENT(IN) :: w,b REAL(DP), DIMENSION(:), INTENT(OUT) :: x END SUBROUTINE svbksb_dp !BL SUBROUTINE svbksb_sp(u,w,v,b,x) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: u,v REAL(SP), DIMENSION(:), INTENT(IN) :: w,b REAL(SP), DIMENSION(:), INTENT(OUT) :: x END SUBROUTINE svbksb_sp END INTERFACE INTERFACE svdcmp SUBROUTINE svdcmp_dp(a,w,v) USE nrtype REAL(DP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(DP), DIMENSION(:), INTENT(OUT) :: w REAL(DP), DIMENSION(:,:), INTENT(OUT) :: v END SUBROUTINE svdcmp_dp !BL SUBROUTINE svdcmp_sp(a,w,v) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:), INTENT(OUT) :: w REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v END SUBROUTINE svdcmp_sp END INTERFACE INTERFACE SUBROUTINE svdfit(x,y,sig,a,v,w,chisq,funcs) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x,y,sig REAL(SP), DIMENSION(:), INTENT(OUT) :: a,w REAL(SP), DIMENSION(:,:), INTENT(OUT) :: v REAL(SP), INTENT(OUT) :: chisq INTERFACE FUNCTION funcs(x,n) USE nrtype REAL(SP), INTENT(IN) :: x INTEGER(I4B), INTENT(IN) :: n REAL(SP), DIMENSION(n) :: funcs END FUNCTION funcs END INTERFACE END SUBROUTINE svdfit END INTERFACE INTERFACE SUBROUTINE svdvar(v,w,cvm) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(IN) :: v REAL(SP), DIMENSION(:), INTENT(IN) :: w REAL(SP), DIMENSION(:,:), INTENT(OUT) :: cvm END SUBROUTINE svdvar END INTERFACE INTERFACE FUNCTION toeplz(r,y) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: r,y REAL(SP), DIMENSION(size(y)) :: toeplz END FUNCTION toeplz END INTERFACE INTERFACE SUBROUTINE tptest(data1,data2,t,prob) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 REAL(SP), INTENT(OUT) :: t,prob END SUBROUTINE tptest END INTERFACE INTERFACE SUBROUTINE tqli(d,e,z) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: d,e REAL(SP), DIMENSION(:,:), OPTIONAL, INTENT(INOUT) :: z END SUBROUTINE tqli END INTERFACE INTERFACE SUBROUTINE trapzd(func,a,b,s,n) USE nrtype REAL(SP), INTENT(IN) :: a,b REAL(SP), INTENT(INOUT) :: s INTEGER(I4B), INTENT(IN) :: n INTERFACE FUNCTION func(x) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: x REAL(SP), DIMENSION(size(x)) :: func END FUNCTION func END INTERFACE END SUBROUTINE trapzd END INTERFACE INTERFACE SUBROUTINE tred2(a,d,e,novectors) USE nrtype REAL(SP), DIMENSION(:,:), INTENT(INOUT) :: a REAL(SP), DIMENSION(:), INTENT(OUT) :: d,e LOGICAL(LGT), OPTIONAL, INTENT(IN) :: novectors END SUBROUTINE tred2 END INTERFACE ! On a purely serial machine, for greater efficiency, remove ! the generic name tridag from the following interface, ! and put it on the next one after that. INTERFACE tridag RECURSIVE SUBROUTINE tridag_par(a,b,c,r,u) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,c,r REAL(SP), DIMENSION(:), INTENT(OUT) :: u END SUBROUTINE tridag_par END INTERFACE INTERFACE SUBROUTINE tridag_ser(a,b,c,r,u) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a,b,c,r REAL(SP), DIMENSION(:), INTENT(OUT) :: u END SUBROUTINE tridag_ser END INTERFACE INTERFACE SUBROUTINE ttest(data1,data2,t,prob) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 REAL(SP), INTENT(OUT) :: t,prob END SUBROUTINE ttest END INTERFACE INTERFACE SUBROUTINE tutest(data1,data2,t,prob) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 REAL(SP), INTENT(OUT) :: t,prob END SUBROUTINE tutest END INTERFACE INTERFACE SUBROUTINE twofft(data1,data2,fft1,fft2) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: data1,data2 COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: fft1,fft2 END SUBROUTINE twofft END INTERFACE INTERFACE FUNCTION vander(x,q) USE nrtype REAL(DP), DIMENSION(:), INTENT(IN) :: x,q REAL(DP), DIMENSION(size(x)) :: vander END FUNCTION vander END INTERFACE INTERFACE SUBROUTINE vegas(region,func,init,ncall,itmx,nprn,tgral,sd,chi2a) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: region INTEGER(I4B), INTENT(IN) :: init,ncall,itmx,nprn REAL(SP), INTENT(OUT) :: tgral,sd,chi2a INTERFACE FUNCTION func(pt,wgt) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: pt REAL(SP), INTENT(IN) :: wgt REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE vegas END INTERFACE INTERFACE SUBROUTINE voltra(t0,h,t,f,g,ak) USE nrtype REAL(SP), INTENT(IN) :: t0,h REAL(SP), DIMENSION(:), INTENT(OUT) :: t REAL(SP), DIMENSION(:,:), INTENT(OUT) :: f INTERFACE FUNCTION g(t) USE nrtype REAL(SP), INTENT(IN) :: t REAL(SP), DIMENSION(:), POINTER :: g END FUNCTION g !BL FUNCTION ak(t,s) USE nrtype REAL(SP), INTENT(IN) :: t,s REAL(SP), DIMENSION(:,:), POINTER :: ak END FUNCTION ak END INTERFACE END SUBROUTINE voltra END INTERFACE INTERFACE SUBROUTINE wt1(a,isign,wtstep) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: a INTEGER(I4B), INTENT(IN) :: isign INTERFACE SUBROUTINE wtstep(a,isign) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: a INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE wtstep END INTERFACE END SUBROUTINE wt1 END INTERFACE INTERFACE SUBROUTINE wtn(a,nn,isign,wtstep) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: a INTEGER(I4B), DIMENSION(:), INTENT(IN) :: nn INTEGER(I4B), INTENT(IN) :: isign INTERFACE SUBROUTINE wtstep(a,isign) USE nrtype REAL(SP), DIMENSION(:), INTENT(INOUT) :: a INTEGER(I4B), INTENT(IN) :: isign END SUBROUTINE wtstep END INTERFACE END SUBROUTINE wtn END INTERFACE INTERFACE FUNCTION wwghts(n,h,kermom) USE nrtype INTEGER(I4B), INTENT(IN) :: n REAL(SP), INTENT(IN) :: h REAL(SP), DIMENSION(n) :: wwghts INTERFACE FUNCTION kermom(y,m) USE nrtype REAL(DP), INTENT(IN) :: y INTEGER(I4B), INTENT(IN) :: m REAL(DP), DIMENSION(m) :: kermom END FUNCTION kermom END INTERFACE END FUNCTION wwghts END INTERFACE INTERFACE SUBROUTINE zbrac(func,x1,x2,succes) USE nrtype REAL(SP), INTENT(INOUT) :: x1,x2 LOGICAL(LGT), INTENT(OUT) :: succes INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE zbrac END INTERFACE INTERFACE SUBROUTINE zbrak(func,x1,x2,n,xb1,xb2,nb) USE nrtype INTEGER(I4B), INTENT(IN) :: n INTEGER(I4B), INTENT(OUT) :: nb REAL(SP), INTENT(IN) :: x1,x2 REAL(SP), DIMENSION(:), POINTER :: xb1,xb2 INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END SUBROUTINE zbrak END INTERFACE INTERFACE FUNCTION zbrent(func,x1,x2,tol) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,tol REAL(SP) :: zbrent INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION zbrent END INTERFACE INTERFACE SUBROUTINE zrhqr(a,rtr,rti) USE nrtype REAL(SP), DIMENSION(:), INTENT(IN) :: a REAL(SP), DIMENSION(:), INTENT(OUT) :: rtr,rti END SUBROUTINE zrhqr END INTERFACE INTERFACE FUNCTION zriddr(func,x1,x2,xacc) USE nrtype REAL(SP), INTENT(IN) :: x1,x2,xacc REAL(SP) :: zriddr INTERFACE FUNCTION func(x) USE nrtype REAL(SP), INTENT(IN) :: x REAL(SP) :: func END FUNCTION func END INTERFACE END FUNCTION zriddr END INTERFACE INTERFACE SUBROUTINE zroots(a,roots,polish) USE nrtype COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: roots LOGICAL(LGT), INTENT(IN) :: polish END SUBROUTINE zroots END INTERFACE END MODULE nr module stat_elementaires use constants_module; use nrtype use mod_mat implicit none private private ::autocovmult public :: tri, moy, sig, med, quart, autocov, autocor, autocov0, autocor0 !tri(x) est le fichier x trie par ordre croissant contains function autocovmult(x,k) result(autocov_result) real (SP), intent(in), dimension(:,:) :: x integer(I4B), intent(in) :: k real (SP) :: autocov_result(0:k,1:size(x,2),1:size(x,2)),m(1:size(x,2)) integer(I4B) :: n, d, i,t n=size(x,1) d=size(x,2) m(1:d)=sum(x,1)/n autocov_result=zero do i=0,k do t=1,n-i autocov_result(i,1:d,1:d)=autocov_result(i,1:d,1:d)+prod(x(t,1:d),x(t+i,1:d)) end do autocov_result(i,1:d,1:d)=(autocov_result(i,1:d,1:d)/n)-prod(m,m) end do return end function autocovmult function autocov(x,k) result(autocov_result) real (SP), intent(in), dimension(:) :: x integer(I4B), intent(in) :: k real (SP) :: autocov_result integer(I4B) :: n, kk n=size(x) kk=k if (k < 0)then kk=-k end if autocov_result=dot_product(x(1:n-kk)-moy(x(1:n)),x(kk+1:n)-moy(x(1:n)))/n return end function autocov function autocor(x,k) result(autocor_result) real (SP), intent(in), dimension(:) :: x integer(I4B), intent(in) :: k real (SP) :: autocor_result integer(I4B) :: kk kk=k if (k < 0) then kk=-k end if autocor_result=autocov(x,kk)/autocov(x,0) return end function autocor function autocov0(x,k) result(autocov_result) real (SP), intent(in), dimension(:) :: x integer(I4B), intent(in) :: k real (SP) :: autocov_result integer(I4B) :: n, kk n=size(x) kk=k if (k < 0)then kk=-k end if autocov_result=dot_product(x(1:n-kk),x(kk+1:n))/n return end function autocov0 function autocor0(x,k) result(autocor_result) real (SP), intent(in), dimension(:) :: x integer(I4B), intent(in) :: k real (SP) :: autocor_result integer(I4B) :: kk kk=k if (k < 0) then kk=-k end if autocor_result=autocov0(x,kk)/autocov0(x,0) return end function autocor0 function tri(x) result(tri_result) !tri(x) est le fichier x trie par ordre croissant real (SP), intent(in),dimension(:) :: x real (SP),dimension(size(x)) :: tri_result integer(I4B) :: n, l, ir, i, j real (SP) :: rra n=size(x) tri_result=x if(n<=1)return l=n/2+1 ir=n boucle1 : do if (l > 1) then l=l-1 rra=tri_result(l) else rra=tri_result(ir) tri_result(ir)=tri_result(1) ir=ir-1 if (ir == 1) then tri_result(1)=rra return end if end if i=l j=l+l boucle2 : do if (j <= ir) then if (j < ir ) then if (tri_result(j) < tri_result(j+1))then j=j+1 end if end if if( rra < tri_result(j)) then tri_result(i)=tri_result(j) i=j j=j+j else j=ir+1 end if else exit boucle2 end if end do boucle2 tri_result(i)=rra end do boucle1 return end function tri function med(x) result(med_result) real (SP), intent(in),dimension(:) :: x real (SP) :: med_result integer(I4B) :: n,id real (SP) :: rest real (SP),dimension(size(x)) :: xx n=size(x) xx=tri(x) id=int(1.0_sp*n/2) rest=1.0_sp*n/2 -id if (rest == 0.00_sp) then med_result=(xx(id)+xx(id+1))/2 else med_result=xx(id+1) end if return end function med subroutine quart(x,mmin,q1,q2,q3,mmax) real (SP), intent(in),dimension(:) :: x real (SP), intent(out) :: mmin, q1, q2, q3, mmax integer(I4B) :: n real (SP) :: rest real (SP), dimension(size(x)) :: xx n=size(x) xx=tri(x) mmin=xx(1) mmax=xx(n) q2=med(xx) rest=1.0_sp*n/2 -int(1.0_sp*n/2) if (rest == 0.00_sp) then q1=med(xx(1:int(1.0_sp*n/2))) q3=med(xx(int(1.0_sp*n/2)+1:n)) else q1=med(xx(1:int(1.0_sp*n/2))) q3=med(xx(int(1.0_sp*n/2)+2:n)) end if return end subroutine quart function moy(x) result(moy_result) real (SP), intent(in), dimension (:) :: x real (SP) :: moy_result integer(I4B) :: n n=size(x) moy_result=sum(x)/n return end function moy function sig(x) result(sig_result) real (SP), intent(in),dimension(:) :: x real (SP) :: sig_result integer(I4B) :: n n=size(x) sig_result=(dot_product(x,x)/n)-(sum(x)/n)**2 sig_result=sqrt(sig_result) return end function sig end module stat_elementaires module fdr ! quelque fonctions de repartitions ! Phi(x):=fonction de repartition d'une normale standard au point x ! Phiinv(p):=Phi^{-1}(p) ! Khi2(x,nu):=fonction de repartition au point x d'une khi2 a nu degres de libertes ! Khi2inv(p,nu):=Khi2inv^{-1}(p,nu) use constants_module; use nrtype implicit none private public :: norm,Khi2b, gammp,gammln,fKhi2, gamma, Khi2inv, Khi2, Phi, Phiinv, invfdr, invfdr2 contains function norm(x) ! norm(x)=densite d'une normale standard au point x implicit none real (SP), intent(in) :: x real (SP) :: norm norm=(un/(sqrt(deux*pi)))*exp(-(x**2)/deux) return end function norm function Khi2(x,nu) ! Khi2(nu,x):=fonction de repartition au point x d'une khi2 a nu degres de libertes implicit none integer(I4B), intent(in) :: nu real (SP), intent(in) :: x real (SP) :: Khi2 if(x<=zero)then Khi2=zero else Khi2=gammp(un*nu/2,x/2) end if return end function Khi2 function gammp(a,x) real(SP),intent(in) :: a,x real(SP) :: gammp,gammcf,gln if(x.lt.zero.or.a.le.zero)then write(8,*)"pb dans gammp" end if if(x.lt.a+un)then call gser(gammp,a,x,gln) else call gcf(gammcf,a,x,gln) gammp=un-gammcf end if return contains subroutine gser(gamser,a,x,gln) real(SP),intent(in) :: a,x real(SP),intent(out) :: gamser,gln real(SP), parameter :: eps=0.0000003 integer(I4B), parameter :: itmax=100 real(SP) :: ap,sum,del integer(I4B) :: n gln=gammln(a) if(x.le.zero)then if(x.lt.zero)write(8,*)"pb dans gser" gamser=zero return end if ap=a sum=un/a del=sum do n=1,itmax ap=ap+un del=del*x/ap sum=sum+del if(abs(del).lt.abs(sum)*eps)go to 1 end do write(8,*)"a too large, itmax too small" 1 gamser=sum*exp(-x+a*log(x)-gln) return end subroutine gser subroutine gcf(gammcf,a,x,gln) real(SP),intent(in) :: a,x real(SP),intent(out) :: gammcf,gln real(SP), parameter :: eps=0.0000003 integer(I4B), parameter :: itmax=100 real(SP) :: gold,a0,a1,b0,b1,fac,an,ana,anf,g integer(I4B) :: n gln=gammln(a) gold=zero a0=un a1=x b0=zero b1=un fac=un do n=1,itmax an=float(n) ana=an-a a0=(a1+a0*ana)*fac b0=(b1+b0*ana)*fac anf=an*fac a1=x*a0+anf*a1 b1=x*b0+anf*b1 if(a1.ne.zero)then fac=un/a1 g=b1*fac if(abs((g-gold)/g).lt.eps)go to 1 gold=g end if end do write(8,*)"a too large, itmax too small" 1 gammcf=exp(-x+a*alog(x)-gln)*g return end subroutine gcf end function gammp function gammln(xx) ! fonction gamma d'Euler real(SP) :: cof(6),stp,half,one,fpf,x,tmp,ser real (SP), intent(in) :: xx real(SP) :: gammln integer(I4B) :: j cof=(/76.18009173_sp,-86.50532033_sp,24.01409822_sp,& -1.231739516_sp,.00120858003_sp,-.00000536382_sp/) stp=2.50662827465_sp half=0.5_sp one=1.0_sp fpf=5.5_sp x=xx-one tmp=x+fpf tmp=(x+half)*log(tmp)-tmp ser=one do j=1,6 x=x+one ser=ser+cof(j)/x end do gammln=tmp+log(stp*ser) return end function gammln function fKhi2(x,m) ! densite du khi2 a m degres de liberte integer(I4B), intent(in) :: m real (SP), intent(in) :: x real(SP) :: fKhi2 fKhi2=zero if(x>zero)fKhi2=exp(-x/2)*(x**(demi*m-un))/(gamma(demi*m)*(deux**(demi*m))) return end function fKhi2 function gamma(xx) ! fonction gamma d'Euler real(SP) :: cof(6),stp,half,one,fpf,x,tmp,ser real (SP), intent(in) :: xx real(SP) :: gamma integer(I4B) :: j cof=(/76.18009173_sp,-86.50532033_sp,24.01409822_sp,& -1.231739516_sp,.00120858003_sp,-.00000536382_sp/) stp=2.50662827465_sp half=0.5_sp one=1.0_sp fpf=5.5_sp x=xx-one tmp=x+fpf tmp=(x+half)*log(tmp)-tmp ser=one do j=1,6 x=x+one ser=ser+cof(j)/x end do gamma=exp(tmp+log(stp*ser)) return end function gamma function invfdr2(f,nu,p) ! invfdr=f^{-1}(p) ou f est une fdr qui possede le parametre nu implicit none integer(I4B), intent(in) :: nu real (SP), intent(in) :: p real (SP) :: invfdr2 real (SP) :: xmin, xmax, middle real (SP), parameter :: tol=.0001 interface function f(x,n) use constants_module; use nrtype implicit none integer(I4B), intent(in) :: n real (SP), intent(in) :: x real (SP) :: f end function f end interface xmin=-100*un xmax=100*un if (f(xmin,nu) > p) then invfdr2=-huge(1.) return end if if (f(xmax,nu) < p) then invfdr2=huge(1.) return end if do if (xmax-xmin < tol) exit middle=(xmax+xmin)/2 if (f(middle,nu) == p) then invfdr2=middle return end if if (f(middle,nu) > p) xmax=middle if (f(middle,nu) < p) xmin=middle end do invfdr2=(xmax+xmin)/2 return end function invfdr2 function Khi2b(x,nu) ! Khi2(nu,x):=fonction de repartition au point x d'une khi2 a nu degres de libertes implicit none integer(I4B), intent(in) :: nu real (SP), intent(in) :: x real (SP) :: Khi2b real (SP) :: dum if (nu < 1 ) then write(8,*)"pb dans Khi2" return end if if (nu == 1) Khi2b=2*Phi(sqrt(x))-un if (nu == 2) Khi2b=un-exp(-x/2) if (nu > 2) then dum=sqrt(9*un*nu/2)*((x/nu)**(un/3)+2./(9*un*nu)-un) Khi2b=Phi(dum) end if return end function Khi2b function Khi2inv(prob,nu) ! valeur ayant une proba prob d'etre plus petite qu'un Khi2 a nu degres de libertes ! Khi2inv(proba,nu):=Khi2^{-1}(prob,nu) implicit none integer(I4B), intent(in) :: nu real (SP), intent(in) :: prob real (SP) :: Khi2inv Khi2inv = invfdr2(Khi2,nu,prob) return end function Khi2inv function invfdr(f,p) ! invfdr=f^{-1}(p) ou f est une fdr implicit none real (SP), intent(in) :: p real (SP) :: invfdr real (SP) :: xmin, xmax, middle real (SP), parameter :: tol=.0001 interface function f(x) use constants_module; use nrtype implicit none real (SP), intent(in) :: x real (SP) :: f end function f end interface xmin=-10. xmax=10. if (f(xmin) > p) then invfdr=-huge(1.) return end if if (f(xmax) < p) then invfdr=huge(1.) return end if do if (xmax-xmin < tol) exit middle=(xmax+xmin)/2 if (f(middle) == p) then invfdr=middle return end if if (f(middle) > p) xmax=middle if (f(middle) < p) xmin=middle end do invfdr=(xmax+xmin)/2 return end function invfdr function Phi(x) ! fonction de repartition d'une normale standard au point x implicit none real (SP), intent(in) :: x real (SP) :: Phi real (SP), parameter :: b1=0.319381530, & b2=-0.356563782, & b3=1.781477937, & b4=-1.821255978, & b5=1.330274429 real (SP) :: t,dum t=1./(1.+0.2316419*abs(x)) dum=1./(sqrt(2.*pi)) dum=dum*exp(-.5*x**2) dum=dum*(b1*t+b2*t**2+b3*t**3+b4*t**4+b5*t**5) Phi=1.-dum if (x < 0.0) Phi=1.-Phi return end function Phi function Phiinv(prob) ! valeur ayant une proba prob d'etre plus petite qu'une N(0,1) ! gaussd(prob)=F^{-1}(prob) ou F est la fdr d'une N(0,1) real(SP) :: Phiinv real(SP), intent(in) :: prob Phiinv = invfdr(Phi,prob) return end function Phiinv end module fdr module random ! random_normal() est la realisation d'une N(0,1) ! random_exponential( ) est la realisation d'une exponentielle de parametre 1 ! random_t(m) est la realisation d'une Student de parametre m ! subroutine chaine_Markov(p,p0,d,n,delta) ! genere une chaine de Markov ! delta(-49) est tiree selon la loi initiale p0 ! delta(-48),....,delta(n) est realisation d'une chaine de Markov ! a espace d'etats {1,...,d} et de matrice de transition p ! p(i,j)est la proba d'etre en j sachant que l'on etait en i la fois d'avant ! ! ar(a,n) donne une simulation de taille n d'un AR(p) x(t)=a(1)*x(t-1)+...+a(p)*x(t-p)+eta(t) ou eta(t) iid N(0,1) ! 500 valeurs initiales sont supprimees use constants_module; use nrtype use stat_elementaires use fdr implicit none private public :: random_normal, random_exponential, random_t, chaine_Markov,ar,simarma, faible, & garch, estimcov,estimcoviid contains subroutine simarma(a,b,sigma,xmoy,epsil,x) ! x(t)=a(1)x(t-1)+...+a(p)x(t-p) ! +sigma*epsilon(t)-b(1)sigma*epsilon(t-1)-....-b(q)sigma*epsilon(t-q) ! +constante ! pour t=1,...,n ! valeur initiale nulles implicit none REAL(SP), intent(in) :: a(:), b(:), xmoy, sigma, epsil(:) REAL(SP), intent(out) :: x(:) ! Local variables integer(I4B) :: i,j,p,q,n n=size(x) p=size(a) q=size(b) do i=1,n x(i)=sigma*epsil(i) do j=1,p if(i-j>=1)x(i)=x(i)+a(j)*x(i-j) end do do j=1,q if(i-j>=1)x(i)=x(i)-b(j)*sigma*epsil(i-j) end do end do x(1:n)=x(1:n)+xmoy return end subroutine simarma subroutine faible(epsil,m,n) ! si m=0 ! epsil(t)=eta(t) ! avec eta(t) suite i.i.d. de loi normale standard implicit none integer(I4B), intent(in) :: n,m REAL(SP), intent(out) :: epsil(:) ! Local variables integer(I4B) :: i epsil=zero select case (m) case (0) do i=1,n epsil(i)=random_normal() end do case default end select return end subroutine faible FUNCTION garch(omega,alpha,beta,m,n) ! The function garch(alpha,beta,n) returns a simulation of the garch(p,q) model ! epsilon(t)=sigma(t)*eta(t) ou eta(t) bruit blanc genere par faible(epsilon,m,n) ! sigma(t)**2=omega+alpha(1)*epsilon(t-1)**2+...+alpha(q)*epsilon(t-q)**2+beta(1)*sigma(t-1)**2+...+beta(p)*sigma(t-p)**2 ! nvalint=5000 initial values are suppressed IMPLICIT NONE REAL(SP), intent(in) :: omega,alpha(:),beta(:) integer(I4B), intent(in) :: n,m REAL(SP) :: garch(1:n) ! Local variables integer(I4B), parameter :: nvalinit=5000 REAL(SP) :: eta(1:n+nvalinit),h(1:n+nvalinit), x(1:n+nvalinit) integer(I4B) :: t, i, p,q q=size(alpha) p=size(beta) call faible(eta(1:n+nvalinit),m,n+nvalinit) DO t=1,n+nvalinit h(t)=omega do i=1,q if(t-i>=1)h(t)=h(t)+alpha(i)*(x(t-i)**2) end do do i=1,p if(t-i>=1)h(t)=h(t)+beta(i)*h(t-i) end do x(t)=sqrt(h(t))*eta(t) END DO garch(1:n)=x(nvalinit+1:n+nvalinit) END FUNCTION garch subroutine estimcoviid(x,lag,rhom,sigmam,qm,pval) implicit none real(SP), intent(in) :: x(:) integer(I4B), intent(in) :: lag real(SP), intent(out) :: rhom(:),sigmam(:),qm,pval ! local variables integer(I4B) :: n,i n=size(x) qm=zero do i=1,lag rhom(i)=autocor(x(1:n),i) sigmam(i)=un qm=qm+n*(((n+2)*un)/(n-i))*(rhom(i)**2) ! qm=qm+n*(rhom(i)**2) end do pval=un-khi2(qm,lag) return end subroutine estimcoviid subroutine estimcov(x,lag,rhom,sigmam,qm,pval) implicit none real(SP), intent(in) :: x(:) integer(I4B), intent(in) :: lag real(SP), intent(out) :: rhom(:),sigmam(:),qm,pval ! local variables real(SP) :: gamma0 integer(I4B) :: n,i n=size(x) gamma0=autocov(x(1:n),0) qm=zero do i=1,lag rhom(i)=autocor(x(1:n),i) sigmam(i)=sum((x(1:n-i)*x(i+1:n))**2)/(n-i) ! sigmam(i)=autocov(x**2,i)+gamma0**2 sigmam(i)=sigmam(i)/gamma0**2 qm=qm+(rhom(i)**2)/sigmam(i) end do qm=qm*n pval=un-khi2(qm,lag) return end subroutine estimcov FUNCTION ar(a,n) ! The function ar(a) returns a simulation of the ar model ! x(t)=a(1)*x(t-1)+...+a(p)*x(t-p)+eta(t) ou eta(t) iid N(0,1) ! nvalint=500 initial values are suppressed IMPLICIT NONE REAL(SP), intent(in) :: a(:) integer(I4B), intent(in) :: n REAL(SP) :: ar(1:n) ! Local variables integer(I4B), parameter :: nvalinit=500 REAL(SP) :: eta(1:n+nvalinit), x(1:n+nvalinit) integer(I4B) :: t, i, p p=size(a) DO t=1,n+nvalinit eta(t)=random_normal() x(t)=eta(t) do i=1,p if(t-i>=1)x(t)=x(t)+a(i)*x(t-i) end do END DO ar(1:n)=x(nvalinit+1:n+nvalinit) END FUNCTION ar FUNCTION random_normal() ! The function random_normal() returns a normally distributed pseudo-random ! number with zero mean and unit variance. IMPLICIT NONE REAL(SP) :: random_normal ! Local variables REAL(SP) :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, & r1 = 0.27597, r2 = 0.27846, u, v, x, y, q REAL(SP), parameter :: half = 0.5 ! Generate P = (u,v) uniform in rectangle enclosing acceptance region DO CALL RANDOM_NUMBER(u) CALL RANDOM_NUMBER(v) v = 1.7156 * (v - half) ! Evaluate the quadratic form x = u - s y = ABS(v) - t q = x**2 + y*(a*y - b*x) ! Accept P if inside inner ellipse IF (q < r1) EXIT ! Reject P if outside outer ellipse IF (q > r2) CYCLE ! Reject P if outside acceptance region IF (v**2 < -4.0*LOG(u)*u**2) EXIT END DO ! Return ratio of P's coordinates as the normal deviate random_normal = v/u RETURN END FUNCTION random_normal FUNCTION random_exponential( ) ! Adapted from Fortran 77 code from the book: ! Dagpunar, J. 'Principles of random variate generation' ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 ! FUNCTION GENERATES A RANDOM VARIATE IN [0,INFINITY) FROM ! A NEGATIVE EXPONENTIAL DlSTRIBUTION WlTH DENSITY PROPORTIONAL ! TO EXP(-random_exponential), USING INVERSION. IMPLICIT NONE REAL (SP) :: random_exponential ! Local variable REAL (SP) :: r REAL (SP), parameter :: zero=0.0 DO CALL RANDOM_NUMBER(r) IF (r > zero) EXIT END DO random_exponential = -LOG(r) RETURN END FUNCTION random_exponential FUNCTION random_t(m) ! Adapted from Fortran 77 code from the book: ! Dagpunar, J. 'Principles of random variate generation' ! Clarendon Press, Oxford, 1988. ISBN 0-19-852202-9 ! FUNCTION GENERATES A RANDOM VARIATE FROM A ! T DISTRIBUTION USING KINDERMAN AND MONAHAN'S RATIO METHOD. ! M = DEGREES OF FREEDOM OF DISTRIBUTION ! (1 <= 1NTEGER) IMPLICIT NONE REAL(SP) :: random_t integer(I4B), INTENT(IN) :: m ! Local variables REAL (SP) :: s, c, a, f, g, r, x, v REAL (SP), parameter :: one=1.0, two=2.0, three = 3.0, four = 4.0, quart = 0.25, & five = 5.0, sixteen =16.0, zero=0.0 integer(I4B) :: mm = 0 mm=0 !SAVE s, c, a, f, g IF (m < 1) THEN write(8,*)"IMPERMISSIBLE DEGREES OF FREEDOM" STOP END IF IF (m .NE. mm) THEN ! Initialization, if necessary s = m c = -quart*(s + one) a = four/(one + one/s)**c f = sixteen/a IF (m > 1) THEN g = s - one g = ((s + one)/g)**c*SQRT((s+s)/g) ELSE g = one END IF mm = m END IF DO CALL RANDOM_NUMBER(r) IF (r <= zero) CYCLE CALL RANDOM_NUMBER(v) x = (two*v - one)*g/r v = x*x IF (v > five - a*r) THEN IF (m >= 1 .AND. r*(v + three) > f) CYCLE IF (r > (one + v/s)**c) CYCLE END IF EXIT END DO random_t = x RETURN END FUNCTION random_t subroutine chaine_Markov(p,p0,d,n,delta) implicit none ! delta(-49) est tiree selon la loi initiale p0 ! delta(-49),....,delta(n) est realisation d'une chaine de Markov ! a espace d'etats {1,...,d} et de matrice de transition p ! p(i,j)est la proba d'etre en j sachant que l'on etait en i la fois d'avant REAL(SP), intent(in) :: p(:,:), p0(:) integer(I4B), intent(out) :: delta(-49:) integer(I4B), intent(in) :: d, n ! Local variables real(SP) :: u, fdr integer(I4B) :: i, j call random_number(u) fdr=zero do i=1,d fdr=fdr+p0(i) delta(-49)=i if (u .lt. fdr) exit end do do i=-48,n call random_number(u) fdr=zero do j=1,d fdr=fdr+p(delta(i-1),j) delta(i)=j if (u .lt. fdr) exit end do end do return end subroutine chaine_Markov end module random module momgarch use constants_module; use nrtype use stat_elementaires use mod_mat use fdr implicit none private private :: durbin public :: verifmom2,verifmom4,calmom2,autocov2,autocor2,autocorpar2,autocorparemp,& portpartiid,portpart,estimcorgarch,estimcor contains subroutine estimcor(x,lag,rho) !calcule les autocorrelations rho(1),....rho(lag) de x(1),...,x(n) implicit none real(sp), intent(in) :: x(:) integer(I4b), intent(in) :: lag real(sp), intent(out) :: rho(:) ! local variables integer(I4b) :: n,i n=size(x) if(lag>=n)then write(8,*)"lag trop grand: lag,n",lag,n return end if do i=1,lag rho(i)=autocor(x(1:n),i) end do return end subroutine estimcor subroutine estimcorgarch(x,lag,rhom,sigmam) !calcule les autocorrelations rhom(1),....rhom(lag) de x(1),...,x(n) !calcule un estimateur sigmam(1:lag,1:lag) de la matrice de !variance de rhom(1:lag) pour un GARCH symétrique implicit none real(sp), intent(in) :: x(:) integer(I4b), intent(in) :: lag real(sp), intent(out) :: rhom(:),sigmam(:,:) ! local variables real(sp) :: gamma0 integer(I4b) :: n,i n=size(x) gamma0=autocov(x(1:n),0) sigmam=zero do i=1,lag rhom(i)=autocor(x(1:n),i) sigmam(i,i)=sum((x(1:n-i)*x(i+1:n))**2)/(n-i) ! sigmam(i,i)=autocov(x**2,i)+gamma0**2 sigmam(i,i)=sigmam(i,i)/gamma0**2 end do return end subroutine estimcorgarch subroutine portpart(x,r,i,sigmarm,qm,pval) implicit none real(sp), intent(in) :: x(:) real(sp), intent(out) :: r(:),sigmarm(:,:),qm,pval real(sp) :: rho(size(r)),sigmarho(size(r),size(r)) integer(I4b), intent(in) :: i ! local variables integer(I4b) :: n,lag n=size(x) lag=size(r) call estimcorgarch(x(1:n),lag,rho(1:lag),sigmarho(1:lag,1:lag)) call autocorparemp(rho(1:lag),sigmarho(1:lag,1:lag),r(1:lag),sigmarm(1:lag,1:lag)) ! qm=dot_product(sqrt(n*un)*r(1:i),matmul(inverse(sigmarm(1:i,1:i)),sqrt(n*un)*r(1:i))) qm=dot_product(sqrt(n*un)*r(1:i),matmul(inverse(sigmarho(1:i,1:i)),sqrt(n*un)*r(1:i))) pval=un-khi2(qm,i) return end subroutine portpart subroutine portpartiid(x,r,sigmarm,qm,pval) implicit none real(sp), intent(in) :: x(:) real(sp), intent(out) :: r(:),sigmarm(:,:),qm,pval real(sp) :: rho(size(r)),sigmarho(size(r),size(r)) ! local variables integer(I4b) :: n,i,lag n=size(x) lag=size(r) call estimcorgarch(x,lag,rho,sigmarho) sigmarho=identity(lag) call autocorparemp(rho(1:lag),sigmarho(1:lag,1:lag),r(1:lag),sigmarm(1:lag,1:lag)) qm=zero do i=1,lag qm=qm+n*(r(i)**2) end do pval=un-khi2(qm,lag) return end subroutine portpartiid function verifmom2(alpha,beta) real(sp), intent(in) :: alpha(:),beta(:) character (len=3) :: verifmom2 if(sum(alpha)+sum(beta).lt.un)then verifmom2='yes' else verifmom2='no' end if return end function verifmom2 function verifmom4(alpha,beta,mu4) real(sp), intent(in) :: alpha(:),beta(:),mu4 character (len=3) :: verifmom4 integer(I4b) :: i,p,q,r real(sp) :: A2((size(alpha)+size(beta))**2,(size(alpha)+size(beta))**2),& B(size(alpha)+size(beta),size(alpha)+size(beta)),& C(size(alpha)+size(beta),size(alpha)+size(beta)) real(sp) :: spectA2 q=size(alpha) p=size(beta) r=p+q if(q==0)then verifmom4='yes' return end if B=zero C=zero if(q>0)B(1,1:q)=alpha(1:q) if(p>0)B(1,q+1:r)=beta(1:p) if(p>0)C(q+1,1:r)=B(1,1:r) if(q>1)C(2:q,1:q-1)=identity(q-1) if(p>1)C(q+2:r,q+1:r-1)=identity(p-1) ! write(8,*)"B=" ! do i=1,r ! write(8,fmt='(x,20(f6.2,x))')B(i,:) ! end do ! write(8,*)"C=" ! do i=1,r ! write(8,fmt='(x,20(f6.2,x))')C(i,:) ! end do A2(1:r**2,1:r**2)=kro(B(1:r,1:r),C(1:r,1:r))+kro(C(1:r,1:r),B(1:r,1:r))+& kro(C(1:r,1:r),C(1:r,1:r))+mu4*kro(B(1:r,1:r),B(1:r,1:r)) ! write(8,*)"A2=" ! do i=1,r**2 ! write(8,fmt='(x,20(f6.2,x))')A2(i,:) ! end do spectA2=spectral_radius(A2) ! write(8,*)"rayon spectral de A2=",spectA2 ! if(p>=1.and.q>=1)then ! write(8,*)"3*alpha(1)**2+beta(1)**2+2*alpha(1)*beta(1)",3*alpha(1)**2+beta(1)**2+2*alpha(1)*beta(1) ! end if if(spectA2.lt.un)then verifmom4='yes' else verifmom4='no' end if return end function verifmom4 subroutine autocorparemp(rhom,sigmarhom,r,sigmarm) ! calcule les autocorrélations partielles empiriques et leur ! matrice de variance en fonction des autocorrelations et de la ! matrice de variance des autocorrelations real(sp), intent(in) :: rhom(:),sigmarhom(:,:) real(sp), intent(out):: r(:),sigmarm(:,:) real(sp):: J(size(rhom),size(rhom)) integer(I4b) :: lag lag=size(rhom) call durbin(rhom(1:lag),r(1:lag)) call Jdurbin(rhom(1:lag),J(1:lag,1:lag)) sigmarm=matmul(J,matmul(sigmarhom,transpose(J))) return contains subroutine Jdurbin(rho,JJ) ! algorithme de Durbin-Levinson real(sp), intent(in) :: rho(:) real(sp), intent(out) :: JJ(:,:) real(sp) :: phi(1:size(rho),1:size(rho)),derphi(1:size(rho),1:size(rho),1:size(rho)) real(sp) :: num,den,dernum,derden integer(I4b) :: m,k,i,j m=size(rho) derphi=zero phi(1,1)=rho(1) derphi(1,1,1)=un do k=2,m num=rho(k)-sum(phi(k-1,1:(k-1))*rho((k-1):1:-1)) den=un-sum(phi(k-1,1:(k-1))*rho(1:(k-1))) phi(k,k)=num/den do j=1,k dernum=-sum(rho(k-1:1:-1)*derphi(k-1,1:k-1,j)) if(j==k)dernum=dernum+un if(jcrit)then exit else if(ifirst==1)then ! write(8,*)"0:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit if(ndetect==0.or. (ndetect>0.and.minq(ndetect)>iq))then ndetect=ndetect+1 ! write(8,*)"1:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit minp(ndetect)=ip minq(ndetect)=iq iqmax=iq ! write(8,*)"A:minp(1:ndetect),minq(1:ndetect)",minp(1:ndetect),minq(1:ndetect) end if ifirst=0 else ! write(8,*)"2:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit if(ndetect==0.or. (ndetect>0.and.minq(ndetect)>iq))then ! write(8,*)"3:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit minq(ndetect)=iq iqmax=iq ! write(8,*)"B:minp(1:ndetect),minq(1:ndetect)",minp(1:ndetect),minq(1:ndetect) end if end if end if end do ! write(8,*)"4:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit end do return end subroutine detectegarch subroutine writedeltanorm(delta) ! ecrit le tableau delta normé de la methode du coin implicit none real(SP), intent(in) :: delta(:,0:) integer(I4B) :: ip,iq,lag lag=size(delta,1) write(8, fmt='(".max(p,q).|.p.",9(".",i1,"..."),16(i2,"...")/(5x,25i5))') (iq,iq=1,lag) do ip=1,lag write(8,fmt='(7x,i2," | ",24f5.1/(5x,25f5.1))' )ip,(delta(ip,iq),iq=1,lag-ip+1) end do return end subroutine writedeltanorm end subroutine coingarch subroutine detecte(table,crit,minp,minq,ndetect) !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC !C AUTOMATIC DETECTION OF THE ORDERS USING STANDARDIZED ARRAY C !C (THE ABSOLUTE VALUE OF THE ELEMENTS OF table COMPARED TO CRIT) C !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC implicit none real(SP), intent(in) :: table(:,:),crit integer(I4B), intent(out) :: minp(:),minq(:),ndetect integer(I4B) :: lag,ip,iq,ifirst,iqmax lag=size(table,1) minp=0 minq=0 ndetect=0 iqmax=lag do ip=0,lag-1 ifirst=1 do iq=min(lag-ip-1,iqmax-1),0,-1 if(maxval(abs(table(ip+1:lag,iq+1:lag)))>crit)then exit else if(ifirst==1)then ! write(8,*)"0:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit if(ndetect==0.or. (ndetect>0.and.minq(ndetect)>iq))then ndetect=ndetect+1 ! write(8,*)"1:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit minp(ndetect)=ip minq(ndetect)=iq iqmax=iq ! write(8,*)"A:minp(1:ndetect),minq(1:ndetect)",minp(1:ndetect),minq(1:ndetect) end if ifirst=0 else ! write(8,*)"2:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit if(ndetect==0.or. (ndetect>0.and.minq(ndetect)>iq))then ! write(8,*)"3:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit minq(ndetect)=iq iqmax=iq ! write(8,*)"B:minp(1:ndetect),minq(1:ndetect)",minp(1:ndetect),minq(1:ndetect) end if end if end if end do ! write(8,*)"4:ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit",ip,iq,maxval(abs(table(ip+1:lag,iq+1:lag))),crit end do return end subroutine detecte subroutine dercorner(rho,derdelta) implicit none real(SP), intent(in) :: rho(:) real(SP), intent(out) :: derdelta(:,0:,:) real(SP) :: delta(1:size(rho),0:size(rho)),ddelta(-size(rho)-1:size(rho)+1,0:size(rho)+1),& derddelta(-size(rho)-1:size(rho)+1,0:size(rho)+1,1:size(rho)),num,den,& dernum(1:size(rho)),derden(1:size(rho)) integer(I4B) :: i,j,lag lag=size(rho) ddelta=zero derddelta=zero ddelta(-lag-1:lag+1,0)=un ddelta(0,1)=un ddelta(1:lag,1)=rho(1:lag) ddelta(-lag:-1,1)=rho(lag:1:-1) do i=1,lag derddelta(i,1,i)=un derddelta(-i,1,i)=un end do do j=1,lag do i=-lag+j,lag-j,1 num=ddelta(i,j)**2-ddelta(i+1,j)*ddelta(i-1,j) dernum(1:lag)=2*ddelta(i,j)*derddelta(i,j,1:lag)-ddelta(i+1,j)*derddelta(i-1,j,1:lag)-& derddelta(i+1,j,1:lag)*ddelta(i-1,j) den=ddelta(i,j-1) derden(1:lag)=derddelta(i,j-1,1:lag) ddelta(i,j+1)=num/den derddelta(i,j+1,1:lag)=(den*dernum(1:lag)-num*derden(1:lag))/den**2 end do end do do i=1,lag do j=0,lag delta(i,j)=ddelta(j,i) derdelta(i,j,1:lag)=derddelta(j,i,1:lag) end do end do return end subroutine dercorner subroutine writedeltanorm(delta) ! ecrit le tableau delta normé de la methode du coin implicit none real(SP), intent(in) :: delta(:,0:) integer(I4B) :: ip,iq,lag lag=size(delta,1) write(8, fmt='(" .p.|.q.",9(".",i1,"..."),16(i2,"...")/(5x,25i5))') (iq,iq=1,lag) do ip=1,lag write(8,fmt='(2x,i2," | ",24f5.1/(5x,25f5.1))' )ip,(delta(ip,iq),iq=1,lag-ip+1) end do return end subroutine writedeltanorm subroutine autocorviaar(x,rho) !********************************************************************** ! estime les autocorrelations par approximation AR ! !********************************************************************** implicit none real(SP), intent(in) :: x(:) ! output parameters ! real(SP), intent(out) :: rho(:) ! Local variables integer(I4b) :: ordremax,ordre real(SP) :: b(0),a0(min(30,int(size(x)/2))),gamma(0:min(30,int(size(x)/2))) integer(I4b) :: i,q,n,lag,lagsuf ordremax=min(30,int(size(x)/2)) n=size(x) lag=size(rho) do i=0,ordremax gamma(i)=autocov(x(1:n),i) end do call durb(n,gamma(0:ordremax),a0(1:ordremax),ordre) q=0 call autocorarma(a0(1:ordre),b(1:q),rho(1:lag),lagsuf) return end subroutine autocorviaar subroutine estsigmarho(x,r,rmax,sigmarho,rho,aicbic) ! estime \Sigma_{\rho} pour lineaire faible par regressions auxiliaires AR(r) real(SP), intent(in) :: x(:) integer(I4b), intent(out) :: r integer(I4b), intent(in) :: rmax,aicbic real(SP), intent(out) :: sigmarho(:,:),rho(:) real(SP) :: moyupsilon(0:size(sigmarho,1)) real(SP) :: sigmagamma(0:size(sigmarho,1),0:size(sigmarho,1)) real(SP) :: upsilon(size(x)-size(sigmarho,1),0:size(sigmarho,1)) real(SP) :: upsiloncentre(size(x)-size(sigmarho,1),0:size(sigmarho,1)) real(sp):: J(size(rho),size(rho)+1),gamma0 integer(I4b) :: lag,n,i n=size(x) lag=size(sigmarho,1) moyupsilon(0:lag)=zero do i=1,n-lag upsilon(i,0:lag)=x(i)*x(i:i+lag) moyupsilon(0:lag)=moyupsilon(0:lag)+upsilon(i,0:lag)/(n-lag) end do do i=1,n-lag upsiloncentre(i,0:lag)=upsilon(i,0:lag)-moyupsilon(0:lag) end do call determinep(upsiloncentre(1:(n-lag),0:lag),r,rmax,aicbic) if(r >rmax)write(*,*)"rmax trop petit" call estR(upsiloncentre(1:(n-lag),0:lag),r,sigmagamma(0:lag,0:lag)) call estimcor(x(1:n),lag,rho(1:lag)) gamma0=autocov(x(1:n),0) J=zero do i=1,lag J(i,1)=-rho(i)/gamma0 J(i,i+1)=un/gamma0 end do sigmarho(1:lag,1:lag)=matmul(J(1:lag,1:lag+1),matmul(sigmagamma(0:lag,0:lag),transpose(J(1:lag,1:lag+1)))) return end subroutine estsigmarho subroutine estR(upsilon,p,R) ! estime R_{m,m} en utilisant une regression auxiliaire real(SP), intent(in) :: upsilon(:,:) integer(I4b), intent(in) :: p real(SP), intent(out) :: R(:,:) real(SP) :: phi(size(upsilon,2),size(upsilon,2)),invphi(size(upsilon,2),size(upsilon,2)),& s2(size(upsilon,2),size(upsilon,2)),a(p,size(upsilon,2),size(upsilon,2)) integer(I4b) :: n,i,j,m n=size(upsilon,1) m=size(upsilon,2) if (p>0) then call estarmultmod(upsilon(1:n,1:m),a,s2) phi(1:m,1:m)=identity(m) do j=1,p phi(1:m,1:m)=phi(1:m,1:m)-a(j,1:m,1:m) end do invphi(1:m,1:m)=inverse(phi(1:m,1:m)) R(1:m,1:m)=matmul(matmul(invphi(1:m,1:m),s2(1:m,1:m)),transpose(invphi(1:m,1:m))) else R(1:m,1:m)=zero do i=1,n R(1:m,1:m)=R(1:m,1:m)+prod(upsilon(i,1:m),upsilon(i,1:m)) end do R(1:m,1:m)=R(1:m,1:m)/n end if return end subroutine estR subroutine estarmultmod(x,a,s2) ! estime le modele d-multivarie x(t)=a(1)*x(t-1)+...+a(p)*x(t-p)+e(t) ou var(e(t))=s2 ! algorithme de Whittle (multivariate Durbin-Levinson algorithm) real(SP), intent(in) :: x(1:,1:) real(SP), intent(out) :: a(:,:,:), s2(:,:) real(SP) :: Gamma(0:size(a,1),size(x,2),size(x,2)) real(SP) :: Delta(0:size(a,1),size(x,2),size(x,2)),Deltatilde(0:size(a,1),size(x,2),size(x,2)) real(SP) :: V(0:size(a,1),size(x,2),size(x,2)),Vtilde(0:size(a,1),size(x,2),size(x,2)) real(SP) :: Phi(0:size(a,1),0:size(a,1),size(x,2),size(x,2)) real(SP) :: Phitilde(0:size(a,1),0:size(a,1),size(x,2),size(x,2)) integer(I4b) :: k,nn,d,p,j d=size(x,2) p=size(a,1) do j=0,p Gamma(j,1:d,1:d)=autocovmult(x,j) end do ! write(8,*)Gamma V(0,1:d,1:d)=Gamma(0,1:d,1:d) Vtilde(0,1:d,1:d)=V(0,1:d,1:d) Delta(0,1:d,1:d)=Gamma(1,1:d,1:d) Deltatilde(0,1:d,1:d)=transpose(Delta(0,1:d,1:d)) do nn=1,p Phi(nn,nn,1:d,1:d)=matmul(Delta(nn-1,1:d,1:d),inverse(Vtilde(nn-1,1:d,1:d))) Phitilde(nn,nn,1:d,1:d)=matmul(Deltatilde(nn-1,1:d,1:d),inverse(V(nn-1,1:d,1:d))) do k=1,nn-1 Phi(nn,k,1:d,1:d)=Phi(nn-1,k,1:d,1:d)-matmul(Phi(nn,nn,1:d,1:d),Phitilde(nn-1,nn-k,1:d,1:d)) Phitilde(nn,k,1:d,1:d)=Phitilde(nn-1,k,1:d,1:d)-matmul(Phitilde(nn,nn,1:d,1:d),Phi(nn-1,nn-k,1:d,1:d)) end do if(nn= 0" stop end if write(8,*)"Pour une realisation de taille n= ",n write (8,*) write (8,*) write(8,*)"Lag borne inf centre borne sup" write(8,*)" " do i=1,lag binf=rho(i)-1.96*sqrt(var(i,i)/n) bsup=rho(i)+1.96*sqrt(var(i,i)/n) write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,binf,rho(i),bsup end do write(8,fmt='(1x,40(2x,f7.3,a1))')((rho(i)-1.96*sqrt(var(i,i)/n),virgul),i=1,lag) write(8,fmt='(1x,40(2x,f7.3,a1))')((rho(i),virgul),i=1,lag) write(8,fmt='(1x,40(2x,f7.3,a1))')((rho(i)+1.96*sqrt(var(i,i)/n),virgul),i=1,lag) write (8,*) write (8,*) write (8,*) " ____________________________________________________________________________________________" write (8,*) " Intervalle de confiance a 95% pour les autocorrelations partielles empiriques " write (8,*) " Formule de bartlett generalisee " write (8,*) " ____________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Pour une realisation de taille n= ",n call autocorparemp(rho(1:lag),var(1:lag,1:lag),r(1:lag),varr(1:lag,1:lag)) ! write(8,*)"Fonction d'autocorrelation partielle de X:" ! write(8,fmt='(1x,48(f7.3,1x))')r(1:lag) ! write (8,*) write (8,*) write (8,*) write(8,*)"Lag borne inf centre borne sup" write(8,*)" " do i=1,lag binf=r(i)-1.96*sqrt(varr(i,i)/n) bsup=r(i)+1.96*sqrt(varr(i,i)/n) write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,binf,r(i),bsup end do write(8,fmt='(1x,40(2x,f7.3,a1))')((r(i)-1.96*sqrt(varr(i,i)/n),virgul),i=1,lag) write(8,fmt='(1x,40(2x,f7.3,a1))')((r(i),virgul),i=1,lag) write(8,fmt='(1x,40(2x,f7.3,a1))')((r(i)+1.96*sqrt(varr(i,i)/n),virgul),i=1,lag) return end subroutine momtheo subroutine statioinver(a,b,statio,inver) use nrtype real (SP), intent(in), dimension(:) :: a,b character (len=3), intent(out) :: statio,inver integer(I4B) :: p,q,p0,q0,i complex(sp) :: rootsa(size(a)),rootsb(size(b)) complex(sp) :: aa(size(a)+1),bb(size(b)+1) logical(lgt) :: polish polish=.true. p=size(a) q=size(b) p0=p do i=1,p if(a(p-i+1)==zero)then p0=p-1 else exit end if end do q0=q do i=1,q if(b(q-i+1)==zero)then q0=q-1 else exit end if end do if(p0>=1)then aa(1)=cmplx(un) aa(2:p0+1)=-cmplx(a(1:p0)) call zroots(aa(1:p0+1),rootsa(1:p0),polish) write(8,*)"racines du polynome AR:",rootsa if(minval(abs(rootsa(1:p0)))>un)then statio='yes' else statio='no' end if else statio='yes' end if if(q0>=1)then bb(1)=cmplx(un) bb(2:q0+1)=-cmplx(b(1:q0)) call zroots(bb(1:q0+1),rootsb(1:q0),polish) write(8,*)"racines du polynome MA:",rootsb if(minval(abs(rootsb(1:q0)))>un)then inver='yes' else inver='no' end if else inver='yes' end if return end subroutine statioinver subroutine readpara(omega,mu4,lag,sigma) !lit les parametres ARMA-GARCH implicit none real(SP), intent(out) :: omega,mu4,sigma integer(I4B), intent(out) :: lag integer(I4B) :: i,p,q,pp,qq character (len=3) :: statio, inver character (len=5), parameter :: xdet=' X(t)',Zdet=' Z(t)',hdet='h(t)=',Z2det='Z^2(t' character (len=4), parameter :: xde=' X(t',hde=' h(t',Zde=' Z(t',etc='+...' character (len=1), parameter :: par=')',minus='-',plus='+',virgul=',' character (len=5), parameter :: egalZdet='=Z(t)',egalXdet='=X(t)' ! lecture des parametres du modele ARMA(p,q) write(*,*)"Modele ARMA(p,q)" write(*,*)"X(t)-a(1)X(t-1)-....-a(p)X(t-1)=Z(t)-b(1)-.......-b(q)Z(t-q)" write(*,*)"p=?" read(7,*) p if(p<0)then write(8,*)"p doit etre positif ou nul" stop end if allocate(ar0(p)) do i=1,p write(*,*)"a(",i,")=?" read(7,*)ar0(i) end do write(*,*)"q=?" read(7,*) q if(q<0)then write(8,*)"q doit etre positif ou nul" stop end if allocate(bma0(q)) do i=1,q write(*,*)"b(",i,")=?" read(7,*)bma0(i) end do write(*,*)"nombre d'autocorrélations a calculer ?" read(7,*) lag if(lag<=0)then write(8,*)"ce nombre doit etre positif" stop end if ! ! lecture des parametres du modele GARCH(p,q) write(*,*)"Modele GARCH(P,Q)" write(*,*)"Z(t)=h(t)^{1/2}e(t), e(t) bb iid centre reduit symetrique" write(*,*)"h(t)=omega+alpha(1)*Z(t-1)^2+...+alpha(Q)*Z(t-Q)^2+beta(1)*h(t-1)+...+beta(P)*h(t-P)" write(*,*)"Ee(t)^4=?" read(7,*)mu4 if(mu4=1" stop end if write(*,*)"omega=?" read(7,*)omega if(omega<=zero)then write(8,*)"ce nombre doit etre >0" stop end if write(*,*)"Q=" read(7,*) qq if(qq<=0)then write(8,*)"Q doit etre > 0" stop end if allocate(alpha0(qq)) do i=1,qq write(*,*)"alpha(",i,")=?" read(7,*)alpha0(i) if(alpha0(i)= 0" stop end if end do write(*,*)"P=" read(7,*) pp if(pp<0)then write(8,*)"P doit etre >= 0" stop end if allocate(beta0(pp)) do i=1,pp write(*,*)"beta(",i,")=?" read(7,*)beta0(i) if(beta0(i)= 0" stop end if end do ! ! ecriture des parametres du modele ARMA(p,q)-GARCH(pp,qq) write (8,*) write (8,*) write (8,*) " ______________________________________________" write (8,*) " Modele ARMA-GARCH " write (8,*) " ______________________________________________" write (8,*) write (8,*) write (8,*) " Paramètres ARMA" write (8,*) write (8,*) ! write(8,fmt='(10x,a5,10(sp,f6.3,a4,a1,ss,i1,a1))')xdet,((-un*ar0(i),xde,minus,i,par),i=1,p) write(8,fmt='(10x,a5,10(sp,f6.3,a4,a1,ss,i1,a1))')egalZdet,((-un*bma0(i),Zde,minus,i,par),i=1,q) write (8,*) sigma=sqrt(omega/(un-sum(alpha0(1:qq))-sum(beta0(1:pp)))) write(8,*)"La variance de Z(t) vaut", sigma**2 write (8,*) ! vérification des condition de stationnarité et inversibilite call statioinver(ar0(1:p),bma0(1:q),statio,inver) ! write(8,*)"polynome AR inversible : ",statio ! write(8,*)"polynome MA inversible : ",inver if(statio=='no')then write(8,*)"le polynome AR n'est pas inversible " stop end if if(inver=='no')then write(8,*)"le polynome MA n'est pas inversible " stop end if write (8,*) write (8,*) write(8,*) " Paramètres GARCH" write (8,*) write(8,*) " Z(t)=h(t)^{1/2}eta(t), eta(t) iid, Ee(t)=Ee(t)^3=0, Ee(t)^2=1, Ee(t)^4=",mu4 write (8,*) write(8,fmt='(10x,a5,sp,f6.3,10(sp,f6.3,a5,a1,ss,i1,a1))')hdet,omega,((alpha0(i),Z2det,minus,i,par),i=1,qq) write(8,fmt='(15x,10(sp,f6.3,a4,a1,ss,i1,a1))')((beta0(i),hde,minus,i,par),i=1,pp) return end subroutine readpara subroutine estbart(x,gamma,var,lagmax,ordre,epsil,rho,lagsuf) !********************************************************************** ! ! var= estimation de la matrice de variance asymptotique des lag premieres ! autocorrelations empiriques d'un modele lineaire fort ! implicit none real(SP), intent(in) :: x(:),gamma(0:) integer(I4b), intent(in) :: lagmax ! output parameters ! integer(I4b), intent(out) :: ordre,lagsuf real(SP), intent(out) :: var(:,:),epsil(:),rho(0:) ! Local variables real(SP) :: b(0),a0(size(gamma)-1) integer(I4b) :: i,lag,n,q,ordremax lag=size(var,1) n=size(x) ordremax=size(gamma)-1 call durb(n,gamma(0:ordremax),a0(1:ordremax),ordre) q=0 call bartlett(a0(1:ordre),b(1:q),var(1:lag,1:lag),rho(0:lagmax),lagsuf) if(lagsuf>lagmax)then write(8,*)"Il Faudrait prendre lagmax>",lagsuf return end if do i=ordre+1,n epsil(i)=x(i)-sum(a0(1:ordre)*x(i-1:i-ordre:-1)) end do return end subroutine estbart subroutine durb(n,gamma,a0,ordre) ! estime le modele x(t)=a(1)*x(t-1)+...+a(p)*x(t-p)+e(t) ou var(e(t))=s2 ! algorithme de Durbin-Levinson real(SP), intent(in) :: gamma(0:) real(SP), intent(out) :: a0(:) integer(I4b), intent(in) :: n integer(I4b), intent(out) :: ordre real(SP) :: v(1:size(gamma)-1),critmin,crit real(SP) :: phi(1:size(gamma)-1,1:size(gamma)-1) integer(I4b) :: m,ordremax ordremax=size(gamma)-1 critmin=log(gamma(0)) ! write(8,*)"critmin",critmin,"v(0)",gamma(0) ordre=0 if(ordremax>0)then phi(1,1)=gamma(1)/gamma(0) v(1)=gamma(0)*(un-phi(1,1)**2) ! crit=log(v(1))+(deux)/n crit=log(v(1))+(log(un*n))/n ! write(8,*)"Step",1,"crit",crit,"v(1)",v(1) if(crit=n)then write(8,*)"ordremax trop grand: ordremax,n=",ordremax,n return end if do i=0,ordremax gamma(i)=autocov(x(1:n),i) end do call estbart(x(1:n),gamma(0:ordremax),var2(1:lag,1:lag),lagmax,ordre,epsil(1:n),rho(0:lagmax),lagsuf) do i=1,lag rrho(i)=autocor(x(1:n),i) end do if(ordremax>=(n-ordre+1))then write(8,*)"ordremax trop grand: ordremax,n-ordre+1=",ordremax,n-ordre+1 return end if if(method==1)then var=var2 ordre2=0 return end if neps2=n-ordre epsil2(1:neps2)=epsil(ordre+1:n)**2 do i=0,ordremax gamma2(i)=autocov(epsil2(1:neps2),i) end do call durb(n-ordre,gamma2(0:ordremax),a20(1:ordremax),ordre2) gameps0=autocov(epsil(ordre+1:n),0) fact=gamma2(0)/(gameps0**2) if(lagmaxlagmax)then write(8,*)"Il faudrait prendre lagmax>",lagsuf return end if lagutil=lagmax if(lagsuf=-lagmax))var(i,j)=var(i,j)+rho2(abs(l))*rho(abs(l+i))*rho(abs(l-j)) end do if(i.ne.j)var(j,i)=var(i,j) end do end do var=fact*var+var2 return end subroutine estbartlettgene subroutine estbartlett2(x,var,lagmax,ordremax,ordre,epsil,lagsuf) !********************************************************************** ! ! var= estimation de la matrice de variance asymptotique des lag premieres ! autocorrelations empiriques d'un modele lineaire fort ! implicit none real(SP), intent(in) :: x(:) integer(I4b), intent(in) :: lagmax,ordremax ! output parameters ! integer(I4b), intent(out) :: ordre,lagsuf real(SP), intent(out) :: var(:,:),epsil(:) ! Local variables real(SP) :: b(0),a0(ordremax),gamma(0:ordremax),rho0(0:lagmax) integer(I4b) :: i,lag,n,q lag=size(var,1) n=size(x) if(ordremax>=n)then write(8,*)"ordremax trop grand: ordremax,n=",ordremax,n return end if do i=0,ordremax gamma(i)=autocov(x(1:n),i) end do call durb(n,gamma(0:ordremax),a0(1:ordremax),ordre) q=0 call bartlett(a0(1:ordre),b(1:q),var(1:lag,1:lag),rho0(0:lagmax),lagsuf) do i=ordre+1,n epsil(i)=x(i)-sum(a0(1:ordre)*x(i-1:i-ordre:-1)) end do write(8,*)"moyenne des carres des residus=",sum(epsil(ordre+1:n)**2)/(n-ordre+1) return end subroutine estbartlett2 subroutine estbartlett(x,var,lagmax) !********************************************************************** ! ! var= estimation de la matrice de variance asymptotique des lag premieres ! autocorrelations empiriques d'un modele lineaire fort ! implicit none real(SP), intent(in) :: x(:) integer(I4b), intent(in) :: lagmax ! output parameters ! real(SP), intent(out) :: var(:,:) ! Local variables real(SP) :: rho(0:lagmax) integer(I4b) :: i,j,lag,l,n lag=size(var,1) n=size(x) call estimcor(x(1:n),lagmax,rho(1:lagmax)) rho(0)=un if(lagmax=-lagmax))var(i,j)=var(i,j)+ rho(abs(l))*rho(abs(l+j-i)) if(l-j-i>=-lagmax)var(i,j)=var(i,j)+rho(abs(l))*rho(abs(l-j-i)) end do if(i.ne.j)var(j,i)=var(i,j) end do end do return end subroutine estbartlett subroutine bartlettgene(a,b,omega,alpha,beta,mu4,var,lagmax,lagsuf1,lagsuf2) !********************************************************************** ! ! var= matrice de variance asymptotique des lag premieres ! autocorrelations empiriques d'un modele lineaire avec bb GARCH ! implicit none real(SP), intent(in) :: a(:),b(:),omega,alpha(:),beta(:),mu4 integer(I4b), intent(in) :: lagmax ! output parameters ! integer(I4b), intent(out) :: lagsuf1,lagsuf2 real(SP), intent(out) :: var(:,:) ! Local variables real(SP) :: rho(0:lagmax),rho2(0:lagmax),var1(size(var,1),size(var,1)),fact,autocoveps2l0,autocovepsl0 integer(I4b) :: p,q,pp,qq,i,j,lag,l,lagutil pp=size(beta) qq=size(alpha) p=size(a) q=size(b) lag=size(var,1) lagsuf1=lagmax+1 lagsuf2=lagmax+1 call bartlett(a(1:p),b(1:q),var1(1:lag,1:lag),rho(0:lagmax),lagsuf1) call autocor2(rho2(1:lagmax),omega,alpha(1:qq),beta(1:pp),mu4,autocoveps2l0,lagsuf2) rho2(0)=un autocovepsl0=omega/(un-sum(alpha(1:qq))-sum(beta(1:pp))) fact=autocoveps2l0/(autocovepsl0**2) if(lagmax=-lagmax))var(i,j)=var(i,j)+rho2(abs(l))*rho(abs(l+i))*rho(abs(l-j)) end do if(i.ne.j)var(j,i)=var(i,j) end do end do var=fact*var+var1 return end subroutine bartlettgene subroutine bartlett(a,b,var,rho,lagsuf) !********************************************************************** ! ! var= matrice de variance asymptotique des lag premieres ! autocorrelations empiriques d'un modele lineaire fort ! implicit none real(SP), intent(in) :: a(:),b(:) ! output parameters ! real(SP), intent(out) :: var(:,:),rho(0:) integer(I4b), intent(out) :: lagsuf ! Local variables integer(I4b) :: lagmax integer(I4b) :: i,j,lag,l,p,q lag=size(var,1) p=size(a) q=size(b) lagmax=size(rho)-1 call autocorarma(a(1:p),b(1:q),rho(1:lagmax),lagsuf) rho(0)=un if(lagmax=-lagmax))var(i,j)=var(i,j)+ rho(abs(l))*rho(abs(l+j-i)) if(l-j-i>=-lagmax)var(i,j)=var(i,j)+rho(abs(l))*rho(abs(l-j-i)) end do if(i.ne.j)var(j,i)=var(i,j) end do end do return end subroutine bartlett subroutine autocorpar(rhom,r) ! calcule les autocorrélations partielles empiriques en fonction des autocorrelations real(sp), intent(in) :: rhom(:) real(sp), intent(out):: r(:) integer(I4B) :: lag lag=size(rhom) call durbin(rhom(1:lag),r(1:lag)) return contains subroutine durbin(rho2,r2) ! algorithme de Durbin-Levinson real(sp), intent(in) :: rho2(:) real(sp), intent(out) :: r2(:) real(sp) :: phi(1:size(r2),1:size(r2)) integer(I4B) :: h,k h=size(r2) phi(1,1)=rho2(1) r2(1)=phi(1,1) do k=2,h phi(k,k)=(rho2(k)-sum(phi(k-1,1:(k-1))*rho2((k-1):1:-1)))/(un-sum(phi(k-1,1:(k-1))*rho2(1:(k-1)))) r2(k)=phi(k,k) phi(k,1:k-1)=phi(k-1,1:k-1)-phi(k,k)*phi(k-1,(k-1):1:-1) end do return end subroutine durbin end subroutine autocorpar subroutine autocorarma(a,b,rho,lagsuf) !********************************************************************** ! ! input parameters ! lag - number of lags to be computed for the autocorrelation function ! a and b - model specification parameters ! z(t)-a(1)z(t-1)-...-a(p)z(t-p)=epsilon(t)-b(1)epsilon(t-1)...-b(q)epsilon(t-q) ! sigma**2=var(epsilon(t)) implicit none real(SP), intent(in) :: a(:),b(:) ! output parameters ! rho(h)=cov(z(t),z(t-h))/var(z(t)) real(SP), intent(out) :: rho(:) integer(I4b), intent(out) :: lagsuf ! Local variables real(SP) :: gamma(0:size(rho)) integer(I4b) :: lag,p,q lag=size(rho) lagsuf=lag+1 p=size(a) q=size(b) call autocovarma(a(1:p),b(1:q),un,gamma(0:lag),lagsuf) if(gamma(0)>zero)then rho(1:lag)=gamma(1:lag)/gamma(0) else write(8,*)"pb dans autocorarma: gamma(0)=",gamma(0) end if return end subroutine autocorarma subroutine autocovarma(a,b,sigma,gamma,lagsuf) !********************************************************************** ! ! input parameters ! lag - number of lags to be computed for the autocovariance function ! a and b - model specification parameters ! z(t)-a(1)z(t-1)-...-a(p)z(t-p)=epsilon(t)-b(1)epsilon(t-1)...-b(q)epsilon(t-q) ! sigma**2=var(epsilon(t)) implicit none real(SP), intent(in) :: a(:),b(:),sigma ! output parameters ! gamma(h)=cov(z(t),z(t-h)) real(SP), intent(out) :: gamma(0:) integer(I4b), intent(out) :: lagsuf ! Local variables REAL(SP), PARAMETER :: TINY=1.0e-5_sp real(SP) :: M(size(a)+1,size(a)+1),c0P(0:(size(a)+size(b))),dum(0:size(a),1) real(SP) :: theta(0:size(b)),phi(0:size(a)),psi(0:size(b)) integer(I4b) :: i,k,j,p,q,lag lag=size(gamma)-1 lagsuf=lag+1 p=size(a) q=size(b) theta(0)=un phi(0)=-un theta(1:q)=-b(1:q) phi(1:p)=a(1:p) M=zero do k=0,p do j=0,p M(k+1,abs(k-j)+1)=M(k+1,abs(k-j)+1)-phi(j) end do end do call marep(psi(0:q),a(1:p),b(1:q)) do k=0,q c0P(k)=(sigma**2)*sum(theta(k:q)*psi(0:q-k)) end do do k=q+1,p c0P(k)=zero end do dum(0:p,1)=c0P(0:p) call gaussj(M(1:p+1,1:p+1),dum(0:p,:)) gamma(0:p)=dum(0:p,1) do k=p+1,q gamma(k)=sum(phi(1:p)*gamma(k-1:k-p:-1))+c0P(k) end do do k=max(p,q+1),lag gamma(k)=dot_product(phi(1:p),gamma(k-1:k-p:-1)) end do do i=lag,0,-1 if(abs(gamma(i))>tiny)then lagsuf=i+1 exit end if end do return end subroutine autocovarma subroutine marep(psi, a,b) !********************************************************************** ! ! input parameters ! lmax - number of lags to be used in ma approximation ! a and b - model specification parameters ! z(t)-a(1)z(t-1)-...-a(p)z(t-p)=epsilon(t)-b(1)epsilon(t-1)...-b(q)epsilon(t-q) ! implicit none real(SP), intent(in) :: a(:), b(:) ! output parameters ! psi - moving average approximation of model ! z(t) = e(t) + psi(1)*e(t-1) +...+ psi(lmax)*e(t-lmax) real(SP), intent(out) :: psi(0:) ! Local variables real(SP) :: theta(0:size(b)),phi(0:size(a)) integer(I4b) :: j,lmax,p,q p=size(a) q=size(b) lmax=size(psi)-1 if(lmax==0)then psi(0)=un return end if theta(0)=un phi(0)=un theta(1:q)=-b(1:q) phi(1:p)=a(1:p) do j=0,q psi(j)=theta(j)+sum(phi(1:min(j,p))*psi(j-1:j-min(j,p):-1)) end do do j=q+1,max(p,q+1)-1 psi(j)=sum(phi(1:min(j,p))*psi(j-1:j-min(j,p):-1)) end do do j=max(p,q+1),lmax psi(j)=sum(phi(1:p)*psi(j-1:j-p:-1)) end do return end subroutine marep SUBROUTINE laguer(a,x,its) USE nrtype; USE nrutil, ONLY : nrerror,poly,poly_term IMPLICIT NONE INTEGER(I4B), INTENT(OUT) :: its COMPLEX(SPC), INTENT(INOUT) :: x COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a REAL(SP), PARAMETER :: EPS=epsilon(1.0_sp) INTEGER(I4B), PARAMETER :: MR=8,MT=10,MAXIT=MT*MR INTEGER(I4B) :: iter,m REAL(SP) :: abx,abp,abm,err COMPLEX(SPC) :: dx,x1,f,g,h,sq,gp,gm,g2 COMPLEX(SPC), DIMENSION(size(a)) :: b,d REAL(SP), DIMENSION(MR) :: frac = & (/ 0.5_sp,0.25_sp,0.75_sp,0.13_sp,0.38_sp,0.62_sp,0.88_sp,1.0_sp /) m=size(a)-1 do iter=1,MAXIT its=iter abx=abs(x) b(m+1:1:-1)=poly_term(a(m+1:1:-1),x) d(m:1:-1)=poly_term(b(m+1:2:-1),x) f=poly(x,d(2:m)) err=EPS*poly(abx,abs(b(1:m+1))) if (abs(b(1)) <= err) RETURN g=d(1)/b(1) g2=g*g h=g2-2.0_sp*f/b(1) sq=sqrt((m-1)*(m*h-g2)) gp=g+sq gm=g-sq abp=abs(gp) abm=abs(gm) if (abp < abm) gp=gm if (max(abp,abm) > 0.0) then dx=m/gp else dx=exp(cmplx(log(1.0_sp+abx),iter,kind=spc)) end if x1=x-dx if (x == x1) RETURN if (mod(iter,MT) /= 0) then x=x1 else x=x-dx*frac(iter/MT) end if end do call nrerror('laguer: too many iterations') END SUBROUTINE laguer SUBROUTINE zroots(a,roots,polish) USE nrtype; USE nrutil, ONLY : assert_eq,poly_term ! USE nr, ONLY : indexx IMPLICIT NONE COMPLEX(SPC), DIMENSION(:), INTENT(IN) :: a COMPLEX(SPC), DIMENSION(:), INTENT(OUT) :: roots LOGICAL(LGT), INTENT(IN) :: polish REAL(SP), PARAMETER :: EPS=1.0e-6_sp INTEGER(I4B) :: j,its,m INTEGER(I4B), DIMENSION(size(roots)) :: indx COMPLEX(SPC) :: x COMPLEX(SPC), DIMENSION(size(a)) :: ad m=assert_eq(size(roots),size(a)-1,'zroots') ad(:)=a(:) do j=m,1,-1 x=cmplx(0.0_sp,kind=spc) call laguer(ad(1:j+1),x,its) if (abs(aimag(x)) <= 2.0_sp*EPS**2*abs(real(x))) & x=cmplx(real(x),kind=spc) roots(j)=x ad(j:1:-1)=poly_term(ad(j+1:2:-1),x) end do if (polish) then do j=1,m call laguer(a(:),roots(j),its) end do end if call indexx_sp(real(roots),indx) roots=roots(indx) return CONTAINS SUBROUTINE indexx_sp(arr,index) USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap IMPLICIT NONE REAL(SP), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50 REAL(SP) :: a INTEGER(I4B) :: n,k,i,j,indext,jstack,l,r INTEGER(I4B), DIMENSION(NSTACK) :: istack n=assert_eq(size(index),size(arr),'indexx_sp') index=arth(1,1,n) jstack=0 l=1 r=n do if (r-l < NN) then do j=l+1,r indext=index(j) a=arr(indext) do i=j-1,l,-1 if (arr(index(i)) <= a) exit index(i+1)=index(i) end do index(i+1)=indext end do if (jstack == 0) RETURN r=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+r)/2 call swap(index(k),index(l+1)) call icomp_xchg_sp(arr,index(l),index(r)) call icomp_xchg_sp(arr,index(l+1),index(r)) call icomp_xchg_sp(arr,index(l),index(l+1)) i=l+1 j=r indext=index(l+1) a=arr(indext) do do i=i+1 if (arr(index(i)) >= a) exit end do do j=j-1 if (arr(index(j)) <= a) exit end do if (j < i) exit call swap(index(i),index(j)) end do index(l+1)=index(j) index(j)=indext jstack=jstack+2 if (jstack > NSTACK) call nrerror('indexx: NSTACK too small') if (r-i+1 >= j-l) then istack(jstack)=r istack(jstack-1)=i r=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i end if end if end do END SUBROUTINE indexx_sp !BL SUBROUTINE icomp_xchg_sp(arr,i,j) USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap IMPLICIT NONE REAL(SP), DIMENSION(:), INTENT(IN) :: arr INTEGER(I4B), INTENT(INOUT) :: i,j INTEGER(I4B) :: swp if (arr(j) < arr(i)) then swp=i i=j j=swp end if END SUBROUTINE icomp_xchg_sp SUBROUTINE indexx_i4b(iarr,index) USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap IMPLICIT NONE INTEGER(I4B), DIMENSION(:), INTENT(IN) :: iarr INTEGER(I4B), DIMENSION(:), INTENT(OUT) :: index INTEGER(I4B), PARAMETER :: NN=15, NSTACK=50 INTEGER(I4B) :: a INTEGER(I4B) :: n,k,i,j,indext,jstack,l,r INTEGER(I4B), DIMENSION(NSTACK) :: istack n=assert_eq(size(index),size(iarr),'indexx_sp') index=arth(1,1,n) jstack=0 l=1 r=n do if (r-l < NN) then do j=l+1,r indext=index(j) a=iarr(indext) do i=j-1,l,-1 if (iarr(index(i)) <= a) exit index(i+1)=index(i) end do index(i+1)=indext end do if (jstack == 0) RETURN r=istack(jstack) l=istack(jstack-1) jstack=jstack-2 else k=(l+r)/2 call swap(index(k),index(l+1)) call icomp_xchg(iarr,index(l),index(r)) call icomp_xchg(iarr,index(l+1),index(r)) call icomp_xchg(iarr,index(l),index(l+1)) i=l+1 j=r indext=index(l+1) a=iarr(indext) do do i=i+1 if (iarr(index(i)) >= a) exit end do do j=j-1 if (iarr(index(j)) <= a) exit end do if (j < i) exit call swap(index(i),index(j)) end do index(l+1)=index(j) index(j)=indext jstack=jstack+2 if (jstack > NSTACK) call nrerror('indexx: NSTACK too small') if (r-i+1 >= j-l) then istack(jstack)=r istack(jstack-1)=i r=j-1 else istack(jstack)=j-1 istack(jstack-1)=l l=i end if end if end do END SUBROUTINE indexx_i4b !BL SUBROUTINE icomp_xchg(iarr,i,j) USE nrtype; USE nrutil, ONLY : arth,assert_eq,nrerror,swap IMPLICIT NONE INTEGER(I4B), INTENT(IN) :: iarr(:) INTEGER(I4B), INTENT(INOUT) :: i,j INTEGER(I4B) :: swp if (iarr(j) < iarr(i)) then swp=i i=j j=swp end if END SUBROUTINE icomp_xchg END SUBROUTINE zroots end module mom_arma program princip use constants_module; use nrtype use mod_mat use mom_arma use random use stat_elementaires use momgarch use paraarmagarch implicit none integer(I4B) :: lag,n,lagmax,lagmax1,ordremax,ordre,ordre2,method integer(I4B) :: p0,pmax,i,aicbic real(SP) :: dum,meanx,sigx real(SP), allocatable :: x(:),var(:,:),varr(:,:),epsilest(:),epsnorm(:),rho(:),r(:) character (len=1), parameter :: virgul=',' open(8,file='identarmagarch.res',status='replace',position='rewind') open(6,file='identarmagarch.don',status='old',position='rewind') ! ! lecture de la serie ! write (8,*) " ___________________________________________" write (8,*) " Statistiques élementaires " write (8,*) " ___________________________________________" write (8,*) write (8,*) n=0 do read(6,*,end=2)dum n=n+1 end do 2 write(8,*)"nombre de données :",n rewind (6) allocate(x(n)) do i=1,n read(6,*)x(i) end do meanx=moy(x(1:n)) write (8,*) write(8,*)"moyenne des données :",meanx sigx=sig(x(1:n)) write (8,*) write(8,*)"écart-type des données :",sigx write(*,*)"nombre d'autocorrelations a prendre en compte ?" read(*,*)lag allocate(epsnorm(n),epsilest(n),var(lag,lag),rho(lag),varr(lag,lag),r(lag)) write(*,*)"methode utilisee pour estimer les covariances entre autocovariances empiriques ? ",& "(2=formule de bartlett generalisee, 3=methode tres robuste) " read(*,*)method pmax=20 ordremax=30 lagmax1=300 lagmax=min(lagmax1,int(n/2)) call estbartlettgene(x(1:n)-moy(x(1:n)),var(1:lag,1:lag),rho(1:lag),lagmax,ordremax,ordre,ordre2,epsilest(1:n),method) call autocorparemp(rho(1:lag),var(1:lag,1:lag),r(1:lag),varr(1:lag,1:lag)) if(method==2)then write (8,*) write (8,*) write (8,*) " ____________________________________________________________________________________________" write (8,*) " Estimation d'un intervalle de confiance à 95% pour les autocorrélations " write (8,*) " Formule de bartlett généralisée " write (8,*) " (pour modèle linéaire avec bruit différence de martingale symétrique) " write (8,*) " ____________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Retard borne inf rho borne sup" write(8,*)" " do i=1,lag write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,rho(i)-1.96*sqrt(var(i,i)/n),& rho(i),rho(i)+1.96*sqrt(var(i,i)/n) end do write (8,*) write (8,*) write (8,*) " _______________________________________________________________________________________________" write (8,*) " Estimation d'un intervalle de confiance à 95% pour les autocorrélations partielles " write (8,*) " _______________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Retard borne inf r borne sup" write(8,*)" " do i=1,lag write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,r(i)-1.96*sqrt(varr(i,i)/n),& r(i),r(i)+1.96*sqrt(varr(i,i)/n) end do write (8,*) write (8,*) write (8,*) " _______________________________________________________________________________________________" write (8,*) " Estimation des ordres ARMA-Méthode du coin " write (8,*) " _______________________________________________________________________________________________" write (8,*) call coin(rho(1:lag),var(1:lag,1:lag),n) end if if(method==3)then write(*,*)"Critere d'information utilise pour selectionner un AR ",& "qui ajuste la serie ? (1 pour AIC, 2 pour BIC) " call estsigmarho(x(1:n)-moy(x(1:n)),p0,pmax,var(1:lag,1:lag),rho(1:lag),aicbic) write (8,*) write (8,*) write (8,*) " ____________________________________________________________________________________________" write (8,*) " Estimation d'un intervalle de confiance à 95% pour les autocorrélations " write (8,*) " par utilisation de la densité spectrale d'une approximation AR " write (8,*) " (pour modéle linéaire faible) " write (8,*) " ____________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Retard borne inf rho borne sup" write(8,*)" " do i=1,lag write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,rho(i)-1.96*sqrt(var(i,i)/n),& rho(i),rho(i)+1.96*sqrt(var(i,i)/n) end do write (8,*) write (8,*) write (8,*) " _______________________________________________________________________________________________" write (8,*) " Estimation d'un intervalle de confiance à 95% pour les autocorrélations partielles " write (8,*) " _______________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Retard borne inf r borne sup" write(8,*)" " do i=1,lag write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,r(i)-1.96*sqrt(varr(i,i)/n),& r(i),r(i)+1.96*sqrt(varr(i,i)/n) end do write (8,*) write (8,*) write (8,*) " _______________________________________________________________________________________________" write (8,*) " Estimation des ordres ARMA-Méthode du coin " write (8,*) " _______________________________________________________________________________________________" write (8,*) call coin(rho(1:lag),var(1:lag,1:lag),n) end if write(*,*)"Critere d'information pour selectionner le modele AR multivarie (1 pour AIC, 2 pour BIC)?" read(*,*)aicbic epsnorm(ordre+1:n)=(epsilest(ordre+1:n)**2-moy(epsilest(ordre+1:n)**2))/sig(epsilest(ordre+1:n)**2) call estsigmarho(epsnorm(ordre+1:n),p0,pmax,var(1:lag,1:lag),rho(1:lag),aicbic) call autocorparemp(rho(1:lag),var(1:lag,1:lag),r(1:lag),varr(1:lag,1:lag)) write (8,*) write (8,*) write (8,*) " ____________________________________________________________________________________________" write (8,*) " Estimation d'un intervalle de confiance à 95% pour les autocorrélations " write (8,*) " des carrés des innovations " write (8,*) " ____________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Retard borne inf rho borne sup" write(8,*)" " do i=1,lag write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,rho(i)-1.96*sqrt(var(i,i)/n),& rho(i),rho(i)+1.96*sqrt(var(i,i)/n) end do write (8,*) write (8,*) write (8,*) " _______________________________________________________________________________________________" write (8,*) " Estimation d'un intervalle de confiance à 95% pour les autocorrélations partielles " write (8,*) " des carrés des innovations " write (8,*) " _______________________________________________________________________________________________" write (8,*) write (8,*) write(8,*)"Retard borne inf r borne sup" write(8,*)" " do i=1,lag write(8,fmt='(1x,i2,15x,f7.3,15x,f7.3,15x,f7.3)')i,r(i)-1.96*sqrt(varr(i,i)/n),& r(i),r(i)+1.96*sqrt(varr(i,i)/n) end do write (8,*) write (8,*) write (8,*) " _______________________________________________________________________________________________" write (8,*) " Estimation des ordres GARCH " write (8,*) " _______________________________________________________________________________________________" write (8,*) call coingarch(rho(1:lag),var(1:lag,1:lag),n) stop end program princip