Sae x Sae x 1: 1. {x (i) 0 0 }N i=1 (x (i) 0 0 p(x 0) ) 2. = 1,, T a d (a) i (i = 1,, N) I, II I. v (i) II. x (i) 1 = f (x (i) 1 1, v(i) (b) i (i = 1,

Size: px

Start display at page:

Download "Sae x Sae x 1: 1. {x (i) 0 0 }N i=1 (x (i) 0 0 p(x 0) ) 2. = 1,, T a d (a) i (i = 1,, N) I, II I. v (i) II. x (i) 1 = f (x (i) 1 1, v(i) (b) i (i = 1,"

くにもとかなり
6 years ago
Views:

1 ( ) 1 : ( ) Sampling/Imporance resampling (SIR) Kiagawa (1993, 1996), Gordon(1993) EnKF EnKF EnKF 1CPU 1core 2 x = f (x 1, v ) y = h (x, w ) (1a) (1b) PF p(x y 1 ) {x (i) 1 }N i=1, p(x y ) {x (i) }N i=1 : 1

2 Sae x Sae x 1: 1. {x (i) 0 0 }N i=1 (x (i) 0 0 p(x 0) ) 2. = 1,, T a d (a) i (i = 1,, N) I, II I. v (i) II. x (i) 1 = f (x (i) 1 1, v(i) (b) i (i = 1,, N) λ (i) /( N ) i=1 λ(i) (c) β (i) = λ (i) ) p(v ) = p(y x (i) 1 ) (d) {x (1) 1,, x(n) 1 } x(i) 1 β(i) {x (1),, x(n) } 2a 2b 2d ( 1) 2

3 3 3.1 run run (run) run 1. (MPI ) 2. 3

4 3. ( ) Forran 90 Forran 90 derived ype 3.2 p(x 0 ) N p(x 0 ) x 0 ( ) 1. 4

5 p(v ) 3.4 x (i) 1 1 ( )

6 3.5.1 x p(y x ) y = H x + w w R x (i) λ (i) 1 ( = p y x (i) 1 = ) [ 1 (2π)l R exp 1 2 ( y H x (i) 1 ) ( ) ] R 1 y H x (i) 1 Poisson β (i) PF 1. l (i) x (K) l (K) 2. ψ (i) = exp(l (i) l (K) ) 3. Ψ = i ψ(i) 4. β (i) = ψ (i) /Ψ subrouine subrouine xlamarr weigh 6

7 1: subrouine calc_weigh ( npoal, xlamarr, weigh ) implici none real (8), dimension ( npoal ) :: xlamarr, weigh real (8) :: xlammax, wsum ineger, inen ( in) :: npoal ineger :: i xlammax = maxval ( xlamarr )! xlamarr (i)! (= ) weigh (1)= exp ( xlamarr (1) - xlammax ) wsum = weigh (1) do i=2, npoal weigh (i)= exp ( xlamarr (i)- xlammax ) wsum = wsum + weigh (i)! 1 do i=1, npoal weigh (i)= weigh (i)/ wsum reurn end subrouine calc_weigh ( ) {x (i) 1 } β(i) N β (i) 7

8 N x (i) 1 x (i) 1 β(i) ( ) Nβ (i) (1 β (i) ) x (i) 1 {x(i) }N i=1} β (i) (1 β (i) )/N N 1 1. Nβ (i) m (i) i m(i) x (i) 1 m(i) 2. N x (i) 1 (Nβ(i) m (i) N i m(i) )/(N i m(i) ) 2 ( 2) 1. {ζ (i) } N i=0, {η (i) } N i=1 ζ(0) = 0, ζ (i) = (0 < ε 1) i j=1 β (j), η (i) = (i ε)/n 2. i = 1,, N ζ (i 1) η (j) < ζ (i) j x (j) := x(i) 1 8

9 (1) η (2) η (3) η (4) η (5) η (6) η (N) η 0 (1) ζ (2) ζ (3) ζ (4) ζ (5) ζ (N) ζ =1 2: Resampling 3: ( ) ( 3) subrouine subrouine 2 9

10 isrc 2: subrouine resample_map ( np, weigh, isrc ) implici none ineger, inen ( in) :: np ineger :: i,j, nzero real (8), parameer :: eps =0. 5 real (8), dimension (np), inen (in) :: weigh real (8), dimension ( np) :: zea, ea ineger, dimension ( np) :: isrc, nsample ineger, dimension ( np) :: izero!! zea zea (1)= weigh (1) do i=2,np -1 zea (i)= weigh (i)+ zea (i -1) zea (np )=1.0!! ea do i=1, np ea (i )=(i-eps )/ np nsample (:)=0! i izero (:)=0! nzero =0 j=1 do i=1, np if(j <= np) hen!! zea (i -1) -- zea (i) do while ( ea (j) < zea (i)) nsample (i)= nsample (i )+1 j=j+1 if(j > np) exi!! if( nsample (i )==0) hen nzero = nzero +1 izero ( nzero )=i do i=1, np isrc (i)=i! 10

11 !! j=1 do i=1, np do while ( nsample (i) > 1)! isrc ( izero (j ))= i! j=j+1 nsample (i)= nsample (i)-1 reurn end subrouine resample_map module sysem_model ineger, parameer :: ndim =50 ype, public :: sae_vecor real (8), dimension ( ndim ) :: a real (8), dimension ( ndim ) :: p end ype sae_vecor end module sysem_model sae vecor derived ype ype ( sae_vecor ) :: x(np) real (8), dimension ( np) :: weigh!! call resample_map (np, weigh, isrc ) do i=1, np if(i /= isrc (i)) hen x(i)=x( isrc (i)) isrc 11

12 Rank 0 Rank 1 Rank 2 Rank 3 4: MPI 4 MPI MPI (Message Passing Inerface) MPI ( 4) MPI MPI subrouine MPI MPI subrouine call MPI Forran C Forran C Forran MPI MPI subrouine MPI 12

13 MPI Web 1) 4.1 MPI MPI MPI subrouine program MPIprogram implici none include mpif.h ineger :: nprocs, irank ineger :: ierr 3: MPI call mpi_ini ( ierr ) call mpi_comm_size ( MPI_COMM_WORLD, nprocs, ierr ) call mpi_comm_rank ( MPI_COMM_WORLD, irank, ierr ) prin *, Num of proc :, nprocs, Rank :, irank call mpi_finalize ( ierr ) sop end program MPIprogram include mpif.h MPI subrouine mpif.h include call mpi ini(ierr) MPI subrouine call MPI ierr subrouine 0 1) hp://accc.riken.jp/hpc/raining/ex.hml 13

