source: trunk/WRF.COMMON/WRFV2/external/RSL_LITE/period.c @ 3022

Last change on this file since 3022 was 11, checked in by aslmd, 15 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
  • Property svn:executable set to *
File size: 16.2 KB
Line 
1#include <stdio.h>
2#include <fcntl.h>
3
4#define STANDARD_ERROR 2
5
6#define STANDARD_OUTPUT 1
7
8#include "mpi.h"
9#include "rsl_lite.h"
10
11#define F_PACK
12
13static int yp_curs, ym_curs, xp_curs, xm_curs ;
14
15RSL_LITE_INIT_PERIOD ( 
16                int * Fcomm0,
17                int * shw0,
18                int * n3dR0, int *n2dR0, int * typesizeR0 , 
19                int * n3dI0, int *n2dI0, int * typesizeI0 , 
20                int * n3dD0, int *n2dD0, int * typesizeD0 , 
21                int * n3dL0, int *n2dL0, int * typesizeL0 , 
22                int * me0, int * np0 , int * np_x0 , int * np_y0 ,
23                int * ips0 , int * ipe0 , int * jps0 , int * jpe0 , int * kps0 , int * kpe0 )
24{
25  int n3dR, n2dR, typesizeR ;
26  int n3dI, n2dI, typesizeI ;
27  int n3dD, n2dD, typesizeD ;
28  int n3dL, n2dL, typesizeL ;
29  int shw ;
30  int me, np, np_x, np_y ;
31  int ips , ipe , jps , jpe , kps , kpe ;
32  int yp, ym, xp, xm ;
33  int nbytes ;
34  int coords[2] ;
35  MPI_Comm comm, *comm0, dummy_comm ;
36
37  comm0 = &dummy_comm ;
38  *comm0 = MPI_Comm_f2c( *Fcomm0 ) ;
39
40  shw = *shw0 ;
41  n3dR = *n3dR0 ; n2dR = *n2dR0 ; typesizeR = *typesizeR0 ;
42  n3dI = *n3dI0 ; n2dI = *n2dI0 ; typesizeI = *typesizeI0 ;
43  n3dD = *n3dD0 ; n2dD = *n2dD0 ; typesizeD = *typesizeD0 ;
44  n3dL = *n3dL0 ; n2dL = *n2dL0 ; typesizeL = *typesizeL0 ;
45  me = *me0 ; np = *np0 ; np_x = *np_x0 ; np_y = *np_y0 ;
46  ips = *ips0-1 ; ipe = *ipe0-1 ; jps = *jps0-1 ; jpe = *jpe0-1 ; kps = *kps0-1 ; kpe = *kpe0-1 ;
47
48#if 0
49  if ( np_y > 1 ) {
50    nbytes = typesizeR*(ipe-ips+1+2*shw)*(shw+1)*(n3dR*(kpe-kps+1)+n2dR) +
51             typesizeI*(ipe-ips+1+2*shw)*(shw+1)*(n3dI*(kpe-kps+1)+n2dI) +
52             typesizeD*(ipe-ips+1+2*shw)*(shw+1)*(n3dD*(kpe-kps+1)+n2dD) +
53             typesizeL*(ipe-ips+1+2*shw)*(shw+1)*(n3dL*(kpe-kps+1)+n2dL) ;
54    MPI_Cart_shift( *comm0, 1, 1, &ym, &yp ) ;
55    if ( yp != MPI_PROC_NULL ) {
56       buffer_for_proc ( yp , nbytes, RSL_RECVBUF ) ;
57       buffer_for_proc ( yp , nbytes, RSL_SENDBUF ) ;
58    }
59    if ( ym != MPI_PROC_NULL ) {
60       buffer_for_proc ( ym , nbytes, RSL_RECVBUF ) ;
61       buffer_for_proc ( ym , nbytes, RSL_SENDBUF ) ;
62    }
63  }
64#endif
65  if ( np_x > 1 ) {
66    nbytes = typesizeR*(jpe-jps+1+2*shw)*(shw+1)*(n3dR*(kpe-kps+1)+n2dR) +
67             typesizeI*(jpe-jps+1+2*shw)*(shw+1)*(n3dI*(kpe-kps+1)+n2dI) +
68             typesizeD*(jpe-jps+1+2*shw)*(shw+1)*(n3dD*(kpe-kps+1)+n2dD) +
69             typesizeL*(jpe-jps+1+2*shw)*(shw+1)*(n3dL*(kpe-kps+1)+n2dL) ;
70
71/*
72 this assumes that the topoology associated with the communicator is periodic
73 the period routines should be called with "local_communicator_periodic", which
74 is set up in module_dm.F for RSL_LITE.  Registry generated code automaticall
75 does this (gen_comms.c for RSL_LITE).
76*/
77
78    MPI_Comm_rank( *comm0, &me ) ;
79    MPI_Cart_coords( *comm0, me, 2, coords ) ;
80    MPI_Cart_shift( *comm0, 1, 1, &xm, &xp ) ;
81    if ( xm != MPI_PROC_NULL && coords[1] == np_x - 1 ) { /* process on right hand side of mesh */
82       buffer_for_proc ( xp , nbytes, RSL_RECVBUF ) ;
83       buffer_for_proc ( xp , nbytes, RSL_SENDBUF ) ;
84    }
85    if ( xp != MPI_PROC_NULL && coords[1] == 0 ) {        /* process on left hand side of mesh */
86       buffer_for_proc ( xm,  nbytes, RSL_RECVBUF ) ;
87       buffer_for_proc ( xm , nbytes, RSL_SENDBUF ) ;
88    }
89  }
90  yp_curs = 0 ; ym_curs = 0 ; xp_curs = 0 ; xm_curs = 0 ;
91}
92
93RSL_LITE_PACK_PERIOD_X ( int* Fcomm0, char * buf , int * shw0 , int * typesize0 , int * xy0 , int * pu0 , char * memord , int * xstag0 ,
94           int *me0, int * np0 , int * np_x0 , int * np_y0 , 
95           int * ids0 , int * ide0 , int * jds0 , int * jde0 , int * kds0 , int * kde0 ,
96           int * ims0 , int * ime0 , int * jms0 , int * jme0 , int * kms0 , int * kme0 ,
97           int * ips0 , int * ipe0 , int * jps0 , int * jpe0 , int * kps0 , int * kpe0 )
98{
99  int me, np, np_x, np_y ;
100  int shw , typesize ;
101  int ids , ide , jds , jde , kds , kde ;
102  int ims , ime , jms , jme , kms , kme ;
103  int ips , ipe , jps , jpe , kps , kpe ;
104  int xstag ;  /* 0 not stag, 1 stag */
105  int xy ;   /* y = 0 , x = 1 */
106  int pu ;   /* pack = 0 , unpack = 1 */
107  register int i, j, k, t ;
108#ifdef crayx1
109  register int i2,i3,i4,i_offset;
110#endif
111  char *p ;
112  int da_buf ;
113  int yp, ym, xp, xm ;
114  int nbytes, ierr ;
115  register int *pi, *qi ;
116  int coords[2], dims[2] ;
117  int js, je, ks, ke, is, ie, wcount ;
118float f ;
119  MPI_Comm comm, *comm0, dummy_comm ;
120
121  comm0 = &dummy_comm ;
122  *comm0 = MPI_Comm_f2c( *Fcomm0 ) ;
123
124  me = *me0 ; np = *np0 ; np_x = *np_x0 ; np_y = *np_y0 ;
125  xstag = *xstag0 ;
126  shw = *shw0 ; typesize = *typesize0 ;
127  ids = *ids0-1 ; ide = *ide0-1 ; jds = *jds0-1 ; jde = *jde0-1 ; kds = *kds0-1 ; kde = *kde0-1 ;
128  ims = *ims0-1 ; ime = *ime0-1 ; jms = *jms0-1 ; jme = *jme0-1 ; kms = *kms0-1 ; kme = *kme0-1 ;
129  ips = *ips0-1 ; ipe = *ipe0-1 ; jps = *jps0-1 ; jpe = *jpe0-1 ; kps = *kps0-1 ; kpe = *kpe0-1 ;
130  xy = *xy0 ;
131  pu = *pu0 ;
132
133  dims[0] = np_x ;
134  dims[1] = np_y ;
135
136/* need to adapt for other memory orders */
137
138#define RANGE(S1,E1,S2,E2,S3,E3,S4,E4) (((E1)-(S1)+1)*((E2)-(S2)+1)*((E3)-(S3)+1)*((E4)-(S4)+1))
139#define IMAX(A) (((A)>ids)?(A):ids)
140#define IMIN(A) (((A)<ide)?(A):ide)
141#define JMAX(A) (((A)>jds)?(A):jds)
142#define JMIN(A) (((A)<jde)?(A):jde)
143
144  da_buf = ( pu == 0 ) ? RSL_SENDBUF : RSL_RECVBUF ;
145
146  if ( np_x > 1 && xy == 1 ) {
147    MPI_Comm_rank( *comm0, &me ) ;
148    MPI_Cart_coords( *comm0, me, 2, coords ) ;
149    MPI_Cart_shift( *comm0, 1, 1, &xm, &xp ) ;
150    if ( coords[1] == np_x - 1 ) {                /* process on right hand edge of domain */
151      p = buffer_for_proc( xp , 0 , da_buf ) ;
152      if ( pu == 0 ) {
153        js = JMAX(jps-shw) ; je = JMIN(jpe+shw) ;
154        ks = kps           ; ke = kpe ;
155        is = ipe-shw       ; ie = ipe-1         ;
156        nbytes = buffer_size_for_proc( xp , da_buf ) ;
157        if ( xp_curs + RANGE( JMAX(jps-shw), JMIN(jpe+shw), kps, kpe, ipe-shw, ipe-1, 1, typesize ) > nbytes ) {
158          fprintf(stderr,"memory overwrite in rsl_lite_pack_period_x, right hand X to %d, %d > %d\n",xp,
159              xp_curs + RANGE( JMAX(jps-shw), JMIN(jpe+shw), kps, kpe, ipe-shw, ipe-1, 1, typesize ), nbytes ) ;
160          MPI_Abort(MPI_COMM_WORLD, 98) ;
161        }
162        if ( typesize == sizeof(long int) ) {
163          F_PACK_LINT ( buf, p+xp_curs, &js, &je, &ks, &ke, &is, &ie,
164                                        &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
165          xp_curs += wcount*typesize ;
166        } else
167        if ( typesize == sizeof(int) ) {
168#ifdef F_PACK
169          F_PACK_INT ( buf, p+xp_curs, &js, &je, &ks, &ke, &is, &ie,
170                                       &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
171          xp_curs += wcount*typesize ;
172#else
173          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
174            for ( k = kps ; k <= kpe ; k++ ) {
175              pi = (int *)(p+xp_curs) ;
176              i = ipe-shw ;
177              qi = (int *)((buf + typesize*( (i-ims) + (ime-ims+1)*(
178                                             (k-kms) + (j-jms)*(kme-kms+1))))) ;
179              for ( i = ipe-shw ; i <= ipe-1 ; i++ ) {
180                *pi++ = *qi++ ;
181                xp_curs += typesize ;
182              }
183            }
184          }
185#endif
186        }
187        else {
188          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
189            for ( k = kps ; k <= kpe ; k++ ) {
190              for ( i = ipe-shw ; i <= ipe-1 ; i++ ) {
191                for ( t = 0 ; t < typesize ; t++ ) {
192                  *(p+xp_curs) = 
193                                 *(buf + t + typesize*(
194                                        (i-ims) + (ime-ims+1)*(
195                                        (k-kms) + (j-jms)*(kme-kms+1))) ) ;
196                  xp_curs++ ;
197                }
198              }
199            }
200          }
201        }
202      } else {
203        js = JMAX(jps-shw) ; je = JMIN(jpe+shw) ;
204        ks = kps           ; ke = kpe ;
205        is = ipe           ; ie = ipe+shw-1+xstag ;
206        if ( typesize == sizeof(long int) ) {
207          F_UNPACK_LINT ( p+xp_curs, buf, &js, &je, &ks, &ke, &is, &ie,
208                                          &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
209          xp_curs += wcount*typesize ;
210        } else
211        if ( typesize == sizeof(int) ) {
212#ifdef F_PACK
213          F_UNPACK_INT ( p+xp_curs, buf, &js, &je, &ks, &ke, &is, &ie,
214                                         &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
215          xp_curs += wcount*typesize ;
216#else
217          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
218            for ( k = kps ; k <= kpe ; k++ ) {
219              pi = (int *)(p+xp_curs) ;
220              i = ipe ;
221              qi = (int *)((buf + typesize*( (i-ims) + (ime-ims+1)*(
222                                             (k-kms) + (j-jms)*(kme-kms+1))))) ;
223              for ( i = ipe ; i <= ipe+shw-1+xstag ; i++ ) {
224                *qi++ = *pi++ ;
225                xp_curs += typesize ;
226              }
227            }
228          }
229#endif
230        }
231        else {
232          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
233            for ( k = kps ; k <= kpe ; k++ ) {
234              for ( i = ipe ; i <= ipe+shw-1+xstag ; i++ ) {
235                for ( t = 0 ; t < typesize ; t++ ) {
236                                 *(buf + t + typesize*(
237                                        (i-ims) + (ime-ims+1)*(
238                                        (k-kms) + (j-jms)*(kme-kms+1))) ) =
239                  *(p+xp_curs) ;
240                  xp_curs++ ;
241                }
242              }
243            }
244          }
245        }
246      }
247    }
248    if ( coords[1] == 0 ) {         /* process on left hand edge of domain */
249      p = buffer_for_proc( xm , 0 , da_buf ) ;
250      if ( pu == 0 ) {
251        js = JMAX(jps-shw) ; je = JMIN(jpe+shw) ;
252        ks = kps           ; ke = kpe ;
253        is = ips           ; ie = ips+shw-1+xstag ;
254        nbytes = buffer_size_for_proc( xm , da_buf ) ;
255        if ( xm_curs + RANGE( JMAX(jps-shw), JMIN(jpe+shw), kps, kpe, ips, ips+shw-1+xstag, 1, typesize ) > nbytes ) {
256          fprintf(stderr,"memory overwrite in rsl_lite_pack,  left hand X to %d , %d > %d\n",xm,
257              xm_curs + RANGE( JMAX(jps-shw), JMIN(jpe+shw), kps, kpe, ips, ips+shw-1+xstag, 1, typesize ), nbytes ) ;
258          MPI_Abort(MPI_COMM_WORLD, 98) ;
259        }
260        if ( typesize == sizeof(long int) ) {
261          F_PACK_LINT ( buf, p+xm_curs, &js, &je, &ks, &ke, &is, &ie,
262                                        &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
263          xm_curs += wcount*typesize ;
264        } else
265        if ( typesize == sizeof(int) ) {
266#ifdef F_PACK
267          F_PACK_INT ( buf, p+xm_curs, &js, &je, &ks, &ke, &is, &ie,
268                                       &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
269          xm_curs += wcount*typesize ;
270#else
271          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
272            for ( k = kps ; k <= kpe ; k++ ) {
273              pi = (int *)(p+xm_curs) ;
274              i = ips ;
275              qi = (int *)((buf + typesize*( (i-ims) + (ime-ims+1)*(
276                                             (k-kms) + (j-jms)*(kme-kms+1))))) ;
277              for ( i = ips ; i <= ips+shw-1+xstag ; i++ ) {
278                *pi++ = *qi++ ;
279                xm_curs += typesize ;
280              }
281            }
282          }
283#endif
284        }
285        else {
286          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
287            for ( k = kps ; k <= kpe ; k++ ) {
288              for ( i = ips ; i <= ips+shw-1+xstag ; i++ ) {
289                for ( t = 0 ; t < typesize ; t++ ) {
290                  *(p+xm_curs) = 
291                                 *(buf + t + typesize*(
292                                        (i-ims) + (ime-ims+1)*(
293                                        (k-kms) + (j-jms)*(kme-kms+1))) ) ;
294                  xm_curs++ ;
295                }
296              }
297            }
298          }
299        }
300      } else {
301        js = JMAX(jps-shw) ; je = JMIN(jpe+shw) ;
302        ks = kps           ; ke = kpe ;
303        is = ips-shw       ; ie = ips-1           ;
304        if ( typesize == sizeof(long int) ) {
305          F_UNPACK_LINT ( p+xm_curs, buf, &js, &je, &ks, &ke, &is, &ie,
306                                          &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
307          xm_curs += wcount*typesize ;
308        } else
309        if ( typesize == sizeof(int) ) {
310#ifdef F_PACK
311          F_UNPACK_INT ( p+xm_curs, buf, &js, &je, &ks, &ke, &is, &ie,
312                                         &jms,&jme,&kms,&kme,&ims,&ime, &wcount ) ;
313          xm_curs += wcount*typesize ;
314#else
315          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
316            for ( k = kps ; k <= kpe ; k++ ) {
317              pi = (int *)(p+xm_curs) ;
318              i = ips-shw ;
319              qi = (int *)((buf + typesize*( (i-ims) + (ime-ims+1)*(
320                                             (k-kms) + (j-jms)*(kme-kms+1))))) ;
321              for ( i = ips-shw ; i <= ips-1 ; i++ ) {
322                *qi++ = *pi++ ;
323                xm_curs += typesize ;
324              }
325            }
326          }
327#endif
328        }
329        else {
330          for ( j = JMAX(jps-shw) ; j <= JMIN(jpe+shw) ; j++ ) {
331            for ( k = kps ; k <= kpe ; k++ ) {
332              for ( i = ips-shw ; i <= ips-1 ; i++ ) {
333                for ( t = 0 ; t < typesize ; t++ ) {
334                                 *(buf + t + typesize*(
335                                        (i-ims) + (ime-ims+1)*(
336                                        (k-kms) + (j-jms)*(kme-kms+1))) ) =
337                  *(p+xm_curs) ;
338                  xm_curs++ ;
339                }
340              }
341            }
342          }
343        }
344      }
345    }
346  }
347}
348
349static MPI_Request yp_recv, ym_recv, yp_send, ym_send ;
350static MPI_Request xp_recv, xm_recv, xp_send, xm_send ;
351
352#if 0
353RSL_LITE_EXCH_PERIOD_Y ( int * Fcomm0, int *me0, int * np0 , int * np_x0 , int * np_y0 )
354{
355  int me, np, np_x, np_y ;
356  int yp, ym, xp, xm, ierr ;
357  MPI_Status stat ;
358  MPI_Comm comm, *comm0, dummy_comm ;
359
360  comm0 = &dummy_comm ;
361  *comm0 = MPI_Comm_f2c( *Fcomm0 ) ;
362#if 1
363  comm = *comm0 ; me = *me0 ; np = *np0 ; np_x = *np_x0 ; np_y = *np_y0 ;
364  if ( np_y > 1 ) {
365    yp = me + np_x ; ym = me - np_x ;
366    if ( yp >= 0 && yp < np ) {
367      ierr=MPI_Irecv ( buffer_for_proc( yp, yp_curs, RSL_RECVBUF ), yp_curs, MPI_CHAR, yp, me, comm, &yp_recv ) ;
368    }
369    if ( ym >= 0 && ym < np ) {
370      ierr=MPI_Irecv ( buffer_for_proc( ym, ym_curs, RSL_RECVBUF ), ym_curs, MPI_CHAR, ym, me, comm, &ym_recv ) ;
371    }
372    if ( yp >= 0 && yp < np ) {
373      ierr=MPI_Isend ( buffer_for_proc( yp, 0,       RSL_SENDBUF ), yp_curs, MPI_CHAR, yp, yp, comm, &yp_send ) ;
374    }
375    if ( ym >= 0 && ym < np ) {
376      ierr=MPI_Isend ( buffer_for_proc( ym, 0,       RSL_SENDBUF ), ym_curs, MPI_CHAR, ym, ym, comm, &ym_send ) ;
377    }
378    if ( yp >= 0 && yp < np ) MPI_Wait( &yp_recv, &stat ) ;
379    if ( ym >= 0 && ym < np ) MPI_Wait( &ym_recv, &stat ) ;
380    if ( yp >= 0 && yp < np ) MPI_Wait( &yp_send, &stat ) ;
381    if ( ym >= 0 && ym < np ) MPI_Wait( &ym_send, &stat ) ;
382  }
383  yp_curs = 0 ; ym_curs = 0 ; xp_curs = 0 ; xm_curs = 0 ;
384#else
385fprintf(stderr,"RSL_LITE_EXCH_Y disabled\n") ;
386#endif
387}
388#endif
389
390RSL_LITE_EXCH_PERIOD_X ( int * Fcomm0, int *me0, int * np0 , int * np_x0 , int * np_y0 )
391{
392  int me, np, np_x, np_y ;
393  int yp, ym, xp, xm, nbytes ;
394  MPI_Status stat ;
395  MPI_Comm comm, *comm0, dummy_comm ;
396  int coords[2], dims[2] ;
397
398  comm0 = &dummy_comm ;
399  *comm0 = MPI_Comm_f2c( *Fcomm0 ) ;
400#if 1
401  comm = *comm0 ; me = *me0 ; np = *np0 ; np_x = *np_x0 ; np_y = *np_y0 ;
402  dims[0] = np_x ;
403  dims[1] = np_y ;
404
405  if ( np_x > 1 ) {
406    MPI_Comm_rank( *comm0, &me ) ;
407    MPI_Cart_coords( *comm0, me, 2, coords ) ;
408    MPI_Cart_shift( *comm0, 1, 1, &xm, &xp ) ;
409    if ( coords[1] == np_x - 1 ) {   /* proc on right hand side of domain */
410      nbytes = buffer_size_for_proc( xp, RSL_RECVBUF ) ;
411      MPI_Irecv ( buffer_for_proc( xp , xp_curs, RSL_RECVBUF ), nbytes, MPI_CHAR, xp, me, comm, &xp_recv ) ;
412    }
413    if ( coords[1] == 0 ) {          /* proc on left hand side of domain */
414      nbytes = buffer_size_for_proc( xm, RSL_RECVBUF ) ;
415      MPI_Irecv ( buffer_for_proc( xm, xm_curs, RSL_RECVBUF ), nbytes, MPI_CHAR, xm, me, comm, &xm_recv ) ;
416    }
417    if ( coords[1] == np_x - 1 ) {   /* proc on right hand side of domain */
418      MPI_Isend ( buffer_for_proc( xp , 0,       RSL_SENDBUF ), xp_curs, MPI_CHAR, xp, xp, comm, &xp_send ) ;
419    }
420    if ( coords[1] == 0 ) {          /* proc on left hand side of domain */
421      MPI_Isend ( buffer_for_proc( xm, 0,       RSL_SENDBUF ), xm_curs, MPI_CHAR, xm, xm, comm, &xm_send ) ;
422    }
423    if ( coords[1] == np_x - 1 ) MPI_Wait( &xp_recv, &stat ) ; 
424    if ( coords[1] == 0        ) MPI_Wait( &xm_recv, &stat ) ; 
425    if ( coords[1] == np_x - 1 ) MPI_Wait( &xp_send, &stat ) ; 
426    if ( coords[1] == 0        ) MPI_Wait( &xm_send, &stat ) ;
427  }
428#else
429fprintf(stderr,"RSL_LITE_EXCH_X disabled\n") ;
430#endif
431  yp_curs = 0 ; ym_curs = 0 ; xp_curs = 0 ; xm_curs = 0 ;
432}
433
Note: See TracBrowser for help on using the repository browser.