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 TASK_FOR_POINT ( i_p , j_p , ids_p, ide_p , jds_p, jde_p , npx_p , npy_p , Px_p, Py_p )
18   int_p i_p , j_p , Px_p , Py_p , ids_p, ide_p , jds_p, jde_p , npx_p , npy_p ;
19 {
20   int i , j , ids, ide, jds, jde, npx, npy ;  /* inputs */
21   int Px, Py ;                                /* output */
22   int idim, jdim ;
23   int rem, a, b ;
24   i = *i_p - 1 ;
25   j = *j_p - 1 ;
26   npx = *npx_p ;
27   npy = *npy_p ;
28   ids = *ids_p - 1 ; ide = *ide_p - 1 ;
29   jds = *jds_p - 1 ; jde = *jde_p - 1 ;
30   idim = ide - ids + 1 ;
31   jdim = jde - jds + 1 ;
32 
33   i = i >= ids ? i : ids ; i = i <= ide ? i : ide ;
34   rem = idim % npx ;
35   a = ( rem / 2 ) * ( (idim / npx) + 1 ) ;
36   b = a + ( npx - rem ) * ( idim / npx ) ;
37   if ( i-ids < a ) {
38     Px = (i-ids) / ( (idim / npx) + 1 ) ;
39   }
40   else if ( i-ids < b ) {
41     Px = ( a / ( (idim / npx) + 1 ) ) + (i-a-ids) / ( ( b - a ) / ( npx - rem ) )     ;
42   }
43   else {
44     Px = ( a / ( (idim / npx) + 1 ) ) + (b-a-ids) / ( ( b - a ) / ( npx - rem ) ) +
45                                         (i-b-ids) / ( ( idim / npx ) + 1 )  ;
46   }
47 
48   j = j >= jds ? j : jds ; j = j <= jde ? j : jde ;
49   rem = jdim % npy ;
50   a = ( rem / 2 ) * ( (jdim / npy) + 1 ) ;
51   b = a + ( npy - rem ) * ( jdim / npy ) ;
52   if ( j-jds < a ) {
53     Py = (j-jds) / ( (jdim / npy) + 1 ) ;
54   }
55   else if ( j-jds < b ) {
56     Py = ( a / ( (jdim / npy) + 1 ) ) + (j-a-jds) / ( ( b - a ) / ( npy - rem ) )     ;
57   }
58   else {
59     Py = ( a / ( (jdim / npy) + 1 ) ) + (b-a-jds) / ( ( b - a ) / ( npy - rem ) ) +
60                                         (j-b-jds) / ( ( jdim / npy ) + 1 )  ;
61   }
62 
63   *Px_p = Px ;
64   *Py_p = Py ;
65 }
66 
67 #if 0
68 main()
69 {
70   int ips[100], ipe[100] ;
71   int jps[100], jpe[100] ;
72   int shw, i , j , ids, ide, jds, jde, npx, npy ;  /* inputs */
73   int Px, Py, P ;                             /* output */
74   printf("i, j, ids, ide, jds, jde, npx, npy\n") ;
75   scanf("%d %d %d %d %d %d %d %d",&i, &j, &ids,&ide,&jds,&jde,&npx,&npy ) ;
76   shw =0 ;
77   for ( i = 0 ; i < 100 ; i++ ) { ips[i] = 9999999 ; ipe[i] = -99999999 ; }
78   for ( i = 0 ; i < 100 ; i++ ) { jps[i] = 9999999 ; jpe[i] = -99999999 ; }
79 #if 1
80   for ( j = jds-shw ; j <= jde+shw ; j++ )
81   {
82   for ( i = ids-shw ; i <= ide+shw ; i++ )
83   {
84 #endif
85   TASK_FOR_POINT ( &i , &j ,
86                    &ids, &ide, &jds, &jde , &npx , &npy ,
87                    &Px, &Py ) ;
88 /*  printf("%3d",P) ; */
89 #if 1
90   }
91 /*  printf("\n") ; */
92   }
93 for ( i = 0 ; i < npx*npy ; i++ ) {
94   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 ) ;
95 }
96 #endif
97 }
98 #endif
99