14 call mpi comm size(mpi COMM WORLD,nprocs,ierr) nprocs MPI COMM WORLD mpif.h call mpi comm rank(mpi COMM WORLD,irank,ierr) Rank 0, 1, 2,... irank rank call mpi finalize(ierr) MPI subrouine call MPI MPI 8 Num of proc : 8 Rank : 0 Num of proc : 8 Rank : 1 Num of proc : 8 Rank : 2 Num of proc : 8 Rank : 3 Num of proc : 8 Rank : 5 Num of proc : 8 Rank : 7 Num of proc : 8 Rank : 4 Num of proc : 8 Rank : mpi gaher subrouine ( 5) mpi gaher 4: mpi gaher program MPIgahTes implici none include mpif.h 14

15 Rank 0 Rank 1 Rank 2 MPI_GATHER Rank 3 5: MPI gaher ineger :: nprocs, irank, iroo ineger :: ierr, isa ineger :: nn ineger :: i ineger, parameer :: n =2 real (8) :: x(n) real (8), allocaable, dimension (:) :: z iroo =0 call mpi_ini ( ierr ) call mpi_comm_size ( MPI_COMM_WORLD, nprocs, ierr ) call mpi_comm_rank ( MPI_COMM_WORLD, irank, ierr ) do i=1,n x(i)=n* irank +i nn=n* nprocs allocae (z (1: nn), sa = isa ) call mpi_gaher (x,n, MPI_REAL8,z,n, MPI_REAL8, iroo,& MPI_COMM_WORLD, ierr ) if( irank == iroo ) hen do i=1, nn wrie (6, (a3,i3,2x,a6,f5.1) ) i:,i, z(i):,z(i) 15

