!Wednesday, August 6, 2003 at 15:25 ! ! simulation d'un garch(p,q), calcule les autocorrelations empiriques, ! estime les bandes de confiance, calcule la stat portmanteau ! ! module constants_module !definition de quelques constantes utiles implicit none integer, public, parameter :: large_int_kind = selected_real_kind(4), & large_real_kind = selected_real_kind(8) real(kind=large_real_kind), public, parameter :: & pi = 3.14159265358979323846_large_real_kind, trois=3.0_large_real_kind, demi=.5_large_real_kind, & two_pi = 2*pi, un=1.0_large_real_kind, deux=2.0_large_real_kind, zero=0.0_large_real_kind end module constants_module 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 (*,*) '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 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 ! 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) ecrirt la matrice A use constants_module USE nrtype; USE nrutil, ONLY : assert_eq implicit none private private :: elmhes, hqr, vech, hh, kk, norm, normv,& invgauss, balanc, adamar, tred2,tqli, choldc public :: identity, inverse, vec, spectral_radius,kro, kromatvec,krovecmat,krovecvec,spectre,gauss,prod,resgauss,& diag,trace,jordan,sqrtmat,sqrtinv,invsim,sqrtmat2,sqrtinv2,writemat,determinant contains ! writemat(A) ecrirt la matrice A subroutine writemat(A) real (kind=large_real_kind), intent(in), dimension(:,:) :: A integer :: i, n n=size(A,1) do i=1,n write(88,fmt='(1x,7(f7.3,1x))') A(i,:) end do return end subroutine writemat function inverse(a) real(kind=large_real_kind), intent(in) :: a(:,:) real(kind=large_real_kind) :: 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 determinant(a) real(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind) :: determinant real(kind=large_real_kind) :: aout(size(a,1),size(a,1)) integer :: i determinant=un call ludcmp(a,aout) do i=1,size(a,1) determinant=determinant*aout(i,i) end do return contains subroutine ludcmp(ain,a) real (kind=large_real_kind), intent(in) :: ain(:,:) real (kind=large_real_kind), intent(out) :: a(:,:) real (kind=large_real_kind), parameter :: tiny=1.0e-20 integer, parameter :: nmax=100 real (kind=large_real_kind) ::indx(size(a,1)),vv(nmax) real (kind=large_real_kind) ::d,aamax,sum,dum integer :: 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(*,*) "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(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind) :: sqrtinv2(size(a,1),size(a,2)) real(kind=large_real_kind) :: p(size(a,1)), b(size(a,1),size(a,2)),sum integer :: 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(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind) :: sqrtmat2(size(a,1),size(a,2)) real(kind=large_real_kind) :: p(size(a,1)), b(size(a,1),size(a,2)) integer :: 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(kind=large_real_kind), intent(in out) :: a(:,:) real(kind=large_real_kind), intent(out) :: p(:) integer, intent(in) :: n integer :: i,j,k real(kind=large_real_kind) :: 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(*,*)"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(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind) :: invsim(size(a,1),size(a,2)) real(kind=large_real_kind) :: P(size(a,1),size(a,2)), Lambda(size(a,1),size(a,2)) integer :: 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(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind) :: sqrtinv(size(a,1),size(a,2)) real(kind=large_real_kind) :: P(size(a,1),size(a,2)), Lambda(size(a,1),size(a,2)) integer :: 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(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind) :: sqrtmat(size(a,1),size(a,2)) real(kind=large_real_kind) :: 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(kind=large_real_kind), intent(in) :: A(:,:) real(kind=large_real_kind), intent(out) :: P(:,:), Lambda(:,:) real(kind=large_real_kind) :: 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(kind=large_real_kind), intent(in out) :: a(:,:) real(kind=large_real_kind), intent(out) :: d(:), e(:) integer :: n,i,l,k,j real(kind=large_real_kind) :: 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 (kind=large_real_kind), intent(in out) :: d(:),e(:) real (kind=large_real_kind), intent(out) :: z(:,:) integer :: n,i,l,iter,m,k real(kind=large_real_kind) :: 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(*,*) '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 (kind=large_real_kind), intent(in), dimension(:) :: x, y real (kind=large_real_kind), dimension(size(x),size(y)) :: prod_result integer :: 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(kind=large_real_kind), intent(in), dimension(:,:) :: a real(kind=large_real_kind), dimension((size(a,1)*(size(a,1)+1))/2) :: vech_result integer :: 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(kind=large_real_kind), intent(in), dimension(:,:) :: a real(kind=large_real_kind), dimension(size(a,1)*size(a,2)) :: vec_result integer :: 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, intent(in) :: n real(kind=large_real_kind), dimension((n*(n+1))/2,n**2) :: kk_result real(kind=large_real_kind), dimension(n**2,(n*(n+1))/2) :: k integer :: i, j k=0.0_large_real_kind do j=1,n do i=j,n k((j-1)*n+i,((j-1)*(2*n-j))/2+i)=1.0_large_real_kind k((i-1)*n+j,((j-1)*(2*n-j))/2+i)=1.0_large_real_kind 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, intent(in) :: n real(kind=large_real_kind), dimension((n*(n+1))/2,n**2) :: hh_result integer :: i, j, id1, id2 hh_result=0.0_large_real_kind 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_large_real_kind 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, intent(in) :: n, m real (kind=large_real_kind), intent(in), dimension(:,:) :: aa, bb real (kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: invgauss_result ! local variables real (kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: a real (kind=large_real_kind), dimension(size(bb,1),size(bb,2)) :: b real (kind=large_real_kind) :: big, dum, pivinv integer, dimension(n) :: ipiv, indxr, indxc integer :: i, j, k, l, ll, irow, icol a=aa b=bb ipiv=0 do i=1,n big=0.0_large_real_kind 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_large_real_kind) then write(unit=4,fmt=*)"matrice singuliere" return end if pivinv=1.0_large_real_kind/a(icol,icol) a(icol,icol)=1.0_large_real_kind 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_large_real_kind 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a real (kind=large_real_kind), dimension(size(a,1),size(a,1)) :: inv_result real (kind=large_real_kind), dimension(size(a,1),1) :: zero zero=0.0_large_real_kind inv_result=invgauss(a,zero,size(a,1),1) return end function inv function balanc(aa,n) result (balanc_result) real(kind=large_real_kind), intent(in), dimension(:,:) :: aa integer, intent(in) :: n real(kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: balanc_result real(kind=large_real_kind), parameter :: rradix=2.0_large_real_kind real(kind=large_real_kind), parameter :: sqrdx=rradix**2 real(kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: a integer :: i, j, last real(kind=large_real_kind) :: c, r, g, f, s a=aa boucle1 : do last=1 do i=1,n c=0.0_large_real_kind r=0.0_large_real_kind 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_large_real_kind .and. r /= 0.0_large_real_kind)then g=r/rradix f=1.0_large_real_kind 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_large_real_kind*s) then last=0 g=1.0_large_real_kind/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 (kind=large_real_kind), intent(in), dimension(:,:) :: aa integer, intent(in) :: n real (kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: elmhes_result real (kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: a integer :: m, i, j real (kind=large_real_kind) :: x, y a=aa a=balanc(a,n) do m=2,n-1 x=0.0_large_real_kind 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_large_real_kind) then do i=m+1,n y=a(i,m-1) if (y /= 0.00_large_real_kind) 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(kind=large_real_kind), intent(in), dimension(:,:) :: aa integer, intent(in) :: n complex(kind=large_real_kind), dimension(size(aa,1)) :: hqr_result real(kind=large_real_kind), dimension(size(aa,1),size(aa,2)) :: a real(kind=large_real_kind), dimension(size(aa,1)) :: wr, wi integer :: m, i, j, k, nn, l, its,ii real(kind=large_real_kind) :: 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_large_real_kind 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_large_real_kind)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_large_real_kind 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_large_real_kind*(y-x) q=p**2+w z=sqrt(abs(q)) x=x+t if (q >= 0.0_large_real_kind) then z=p+sign(z,p) wr(nn)=x+z wr(nn-1)=wr(nn) if (z /= 0.0_large_real_kind) then wr(nn)=x-w/z end if wi(nn)=0.0_large_real_kind wi(nn-1)=0.0_large_real_kind 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(:),large_real_kind) 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_large_real_kind*s y=x w=-0.43750_large_real_kind*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_large_real_kind if (i /= m+2)then a(i,i-3)=0.0_large_real_kind 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_large_real_kind if (k /= nn-1) then r=a(k+2,k-1) end if x=abs(p)+abs(q)+abs(r) if (x /= 0.0_large_real_kind) 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_large_real_kind) 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(:),large_real_kind) 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a real (kind=large_real_kind) :: spectral_result real (kind=large_real_kind), dimension(size(a,1),size(a,2)) :: b integer :: 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a complex(kind=large_real_kind), dimension(size(a,1)) :: spectre_result real (kind=large_real_kind), dimension(size(a,1),size(a,1)) :: b integer :: 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 (kind=large_real_kind), intent(in), dimension(:) :: x, y real (kind=large_real_kind), dimension(size(x)) :: adamar_result integer :: 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, intent(in) :: n, m real (kind=large_real_kind), intent(in), dimension(:,:) :: aa, bb real (kind=large_real_kind),dimension(n,m) :: resgauss_result ! local variables real (kind=large_real_kind) :: big, dum, pivinv integer, dimension(n) :: ipiv, indxr, indxc real (kind=large_real_kind), dimension(n,n) :: a real (kind=large_real_kind), dimension(n,m) :: b integer :: i, j, k, l, ll, irow, icol a=aa b=bb ipiv=0 do i=1,n big=0.0_large_real_kind 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_large_real_kind) then resgauss_result=b write(unit=4,fmt=*)"matrice singuliere" return end if pivinv=1.0_large_real_kind/a(icol,icol) a(icol,icol)=1.0_large_real_kind 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_large_real_kind 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, intent(in) :: n, m real (kind=large_real_kind), intent(in out), dimension(:,:) :: a, b ! local variables real (kind=large_real_kind) :: big, dum, pivinv integer, dimension(n) :: ipiv, indxr, indxc integer :: i, j, k, l, ll, irow, icol ipiv=0 do i=1,n big=0.0_large_real_kind 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_large_real_kind) then write(unit=4,fmt=*)"matrice singuliere" return end if pivinv=1.0_large_real_kind/a(icol,icol) a(icol,icol)=1.0_large_real_kind 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_large_real_kind 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 (kind=large_real_kind), intent(in), dimension(:) :: a real (kind=large_real_kind), intent(in), dimension(:) :: b real (kind=large_real_kind), dimension(size(a,1)*size(b)) :: kro_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a real (kind=large_real_kind), intent(in), dimension(:) :: b real (kind=large_real_kind), dimension(size(a,1)*size(b),size(a,2)) :: kro_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a real (kind=large_real_kind), intent(in), dimension(:) :: b real (kind=large_real_kind), dimension(size(a,1)*size(b),size(a,2)) :: kro_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a, b real (kind=large_real_kind), dimension(size(a,1)*size(b,1),size(a,2)*size(b,2)) :: kro_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:,:) :: a real (kind=large_real_kind) :: trace_result integer :: 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_large_real_kind 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 (kind=large_real_kind), intent(in), dimension(:) :: vec real (kind=large_real_kind), dimension(size(vec),size(vec)) :: diag_result integer :: i,n n=size(vec) diag_result=0.0_large_real_kind 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, intent(in) :: ia real (kind=large_real_kind), dimension(ia,ia) :: identity_result integer :: i identity_result=0.0_large_real_kind do i=1,ia identity_result(i,i)=1.0_large_real_kind end do return end function identity function norm(a) result(norm_result) !!!! ! norm=somme des valeurs absolues des elements de la matrice a !!!! real (kind=large_real_kind), intent(in), dimension(:,:) :: a real (kind=large_real_kind) :: 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 (kind=large_real_kind), intent(in), dimension(:) :: a real (kind=large_real_kind) :: normv_result normv_result=sum(abs(a)) return end function normv end module mod_mat module parabruit ! parametres du sous-programme faible quand m=7 ou m=20 use constants_module implicit none real(kind=large_real_kind), public, save :: cc0,pp0,omega00,alpha00 end module parabruit 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 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 (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: 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, intent(in) :: nu real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: 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(kind=large_real_kind),intent(in) :: a,x real(kind=large_real_kind) :: gammp,gammcf,gln if(x.lt.zero.or.a.le.zero)then write(*,*)"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(kind=large_real_kind),intent(in) :: a,x real(kind=large_real_kind),intent(out) :: gamser,gln real(kind=large_real_kind), parameter :: eps=0.0000003 integer, parameter :: itmax=100 real(kind=large_real_kind) :: ap,sum,del integer :: n gln=gammln(a) if(x.le.zero)then if(x.lt.zero)write(*,*)"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(*,*)"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(kind=large_real_kind),intent(in) :: a,x real(kind=large_real_kind),intent(out) :: gammcf,gln real(kind=large_real_kind), parameter :: eps=0.0000003 integer, parameter :: itmax=100 real(kind=large_real_kind) :: gold,a0,a1,b0,b1,fac,an,ana,anf,g integer :: 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(*,*)"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(kind=large_real_kind) :: cof(6),stp,half,one,fpf,x,tmp,ser real (kind=large_real_kind), intent(in) :: xx real(kind=large_real_kind) :: gammln integer :: j cof=(/76.18009173_large_real_kind,-86.50532033_large_real_kind,24.01409822_large_real_kind,& -1.231739516_large_real_kind,.00120858003_large_real_kind,-.00000536382_large_real_kind/) stp=2.50662827465_large_real_kind half=0.5_large_real_kind one=1.0_large_real_kind fpf=5.5_large_real_kind 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, intent(in) :: m real (kind=large_real_kind), intent(in) :: x real(kind=large_real_kind) :: 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(kind=large_real_kind) :: cof(6),stp,half,one,fpf,x,tmp,ser real (kind=large_real_kind), intent(in) :: xx real(kind=large_real_kind) :: gamma integer :: j cof=(/76.18009173_large_real_kind,-86.50532033_large_real_kind,24.01409822_large_real_kind,& -1.231739516_large_real_kind,.00120858003_large_real_kind,-.00000536382_large_real_kind/) stp=2.50662827465_large_real_kind half=0.5_large_real_kind one=1.0_large_real_kind fpf=5.5_large_real_kind 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, intent(in) :: nu real (kind=large_real_kind), intent(in) :: p real (kind=large_real_kind) :: invfdr2 real (kind=large_real_kind) :: xmin, xmax, middle real (kind=large_real_kind), parameter :: tol=.0001 interface function f(x,n) use constants_module implicit none integer, intent(in) :: n real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: 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, intent(in) :: nu real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: Khi2b real (kind=large_real_kind) :: dum if (nu < 1 ) then write(*,*)"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, intent(in) :: nu real (kind=large_real_kind), intent(in) :: prob real (kind=large_real_kind) :: 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 (kind=large_real_kind), intent(in) :: p real (kind=large_real_kind) :: invfdr real (kind=large_real_kind) :: xmin, xmax, middle real (kind=large_real_kind), parameter :: tol=.0001 interface function f(x) use constants_module implicit none real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: 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 (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: Phi real (kind=large_real_kind), parameter :: b1=0.319381530, & b2=-0.356563782, & b3=1.781477937, & b4=-1.821255978, & b5=1.330274429 real (kind=large_real_kind) :: 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(kind=large_real_kind) :: Phiinv real(kind=large_real_kind), intent(in) :: prob Phiinv = invfdr(Phi,prob) return end function Phiinv end module fdr module parafaint ! parametres de la fonction a integrer use constants_module implicit none integer, public, save :: m0 real(kind=large_real_kind), public, save :: a00, q0 end module parafaint module integral_mod use constants_module use fdr implicit none private public :: qromb, midexp, midpnt, trapzd, integral, integral2, integral3, omega2 contains subroutine qromb(func,a,b,ss) real (kind=large_real_kind), intent (in) :: a,b real (kind=large_real_kind), intent (out) :: ss real (kind=large_real_kind), parameter :: eps=0.000001_large_real_kind integer, parameter :: jmax=10,jmaxp=jmax+1,k=5,km=4 integer :: j,l,it real (kind=large_real_kind) :: s(jmaxp),h(jmaxp),dss interface function func(x) use constants_module implicit none real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: func end function func end interface h(1)=un do j=1,jmax call trapzd(func,a,b,s(j),j,it) if (j.ge.k) then l=j-km call polint(h(l:l+k),s(l:l+k),zero,ss,dss) if (abs(dss).lt.eps*abs(ss)) return end if s(j+1)=s(j) h(j+1)=0.25_large_real_kind*h(j) end do write(*,*)"too many steps." contains subroutine polint(xa,ya,x,y,dy) real (kind=large_real_kind), intent (in) :: xa(:),ya(:),x real (kind=large_real_kind), intent (out) :: y,dy integer :: n,ns,i,m real (kind=large_real_kind) :: c(size(xa)),d(size(xa)),dif,dift,ho,hp,den,w n=size(xa) ns=1 dif=abs(x-xa(1)) do i=1,n dift=abs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift end if c(i)=ya(i) d(i)=ya(i) end do y=ya(ns) ns=ns-1 do m=1,n-1 do i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if(den.eq.zero)write(*,*)"pb dans polint" den=w/den d(i)=hp*den c(i)=ho*den end do if (2*ns.lt.n-m)then dy=c(ns+1) else dy=d(ns) ns=ns-1 end if y=y+dy end do return end subroutine polint end subroutine qromb subroutine midexp(func,aa,s,n,it) real (kind=large_real_kind), intent (in) :: aa real (kind=large_real_kind), intent (in out) :: s integer, intent(in) :: n integer, intent(in out ) :: it real (kind=large_real_kind) :: ddel, del, sum, x,a,b integer :: tnm,j interface function func(x) use constants_module implicit none real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: func end function func end interface b=exp(-aa) a=zero if (n.eq.1) then s=(b-a)*func(-log(demi*(a+b)))/(demi*(a+b)) it=1 else tnm=it del=(b-a)/(trois*tnm) ddel=del+del x=a+demi*del sum=zero do j=1,it sum=sum+func(-log(x))/x x=x+ddel sum=sum+func(-log(x))/x x=x+del end do s=(s+(b-a)*sum/tnm)/trois it=3*it end if return end subroutine midexp subroutine midpnt(func,a,b,s,n,it) real (kind=large_real_kind), intent (in) :: a, b real (kind=large_real_kind), intent (in out) :: s integer, intent(in) :: n integer, intent(in out ) :: it real (kind=large_real_kind) :: ddel, del, sum, x integer :: tnm,j interface function func(x) use constants_module implicit none real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: func end function func end interface if (n.eq.1) then s=(b-a)*func(demi*(a+b)) it=1 else tnm=it del=(b-a)/(trois*tnm) ddel=del+del x=a+demi*del sum=zero do j=1,it sum=sum+func(x) x=x+ddel sum=sum+func(x) x=x+del end do s=(s+(b-a)*sum/tnm)/trois it=3*it end if return end subroutine midpnt ! trapzd(func,a,b,s,n,it): s est l'integrale de func entre a et b subroutine trapzd(func,a,b,s,n,it) real (kind=large_real_kind), intent (in) :: a, b real (kind=large_real_kind), intent (in out) :: s integer, intent(in) :: n integer, intent(in out ) :: it real (kind=large_real_kind) :: del, sum, x integer :: tnm,j interface function func(x) use constants_module implicit none real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: func end function func end interface if (n.eq.1) then s=demi*(b-a)*(func(a)+func(b)) it=1 else tnm=it del=(b-a)/tnm x=a+demi*del sum=zero do j=1,it sum=sum+func(x) x=x+del end do s=demi*(s+(b-a)*sum/tnm) it=2*it end if return end subroutine trapzd recursive function integral3 (a, b, epsil,m,q,a0) & result (integral_result) implicit none real (kind=large_real_kind), intent (in) :: a, b, epsil,q,a0 integer, intent(in) :: m real (kind=large_real_kind) :: integral_result real (kind=large_real_kind) :: h, mid real (kind=large_real_kind) :: one_trapezoid_area, two_trapezoid_area real (kind=large_real_kind) :: left_area, right_area h = b - a mid = (a + b) /2 one_trapezoid_area = h * (Khi2((a-q)/(a0**(2*m)),1)*fKhi2(a,m-1) + & Khi2((b-q)/(a0**(2*m)),1)*fKhi2(b,m-1))/2 two_trapezoid_area = h/2 * (Khi2((a-q)/(a0**(2*m)),1)*fKhi2(a,m-1) + & Khi2((mid-q)/(a0**(2*m)),1)*fKhi2(mid,m-1))/2 + & h/2 * (Khi2((mid-q)/(a0**(2*m)),1)*fKhi2(mid,m-1) + Khi2((b-q)/(a0**(2*m)),1)*fKhi2(b,m-1))/2 if (abs(one_trapezoid_area - two_trapezoid_area) & < trois * epsil) then integral_result = two_trapezoid_area else left_area = integral3(a, mid, epsil/2, m,q,a0) right_area = integral3(mid, b, epsil/2, m,q,a0) integral_result = left_area + right_area end if return end function integral3 recursive function integral2 (a, b, epsil,m) & result (integral_result) implicit none real (kind=large_real_kind), intent (in) :: a, b, epsil integer, intent(in) :: m real (kind=large_real_kind) :: integral_result real (kind=large_real_kind) :: h, mid real (kind=large_real_kind) :: one_trapezoid_area, two_trapezoid_area real (kind=large_real_kind) :: left_area, right_area h = b - a mid = (a + b) /2 one_trapezoid_area = h * (fKhi2(a,m) + fKhi2(b,m)) / deux two_trapezoid_area = h/2 * (fKhi2(a,m) + fKhi2(mid,m)) / deux + & h/2 * (fKhi2(mid,m) + fKhi2(b,m)) / deux if (abs(one_trapezoid_area - two_trapezoid_area) & < trois * epsil) then integral_result = two_trapezoid_area else left_area = integral2(a, mid, epsil/2, m) right_area = integral2(mid, b, epsil/2, m) integral_result = left_area + right_area end if return end function integral2 recursive function integral (a, b, epsil) & result (integral_result) implicit none real (kind=large_real_kind), intent (in) :: a, b, epsil real (kind=large_real_kind) :: integral_result real (kind=large_real_kind) :: h, mid real (kind=large_real_kind) :: one_trapezoid_area, two_trapezoid_area real (kind=large_real_kind) :: left_area, right_area h = b - a mid = (a + b) /2 one_trapezoid_area = h * (omega2(a) + omega2(b)) / deux two_trapezoid_area = h/2 * (omega2(a) + omega2(mid)) / deux + & h/2 * (omega2(mid) + omega2(b)) / deux if (abs(one_trapezoid_area - two_trapezoid_area) & < trois * epsil) then integral_result = two_trapezoid_area else left_area = integral ( a, mid, epsil / 2) right_area = integral ( mid, b, epsil / 2) integral_result = left_area + right_area end if return end function integral function omega2(x) implicit none real (kind=large_real_kind) :: omega2 real (kind=large_real_kind), intent(in) :: x omega2=zero if (abs(x) .lt. un) omega2=un-abs(x) ! if( abs(x) .lt. un) omega2=un ! omega2=omega2**2 return end function omega2 end module integral_mod module stat_elementaires use constants_module use mod_mat implicit none private public :: tri, moy, sig, med, quart, autocov, autocor, autocov0, autocor0, autocovmult !tri(x) est le fichier x trie par ordre croissant contains function autocovmult(x,k) result(autocov_result) real (kind=large_real_kind), intent(in), dimension(:,:) :: x integer, intent(in) :: k real (kind=large_real_kind) :: autocov_result(0:k,1:size(x,2),1:size(x,2)),m(1:size(x,2)) integer :: 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 (kind=large_real_kind), intent(in), dimension(:) :: x integer, intent(in) :: k real (kind=large_real_kind) :: autocov_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:) :: x integer, intent(in) :: k real (kind=large_real_kind) :: autocor_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:) :: x integer, intent(in) :: k real (kind=large_real_kind) :: autocov_result integer :: 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 (kind=large_real_kind), intent(in), dimension(:) :: x integer, intent(in) :: k real (kind=large_real_kind) :: autocor_result integer :: 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 (kind=large_real_kind), intent(in),dimension(:) :: x real (kind=large_real_kind),dimension(size(x)) :: tri_result integer :: n, l, ir, i, j real (kind=large_real_kind) :: 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 (kind=large_real_kind), intent(in),dimension(:) :: x real (kind=large_real_kind) :: med_result integer :: n,id real (kind=large_real_kind) :: rest real (kind=large_real_kind),dimension(size(x)) :: xx n=size(x) xx=tri(x) id=int(1.0_large_real_kind*n/2) rest=1.0_large_real_kind*n/2 -id if (rest == 0.00_large_real_kind) 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 (kind=large_real_kind), intent(in),dimension(:) :: x real (kind=large_real_kind), intent(out) :: mmin, q1, q2, q3, mmax integer :: n real (kind=large_real_kind) :: rest real (kind=large_real_kind), dimension(size(x)) :: xx n=size(x) xx=tri(x) mmin=xx(1) mmax=xx(n) q2=med(xx) rest=1.0_large_real_kind*n/2 -int(1.0_large_real_kind*n/2) if (rest == 0.00_large_real_kind) then q1=med(xx(1:int(1.0_large_real_kind*n/2))) q3=med(xx(int(1.0_large_real_kind*n/2)+1:n)) else q1=med(xx(1:int(1.0_large_real_kind*n/2))) q3=med(xx(int(1.0_large_real_kind*n/2)+2:n)) end if return end subroutine quart function moy(x) result(moy_result) real (kind=large_real_kind), intent(in), dimension (:) :: x real (kind=large_real_kind) :: moy_result integer :: n n=size(x) moy_result=sum(x(1:n))/n return end function moy function sig(x) result(sig_result) real (kind=large_real_kind), intent(in),dimension(:) :: x real (kind=large_real_kind) :: sig_result integer :: 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 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 stat_elementaires use parabruit use fdr use mod_mat implicit none private public :: random_normal, random_exponential, random_t, chaine_Markov,ar,simarma,& faible, garch, estimcov, estimcoviid, LM_ARCH contains subroutine simarma(a,b,p,q,sigma,n,xmoy,epsilon,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 ! 50 valeur initiale sont eliminees ! epsilon est un bruit non correle centre de variance 1 implicit none REAL(kind=large_real_kind), intent(in) :: a(:), b(:), xmoy, sigma, epsilon(-49:) REAL(kind=large_real_kind), intent(out) :: x(:) integer, intent(in) :: p, q, n ! Local variables REAL(kind=large_real_kind) :: xx(-49:n) integer :: i,j do i=-49,-49+p+q xx(i)=sigma*epsilon(i) end do do i=-49+p+q+1,n xx(i)=sigma*epsilon(i) do j=1,p xx(i)=xx(i)+a(j)*xx(i-j) end do do j=1,q xx(i)=xx(i)+b(j)*sigma*epsilon(i-j) end do end do do i=1,n x(i)=xx(i)+xmoy end do return end subroutine simarma subroutine faible(epsilon,m,n) ! si m=0 ! epsilon(t)=eta(t) ! avec eta(t) suite i.i.d. de loi normale standard implicit none integer, intent(in) :: n,m REAL(kind=large_real_kind), intent(out) :: epsilon(:) ! Local variables integer :: i epsilon=zero select case (m) case (0) do i=1,n epsilon(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(kind=large_real_kind), intent(in) :: omega,alpha(:),beta(:) integer, intent(in) :: n,m REAL(kind=large_real_kind) :: garch(1:n) ! Local variables integer, parameter :: nvalinit=5000 REAL(kind=large_real_kind) :: eta(1:n+nvalinit),h(1:n+nvalinit), x(1:n+nvalinit) integer :: 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 LM_ARCH(x,lag,lm,pval) implicit none real(kind=large_real_kind), intent(in) :: x(:) integer, intent(in) :: lag real(kind=large_real_kind), intent(out) :: lm,pval real(kind=large_real_kind) :: Xmat(size(x),lag+1),XtX(lag+1,lag+1),Yvec(size(x)),invXtX(lag+1,lag+1),& dum1(lag+1),dum2(lag+1),num,den real(kind=large_real_kind), parameter :: tol=0.00000001 ! local variables integer :: n,i,j n=size(x) Xmat(1:n,1)=un Yvec(1:n)=x(1:n)**2-moy(x(1:n)**2) ! Yvec(1:n)=x(1:n)**2 do i=n,1,-1 do j=1,lag if((i-j.gt.0).and.((i-j).le.n))then Xmat(i,j+1)=x(i-j)**2 else Xmat(i,j+1)=zero end if end do end do XtX=matmul(transpose(Xmat),Xmat) if(abs(determinant(XtX)).lt.tol)then write(4,*)"matrice non inversible dans LM_ARCH" return end if invXtX=inverse(XtX) dum1=matmul(Yvec,Xmat) dum2=matmul(invXtX,dum1) num=dot_product(dum1,dum2) den=dot_product(Yvec,Yvec) if(abs(den).lt.tol)then write(4,*)"statistique non calculabledans LM_ARCH" return end if lm=n*num/den pval=un-khi2(lm,lag) return end subroutine LM_ARCH subroutine estimcoviid(x,lag,rhom,sigmam,qm,pval) implicit none real(kind=large_real_kind), intent(in) :: x(:) integer, intent(in) :: lag real(kind=large_real_kind), intent(out) :: rhom(:),sigmam(:),qm,pval ! local variables integer :: 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) end do pval=un-khi2(qm,lag) return end subroutine estimcoviid subroutine estimcov(x,lag,rhom,sigmam,qm,pval) implicit none real(kind=large_real_kind), intent(in) :: x(:) integer, intent(in) :: lag real(kind=large_real_kind), intent(out) :: rhom(:),sigmam(:),qm,pval ! local variables real(kind=large_real_kind) :: gamma0 integer :: 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(kind=large_real_kind), intent(in) :: a(:) integer, intent(in) :: n REAL(kind=large_real_kind) :: ar(1:n) ! Local variables integer, parameter :: nvalinit=500 REAL(kind=large_real_kind) :: eta(1:n+nvalinit), x(1:n+nvalinit) integer :: 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(kind=large_real_kind) :: random_normal ! Local variables REAL(kind=large_real_kind) :: s = 0.449871, t = -0.386595, a = 0.19600, b = 0.25472, & r1 = 0.27597, r2 = 0.27846, u, v, x, y, q REAL(kind=large_real_kind), 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 (kind=large_real_kind) :: random_exponential ! Local variable REAL (kind=large_real_kind) :: r REAL (kind=large_real_kind), 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(kind=large_real_kind) :: random_t INTEGER, INTENT(IN) :: m ! Local variables REAL (kind=large_real_kind) :: s, c, a, f, g, r, x, v REAL (kind=large_real_kind), 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 :: mm = 0 mm=0 !SAVE s, c, a, f, g IF (m < 1) THEN write(*,*)"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(kind=large_real_kind), intent(in) :: p(:,:), p0(:) integer, intent(out) :: delta(-49:) integer, intent(in) :: d, n ! Local variables real(kind=large_real_kind) :: u, fdr integer :: 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 formequad use stat_elementaires use constants_module use fdr use parafaint use integral_mod implicit none private private :: fquad public :: fquadsimple,niveau,faint,pvalue contains function pvalue(m,a,q) real(kind=large_real_kind), intent(in) :: a,q integer, intent(in) :: m integer :: mult(2) real(kind=large_real_kind) :: pvalue,d(2) real(kind=large_real_kind),parameter :: tol=0.1 if(m<1)then write(*,*)"m trop petit" return end if if(a**(2*m) < tol)then pvalue=un-khi2(q,m-1) else d(1)=a**(2*m) d(2)=un mult(1)=1 mult(2)=m-1 pvalue=un-fquadsimple(d(1:2),mult(1:2),q) end if return end function pvalue function faint(x) ! fonction a integrer implicit none real (kind=large_real_kind), intent(in) :: x real (kind=large_real_kind) :: faint faint=(un-Khi2(q0-(a00**(2*m0))*x,m0-1))*fkhi2(x,1) return end function faint function niveau(m,alpha) real(kind=large_real_kind), intent(in) :: alpha integer, intent(in) :: m integer :: it,n real(kind=large_real_kind) :: niveau,q,s,s0 real(kind=large_real_kind),parameter :: tol=0.0001 if(m<1)then write(*,*)"m trop petit" return end if q=Khi2inv(un-alpha,m-1) write(8,*)"m-1,un-alpha,q",m-1,un-alpha,q q0=q m0=m it=0 s=zero do n=1,30 s0=s call trapzd(faint,zero,q,s,n,it) if(abs(s-s0)10)exit write(8,*)"n,s",n,s end do niveau=alpha+s return end function niveau function fquadsimple(d,mult,x) real(kind=large_real_kind), intent(in) :: d(:),x integer, intent(in) :: mult(:) real(kind=large_real_kind), parameter :: eps1=.0001_large_real_kind,eps2=.0001_large_real_kind real(kind=large_real_kind), parameter :: tol=.01_large_real_kind real(kind=large_real_kind) :: fquadsimple,y,moy,var real(kind=large_real_kind) :: delta(size(d)),dtri(size(d)) integer :: n,k,i n=size(d) k=n dtri=tri(d) if(minval(d)moy+5*sqrt(var))then fquadsimple=un return end if if(x=2)then call fquad(dtri(1:k),mult(1:k),fquadsimple,eps1,eps2,delta(1:k),x) else y=x/dtri(1) fquadsimple=khi2(y,1) end if if (fquadsimple < zero )then write(8,*) "FAILURE OF THE HIMHOF ALGORITM: fquad <0" write(8,*) "x=",x write(8,*) "vp=",dtri(1:k) end if if (fquadsimple > un)then write(8,*) "FAILURE OF THE HIMHOF ALGORITM: fquad >0" write(8,*) "x=",x write(8,*) "vp=",dtri(1:k) end if if (fquadsimple < zero)fquadsimple=zero if (fquadsimple > un)fquadsimple=un return end function fquadsimple subroutine fquad(d,mult,fd,eps1,eps2,delta,xkrit) ! This subroutine comes from the chapter 9 of: ! 'ON THE THEORY AND APPLICATION OF THE GENERAL LINEAR MODEL' DE ! J. KOERTS ET A.P.J. ABRAHAMSE, ROTTERDAM UNIVERSITY PRESS, 1969. ! PURPOSE : FIND THE DISTRIBUTION FUNCTION F(X) = P(Q<=X) OF ! QUADRATIC FORMS IN NORMAL VARIABLES ! Note: We suggest EPS1 = EPS2 = .0001 ! DESCRIPTION OF PARAMETERS ! D = A VECTOR OF LENGTH NR CONTAINING THE NON-ZERO LATENT ROOTS ! NR = NUMBER OF NON-ZERO LATENT ROOTS ! MULT = A VECTOR OF LENGTH NR CONTAINING THE ORDERS OF MULTIPLICITY ! OF THE CORRESPONDING ROOTS ! FD = VALUE OF THE DISTRIBUTION FUNCTION ! EPS1 = DESIRED TRUNCATION ERROR TU ! EPS2 = DESIRED ACCURACY OF NUMERICAL INTEGRATION ! DELTA = A VECTOR OF LENGTH NR CONTAINING THE NON-CENTRALITY PARAMETERS ! XKRIT = VALUE OF X ! METHOD : THE INTEGRAL DERIVED BY IMHOF IS USED AND NUMERICALLY ! INTEGRATED BY SIMPSON'S RULE ! REFERENCE : 'COMPUTING THE DISTRIBUTION OF QUADRATIC FORMS IN ! NORMAL VARIABLES' BY J.P. IMHOF ! BIOMETRIKA (1961), 48, 3 AND 4, P. 419 real(kind=large_real_kind), intent(in) :: D(:),DELTA(:),EPS1,EPS2,XKRIT integer, intent(in) :: MULT(:) real(kind=large_real_kind), intent(out) :: FD real(kind=large_real_kind)::xk,slam,sdelt,ub,sum,dd,u,tu,teta,rho,fbl,sumk,range,vint2,c1,c2,fd1,step integer ::i,nr,kskip,mh,nit,k nr=size(d) XK=zero SLAM=zero SDELT=zero DO I=1,NR XK=XK+demi*MULT(I) SLAM=SLAM+demi*MULT(I)*ALOG(ABS(D(I))) SDELT=SDELT+demi*DELTA(I) end do ! COMPUTATION OF THE UPPER-BOUND OF THE INTEGRAL GIVEN THE ! TRUNCATION ERROR UB=EXP(-(ALOG(EPS1*XK)+1.14472989_large_real_kind+(SLAM+SDELT))/XK) 7 SUM=zero DO I=1,NR DD=(D(I)*UB)**2 SUM=SUM+DELTA(I)*DD/(un+DD) end do SUM=demi*SUM TU=EXP(1.14472989_large_real_kind+ALOG(XK)+(SLAM+SUM)+XK*ALOG(UB)) IF(un/TU-EPS1) 20,9,9 9 UB=UB+5*un/XK go to 7 ! COMPUTATION OF THE INTEGRAND GIVEN THE VALUE OF U 10 IF(U) 15,15,11 11 TETA=zero SUM=zero RHO=un DO I=1,NR C1=D(I)*U C2=un+C1**2 TETA=TETA+MULT(I)*ATAN(C1)+DELTA(I)*C1/C2 SUM=SUM+DELTA(I)*C1**2/C2 RHO=RHO*(C2**(0.25_large_real_kind*MULT(I))) end do TETA=demi*(TETA-XKRIT*U) RHO=RHO*EXP(demi*SUM) FBL=SIN(TETA)/(RHO*U) go to 18 15 FBL=zero DO I=1,NR 16 FBL=FBL+D(I)*(MULT(I)+DELTA(I)) end do FBL=demi*(FBL-XKRIT) 18 go to (21,22,23,24,25), KSKIP ! EVALUATION OF THE INTEGRAL BY SIMPSON'S RULE 20 RANGE=UB MH=1 U=RANGE*demi KSKIP=1 go to 10 21 SUMK=FBL*RANGE*deux/(3*un) U=zero KSKIP=2 go to 10 22 VINT2=SUMK+FBL*RANGE/(6*un) U=UB KSKIP=3 go to 10 23 VINT2=VINT2+FBL*RANGE/(6*un) FD=demi-.318309886_large_real_kind*VINT2 grandeboucle : DO NIT=1,14 FD1=FD VINT2=(VINT2-SUMK*demi)*demi MH=2*MH STEP=RANGE/MH U=STEP*demi KSKIP=4 go to 10 24 SUMK=FBL ! Modification of the original source for best results on ! many platforms... ! DO K=2,MH k=1 2525 k=k+1 U=U+STEP KSKIP=5 go to 10 25 SUMK=SUMK+FBL ! end do ! Modification of the original source for best results on ! many platforms... if(k.lt.mh) go to 2525 SUMK=SUMK*STEP*2./3. VINT2=VINT2+SUMK FD=demi-.318309886*VINT2 IF( NIT >3 .and. (ABS(FD1-FD) < EPS2)) exit grandeboucle end do grandeboucle 29 RETURN END SUBROUTINE FQUAD end module formequad module portmanteau ! boxpierce(x,m)=n*(\rho^2(1)+...+\rho^2(m)) use constants_module use stat_elementaires use mod_mat implicit none private public :: estsigmarhostrong,Ljungbox,boxpierce,estar1,estarmult,estarmultmod,& estR,estsigmarho,boxpiercemod,determinep contains subroutine estsigmarhostrong(epsilon,a,m,sigmarho) ! estime \Omega_{\rho} pour un AR(1) fort real(kind=large_real_kind), intent(in) :: epsilon(:),a integer, intent(in) :: m real(kind=large_real_kind), intent(out) :: sigmarho(:,:) real(kind=large_real_kind) :: Lambda(m) integer :: i do i=1,m Lambda(i)=-a**(i-1) end do sigmarho(1:m,1:m)=identity(m)-(un-a**2)*prod(Lambda(1:m),Lambda(1:m)) return end subroutine estsigmarhostrong function Ljungbox(x,m) ! Ljung-Box statistic based on the first $m$ autocorrelations real(kind=large_real_kind) :: Ljungbox real(kind=large_real_kind), intent(in) :: x(:) real(kind=large_real_kind) :: s0 integer, intent(in) :: m integer :: i,n n=size(x) s0=sum(x(1:n)**2) Ljungbox=zero do i=1,m Ljungbox=Ljungbox+((n*(n+2)*un)/(n-i))*(sum(x(1:n-i)*x(1+i:n))/s0)**2 end do return end function Ljungbox function boxpiercemod(epsilon,m,Omega) ! modified Box-Pierce statistic based on the first $m$ residuals autocorrelations real(kind=large_real_kind) :: boxpiercemod real(kind=large_real_kind), intent(in) :: epsilon(:),Omega(:,:) real(kind=large_real_kind) :: s0,rho(m) integer, intent(in) :: m integer :: i,n n=size(epsilon) s0=sum(epsilon(1:n)**2)/sqrt(un*n) do i=1,m rho(i)=sum(epsilon(1:n-i)*epsilon(1+i:n))/s0 end do boxpiercemod=dot_product(rho(1:m),matmul(inverse(Omega(1:m,1:m)),rho(1:m))) return end function boxpiercemod subroutine estsigmarho(epsilon,a,p,m,mmax,pmax,sigmarho) ! estime \Omega_{\rho} pour un AR(1) faible par regressions auxiliaires AR(p) real(kind=large_real_kind), intent(in) :: epsilon(:),a integer, intent(out) :: p integer, intent(in) :: m,mmax,pmax real(kind=large_real_kind), intent(out) :: sigmarho(:,:) real(kind=large_real_kind) :: gammat(mmax,mmax),Omega22(m,m),Lambda(mmax),& Omega12(m),Omega11,moyeps,sigeps,epsnorm(1:size(epsilon)) real(kind=large_real_kind) :: upsilon(size(epsilon)-mmax,mmax) integer :: i,n n=size(epsilon) if(mmax 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 real(kind=large_real_kind), intent(in) :: x(:,:) real(kind=large_real_kind), intent(out) :: a(:,:,:), s2(:,:) real(kind=large_real_kind) :: Gamma(0:size(a,1),size(x,2),size(x,2)) real(kind=large_real_kind) :: Delta(0:size(a,1),size(x,2),size(x,2)),Deltatilde(0:size(a,1),size(x,2),size(x,2)) real(kind=large_real_kind) :: V(0:size(a,1),size(x,2),size(x,2)),Vtilde(0:size(a,1),size(x,2),size(x,2)) real(kind=large_real_kind) :: Phi(0:size(a,1),0:size(a,1),size(x,2),size(x,2)) real(kind=large_real_kind) :: Phitilde(0:size(a,1),0:size(a,1),size(x,2),size(x,2)) integer :: k,nn,n,d,p,j n=size(x,1) d=size(x,2) p=size(a,1) Gamma=autocovmult(x,p) 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