wrfv2_fire/external/RSL_LITE/task_for_point.c

   1 #ifndef MS_SUA
   2 # include <stdio.h>
   3 #endif
   4 #include "rsl_lite.h"
   5
   6 /* updated 20051021, new algorithm distributes the remainder, if any, at either ends of the dimension
   7    rather than the first remainder number of processors in the dimension. Idea is that the processes
   8    on the ends have less work because they're boundary processes.  New alg works like this:
   9                      a                         b
  10          + + + + + + o o o o o o o o o o o o o + + + + + +
  11
  12    + represents a process with an extra point (npoints is n/p+1), o processors that don't (n/p)
  13    a and b are the starting process indices in the dimension of the new section of o or x.
  14    JM
  15 */
  16
  17 static char tfpmess[1024] ;
  18
  19 TASK_FOR_POINT ( i_p , j_p , ids_p, ide_p , jds_p, jde_p , npx_p , npy_p , Px_p, Py_p , minx_p, miny_p, ierr_p )
  20   int_p i_p , j_p , Px_p , Py_p , ids_p, ide_p , jds_p, jde_p , npx_p , npy_p, minx_p, miny_p, ierr_p ;
  21 {
  22   int i , j , ids, ide, jds, jde, npx, npy, minx, miny ;  /* inputs */
  23   int Px, Py ;                                            /* output */
  24   int idim, jdim ;
  25   int rem, a, b ;
  26   i = *i_p - 1 ;
  27   j = *j_p - 1 ;
  28   npx = *npx_p ;
  29   npy = *npy_p ;
  30   minx = *minx_p ;
  31   miny = *miny_p ;
  32   ids = *ids_p - 1 ; ide = *ide_p - 1 ;
  33   jds = *jds_p - 1 ; jde = *jde_p - 1 ;
  34   idim = ide - ids + 1 ;
  35   jdim = jde - jds + 1 ;
  36
  37   *ierr_p = 0 ;
  38
  39   /* begin: jm for Peter Johnsen -- noticed problem with polar filters in gwrf
  40      if the number of processors exceeds number of vertical levels */
  41   if ( npx > idim ) { npx = idim ; }
  42   if ( npy > jdim ) { npy = jdim ; }
  43
  44   /* begin: wig; 10-Mar-2008
  45     Check that the number of processors is not so high that the halos begin to overlap.
  46     If they do, then reduce the number of processors allowed for that dimension.
  47   */
  48   tfpmess[0] = '\0' ;
  49   if ( idim / npx < minx ) {
  50     npx = idim/minx ;
  51     if (npx < 1) { npx = 1 ;}
  52     if (npx != *npx_p) {
  53       sprintf(tfpmess,"RSL_LITE: TASK_FOR_POINT LIMITING PROCESSOR COUNT IN X-DIRECTION TO %d %d\n", npx,*npx_p) ;
  54       *ierr_p = 1 ;
  55     }
  56   }
  57   if ( jdim / npy < miny ) {
  58     npy = jdim/miny ;
  59     if (npy < 1) { npy = 1 ;}
  60     if (npy != *npy_p) {
  61       sprintf(tfpmess,"RSL_LITE: TASK_FOR_POINT LIMITING PROCESSOR COUNT IN Y-DIRECTION TO %d %d\n", npy,*npy_p) ;
  62       *ierr_p = 1 ;
  63     }
  64   }
  65   /* end: wig */
  66
  67   i = i >= ids ? i : ids ; i = i <= ide ? i : ide ;
  68   rem = idim % npx ;
  69   a = ( rem / 2 ) * ( (idim / npx) + 1 ) ;
  70   b = a + ( npx - rem ) * ( idim / npx ) ;
  71   if ( i-ids < a ) {
  72     Px = (i-ids) / ( (idim / npx) + 1 ) ;
  73   }
  74   else if ( i-ids < b ) {
  75     Px = ( a / ( (idim / npx) + 1 ) ) + (i-a-ids) / ( ( b - a ) / ( npx - rem ) )     ;
  76   }
  77   else {
  78     Px = ( a / ( (idim / npx) + 1 ) ) + (b-a-ids) / ( ( b - a ) / ( npx - rem ) ) +
  79                                         (i-b-ids) / ( ( idim / npx ) + 1 )  ;
  80   }
  81
  82   j = j >= jds ? j : jds ; j = j <= jde ? j : jde ;
  83   rem = jdim % npy ;
  84   a = ( rem / 2 ) * ( (jdim / npy) + 1 ) ;
  85   b = a + ( npy - rem ) * ( jdim / npy ) ;
  86   if ( j-jds < a ) {
  87     Py = (j-jds) / ( (jdim / npy) + 1 ) ;
  88   }
  89   else if ( j-jds < b ) {
  90     Py = ( a / ( (jdim / npy) + 1 ) ) + (j-a-jds) / ( ( b - a ) / ( npy - rem ) )     ;
  91   }
  92   else {
  93     Py = ( a / ( (jdim / npy) + 1 ) ) + (b-a-jds) / ( ( b - a ) / ( npy - rem ) ) +
  94                                         (j-b-jds) / ( ( jdim / npy ) + 1 )  ;
  95   }
  96
  97   *Px_p = Px ;
  98   *Py_p = Py ;
  99 }
 100
 101 TASK_FOR_POINT_MESSAGE()
 102 {
 103   fprintf(stderr,"%s\n",tfpmess) ;
 104 }
 105
 106 #if 0
 107 main()
 108 {
 109   int ips[100], ipe[100] ;
 110   int jps[100], jpe[100] ;
 111   int shw, i , j , ids, ide, jds, jde, npx, npy ;  /* inputs */
 112   int Px, Py, P ;                             /* output */
 113   printf("i, j, ids, ide, jds, jde, npx, npy\n") ;
 114   scanf("%d %d %d %d %d %d %d %d",&i, &j, &ids,&ide,&jds,&jde,&npx,&npy ) ;
 115   shw =0 ;
 116   for ( i = 0 ; i < 100 ; i++ ) { ips[i] = 9999999 ; ipe[i] = -99999999 ; }
 117   for ( i = 0 ; i < 100 ; i++ ) { jps[i] = 9999999 ; jpe[i] = -99999999 ; }
 118 #if 1
 119   for ( j = jds-shw ; j <= jde+shw ; j++ )
 120   {
 121   for ( i = ids-shw ; i <= ide+shw ; i++ )
 122   {
 123 #endif
 124   TASK_FOR_POINT ( &i , &j ,
 125                    &ids, &ide, &jds, &jde , &npx , &npy ,
 126                    &Px, &Py ) ;
 127   printf("(%3d %3d)   ",Px,Py) ;
 128 #if 1
 129   }
 130   printf("\n") ;
 131   }
 132 /*  for ( i = 0 ; i < npx*npy ; i++ ) { */
 133 /*    fprintf(stderr,"%3d. ips %d ipe %d (%d) jps %d jpe %d (%d)\n", i, ips[i], ipe[i], ipe[i]-ips[i]+1, jps[i], jpe[i], jpe[i]-jps[i]+1 ) ; */
 134 /*  } */
 135 #endif
 136 }
 137 #endif