16 end if deallocae (z, sa = isa ) call mpi_finalize ( ierr ) sop end program MPIgahTes call mpi gaher(x,n,mpi REAL8,z,n,MPI REAL8,iroo,... rank iroo 1 x: 2 n: 3 MPI REAL8: 8 ( MPI DOUBLE PRECISION ) MPI REAL4 4 MPI INTEGER 4 z: 2 5 n: MPI REAL8: 3 7 iroo: ( ) rank 4.3 ( ) mpi bcas subrouine mpi bcas 5: mpi bcas program MPIbroadcas implici none 16

17 include mpif.h ineger :: nprocs, irank, iroo ineger :: ierr ineger :: i ineger, parameer :: n =5 real (8) :: x(n) iroo =0 call mpi_ini ( ierr ) call mpi_comm_size ( MPI_COMM_WORLD, nprocs, ierr ) call mpi_comm_rank ( MPI_COMM_WORLD, irank, ierr ) if( irank == iroo ) do i=1,n x(i)=i else x (:)= -1.0 hen wrie (6, (a16,i3,2x,a3,5 f6.1) ) x:,(x(i),i=1,n) Before... Rank :,irank,& call mpi_bcas (x,n, MPI_REAL8, iroo, MPI_COMM_WORLD, ierr ) wrie (6, (a16,i3,2x,a3,5 f6.1) ) x:,(x(i),i=1,n) Afer... Rank :,irank,& call mpi_finalize ( ierr ) sop end program MPIbroadcas call mpi bcas(x,n,mpi REAL8,iroo,MPI COMM WORLD,ierr) rank iroo rank iroo 1 x: 2 n: 4 z: rank 4 17

18 Before... Rank : 0 x: Before... Rank : 1 x: Before... Rank : 2 x: Before... Rank : 3 x: Afer... Rank : 0 x: Afer... Rank : 1 x: Afer... Rank : 2 x: Afer... Rank : 3 x: mpi send mpi recv subrouine 6: MPI program MPIcomm implici none include mpif.h ineger :: mpisa ( MPI_STATUS_SIZE ) ineger :: nprocs, irank, iroo ineger :: ierr, isa, iag ineger :: isrc, ides ineger :: nn,i ineger, parameer :: n =2 real (8) :: x(n) isrc =0 ides =1 x (:)=0.0 call mpi_ini ( ierr ) call mpi_comm_size ( MPI_COMM_WORLD, nprocs, ierr ) call mpi_comm_rank ( MPI_COMM_WORLD, irank, ierr ) if( irank == isrc ) hen do i=1,n x(i)= dble (i) if( irank == ides ) hen wrie (6, (a10,2 f6.1) ) Before...,(x(i),i=1,n) 18

19 iag =0 if( irank == isrc ) hen call mpi_send (x,n, MPI_REAL8, ides,iag, MPI_COMM_WORLD,& ierr ) if( irank == ides ) hen call mpi_recv (x,n, MPI_REAL8,isrc,iag, MPI_COMM_WORLD,& mpisa, ierr ) if( irank == ides ) hen wrie (6, (a10,2 f6.1) ) Afer...,(x(i),i=1,n) call mpi_finalize ( ierr ) sop end program MPIcomm call mpi send(x,n,mpi REAL8,ides,iag,MPI COMM WORLD,ierr) rank ides 1 x: 2 n: 4 ides: rank 5 iag: ( ) call mpi recv(x,n,mpi REAL8,isrc,iag,MPI COMM WORLD,mpisa,ierr) Rank isrc 1 x: 2 n: 4 ides: rank 5 iag: mpi send 7 mpisa: MPI STATUS SIZE MPI STATUS SIZE mpif.h 19

20 mpi send mpi recv mpi isend mpi irecv subrouine mpi isend mpi irecv mpi wai subrouine call 4.5 (derived ype) derived ype mpi ype sruc mpi send mpi recv 1 ( )call MPI subrouine Forran MPI Inel Forran compiler 10.1 MPI Forran module ramparams use grid_parameers public 7: Forran (derived ype) ype, public :: drif_params sequence real :: ddp 20

21 real :: poen (nr,np) real :: vcl (nr,np),p1(nr,np) real :: pocoef (0:24,0: maximeseq ) end ype drif_params ype, public :: bound_params sequence real :: fb(np,nw,nk) real :: dens_par (8), empe_par (8) ineger :: ib0 ( np) end ype bound_params ype, public :: ram_daa sequence real :: f2(nr,np,nw,nk) real :: ype ( drif_ params ) :: d_ pars real :: hdens_fac (0: maximeseq ) ype ( bound_ params ) :: b_ pars ineger :: nid end ype ram_daa end module ramparams sequence Forran Forran MPI MPI Inel compiler subrouine 8: MPI subrouine gen_ype ( model_ype ) ineger, parameer :: nblock =13, ncomblock =12! In order o preserve he id for each insance,! he las block is no copied in resampling. ineger, dimension ( nblock ) :: block_ len, block_ displ, ypes ineger :: model_ype ineger, parameer :: bye1rd =8, bye1id =4! Block lengh for each elemen block_len (1)= nr*np*iw*ik block_len (2)=1 block_len (3)=1 block_len (4)= nr*np block_len (5)= nr*np block_len (6)= nr*np block_len (7)=25*( maximeseq +1) block_len (8)= maximeseq +1 21

22 block_len (9)= np*iw*ik block_len (10)=8 block_len (11)=8 block_len (12)= np block_len (13)=1! Type for each elemen do i =1,11 ypes (i)= MPI_REAL8 do i =12,13 ypes (i)= MPI_INTEGER! Block displacemen block_ displ ( 1) =0* bye1rd do i=2, nblock selec case ( ypes (i -1)) case ( MPI_REAL8 ) block_displ (i)= block_displ (i -1)+ bye1rd * block_len (i -1) case ( MPI_INTEGER ) block_displ (i)= block_displ (i -1)+ bye1id * block_len (i -1) case ( MPI_CHARACTER ) block_displ (i)= block_displ (i -1)+ block_len (i -1) end selec call mpi_ ype_ sruc ( ncomblock, block_ len, block_ displ, ypes,& model_ype, ierr ) call mpi_ ype_ commi ( model_ ype, ierr )! call mpi_ ype_ exen ( model_ ype, iex, ierr ) reurn end subrouine gen_ype call mpi ype sruc 1 ncomblock: 2 block len: ncomblock 3 block displ: ncomblock 1 22

23 4 ypes: ncomblock MPI REAL8, MPI INTEGER 5 model ype: call mpi ype commi mpi ype sruc model ype call MPI 3 block displ Inel compiler derived ype sequence sequenial block displ mpi send mpi recv ype ( ram_daa ) :: ram!! () 9: MPI if( irank == isrc ) hen call mpi_send (ram,1, model_ype, ides,iag, MPI_COMM_WORLD, ierr ) if( irank == ides ) hen call mpi_recv (ram,1, model_ype,isrc,iag, MPI_COMM_WORLD,& mpisa, ierr ) MPI MPI subrouine subrouine mpi send call 23

24 5 x (i) 1 MPI MPI 5.1 MPI x (i) 1 MPI 24

25 5.2 λ (i) N i=1 λ(i) mpi gaher 1) npoal npoal 1 m rank iroo 10: ineger :: prec_real = MPI_REAL8 real (8) :: xlambda (m) real (8) :: xlamarr ( npoal )!! () call mpi_gaher ( xlambda,m,prec, xlamarr,m,prec, iroo,& MPI_COMM_WORLD, ierr ) xlambda mpi bcas 5.4 ( ) 1) mpi reduce subrouine 25

26 mpi send mpi recv rouine mpi send call rouine mpi recv call isrc rank rank rank rank rank rank mpi send call rank mpi recv call 11: call mpi_comm_size ( MPI_COMM_WORLD, nprocs, ierr ) call mpi_comm_rank ( MPI_COMM_WORLD, irank, ierr ) iag =0 do k=0, nprocs -1! rank do m=1, npeach! ides =m+k* npeach!! ksrc : ides! msrc : ides rank ksrc =( isrc ( ides ) -1)/ npeach msrc = isrc ( ides )- ksrc * npeach rans =. rue.! rank if(k== ksrc ) hen rans =. false.! Rank irank if(k== irank ) hen xp(m)= xp( msrc ) else! m > 1 if( m > 1) hen do mm =1,m -1 if( isrc ( ides )== isrc (mm+k* npeach )) hen rans =. false. if( irank ==k) xp(m)= xp(mm) 26

27 exi! if( rans ) hen! Rank ksrc if( irank == ksrc ) hen call mpi_send (xp( msrc ),1, model_ype,k,iag,& MPI_COMM_WORLD, ierr )! Rank k else if( irank.eq.k) hen call mpi_recv (xp(m),1, model_ype,ksrc,iag,& MPI_COMM_WORLD, mpisa, ierr )! rank iag = iag +1 if( iag > MPI_TAG_UB ) iag =0 27

120802_MPI.ppt

120802_MPI.ppt CPU CPU CPU CPU CPU SMP Symmetric MultiProcessing CPU CPU CPU CPU CPU CPU CPU CPU CPU CPU CPU CPU CP OpenMP MPI MPI CPU CPU CPU CPU CPU CPU CPU CPU CPU CPU MPI MPI+OpenMP CPU CPU CPU CPU CPU CPU CPU CP