wrfv2_fire/external/RSL/RSL/cd3.c

   1 /***********************************************************************
   2
   3                               COPYRIGHT
   4
   5      The following is a notice of limited availability of the code and
   6      Government license and disclaimer which must be included in the
   7      prologue of the code and in all source listings of the code.
   8
   9      Copyright notice
  10        (c) 1977  University of Chicago
  11
  12      Permission is hereby granted to use, reproduce, prepare
  13      derivative works, and to redistribute to others at no charge.  If
  14      you distribute a copy or copies of the Software, or you modify a
  15      copy or copies of the Software or any portion of it, thus forming
  16      a work based on the Software and make and/or distribute copies of
  17      such work, you must meet the following conditions:
  18
  19           a) If you make a copy of the Software (modified or verbatim)
  20              it must include the copyright notice and Government
  21              license and disclaimer.
  22
  23           b) You must cause the modified Software to carry prominent
  24              notices stating that you changed specified portions of
  25              the Software.
  26
  27      This software was authored by:
  28
  29      Argonne National Laboratory
  30      J. Michalakes: (630) 252-6646; email: michalak@mcs.anl.gov
  31      Mathematics and Computer Science Division
  32      Argonne National Laboratory, Argonne, IL  60439
  33
  34      ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES
  35      OF ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT,
  36      AND OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A
  37      CONTRACT WITH THE DEPARTMENT OF ENERGY.
  38
  39                       GOVERNMENT LICENSE AND DISCLAIMER
  40
  41      This computer code material was prepared, in part, as an account
  42      of work sponsored by an agency of the United States Government.
  43      The Government is granted for itself and others acting on its
  44      behalf a paid-up, nonexclusive, irrevocable worldwide license in
  45      this data to reproduce, prepare derivative works, distribute
  46      copies to the public, perform publicly and display publicly, and
  47      to permit others to do so.  NEITHER THE UNITED STATES GOVERNMENT
  48      NOR ANY AGENCY THEREOF, NOR THE UNIVERSITY OF CHICAGO, NOR ANY OF
  49      THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
  50      ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
  51      COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS,
  52      PRODUCT, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD
  53      NOT INFRINGE PRIVATELY OWNED RIGHTS.
  54
  55 ***************************************************************************/
  56
  57 #include <stdio.h>
  58 #include <stdlib.h>
  59 #include <math.h>
  60
  61
  62 #define INDEX_2(A,B,NB)       ( (B) + (A)*(NB) )
  63
  64 mark_line( wrk, m, n, x1, y1, x2, y2, val )
  65   int wrk[] ;
  66   int m, n, x1, y1, x2, y2, val ;
  67 {
  68   int x, y, i, j, k ;
  69   int yz, yz2, dex ;
  70
  71   double SLOPE, X, Y, DX, DY, dY, X1, X2, Y1, Y2 ;
  72
  73   X1 = x1 ;
  74   X2 = x2 ;
  75   DX = X2 - X1 ;
  76   Y1 = y1 ;
  77   Y2 = y2 ;
  78   DY = Y2 - Y1 ;
  79   if ( DX == 0.0 )
  80   {
  81     if ( y2 >= y1 )
  82       for ( i = y1 ; i <= y2 ; i++ )
  83         wrk[ INDEX_2( x1, i, m )] = val ;
  84     else
  85       for ( i = y2 ; i <= y1 ; i++ )
  86         wrk[ INDEX_2( x1, i, m )] = val ;
  87     return ;
  88   }
  89   else
  90   {
  91     SLOPE = DY/DX ;
  92     dY = SLOPE * .5 ;
  93   }
  94
  95   if ( x2 >= X1 )
  96   {
  97     Y = y1 + .5 ;
  98     for ( X = x1+.5 ; X < x2+.5 ; X = X+.5 )
  99     {
 100       j = X + .25 ;
 101       if ( dY >= 0.0 )
 102         for ( i = Y ; i <= (k=Y+dY) ; i++ )     /* k business converts to int */
 103           wrk[ INDEX_2( j, i, m )] = val ;
 104       else
 105         for ( i = Y ; i >= (k=Y+dY) ; i-- )
 106           wrk[ INDEX_2( j, i, m )] = val ;
 107       Y = Y + dY ;
 108     }
 109     wrk[ INDEX_2( x2, y2, m )] = val ;
 110   }
 111   else
 112   {
 113     Y = y2 + .5 ;
 114     for ( X = x2+.5 ; X < x1+.5 ; X = X+.5 )
 115     {
 116       j = X + .25 ;
 117       if ( dY >= 0.0 )
 118         for ( i = Y ; i <= (k=Y+dY) ; i++ )     /* k business converts to int */
 119           wrk[ INDEX_2( j, i, m )] = val ;
 120       else
 121         for ( i = Y ; i >= (k=Y+dY) ; i-- )
 122           wrk[ INDEX_2( j, i, m )] = val ;
 123       Y = Y + dY ;
 124     }
 125     wrk[ INDEX_2( x1, y1, m )] = val ;
 126   }
 127 }
 128
 129 fill_region( wrk, m, n, v, v2 )
 130   int wrk[], m, n, v, v2 ;
 131 {
 132   int x, y ;
 133
 134   for ( y = 0 ; y < m ; y++ )
 135   {
 136     flood( 0, y, v, v2, wrk, m, n ) ;
 137     flood( n-1, y, v, v2, wrk, m, n ) ;
 138   }
 139   for ( x = 0 ; x < n ; x++ )
 140   {
 141     flood( x, 0, v, v2, wrk, m, n ) ;
 142     flood( x, m-1, v, v2, wrk, m, n ) ;
 143   }
 144 }
 145
 146 decomp_region_1( wrk, m, n, py, px )
 147   int wrk[], m, n, py, px ;
 148 {
 149   int *wk ;
 150   int x, y, ncells, nprocs, n_p, n_py, n_px, i, pid, p ;
 151
 152   wk = (int*)malloc(m*n*sizeof(int)) ;
 153
 154   nprocs = px * py ;
 155
 156   for ( x = 0 ; x < m*n ; x++ )
 157   {
 158     wk[x] = -1 ;
 159     if ( wrk[x] != 0 ) ncells++ ;
 160   }
 161
 162   n_p = ncells / nprocs ;
 163
 164   /* divide over py in m dimension first */
 165   pid = -1 ;
 166   i = 0 ;
 167   n_py = ncells / py ;
 168   for ( y = 0 ; y < m ; y++ )
 169   {
 170     for ( x = 0 ; x < n ; x++ )
 171     {
 172       if ( wrk[INDEX_2(x,y,m)] != 0 )   /* only do cells in partition */
 173       {
 174         if ( i % n_py == 0 ) pid++ ;
 175         i++ ;
 176         if ( pid > px-1 ) pid = px-1 ;
 177         wk[INDEX_2(x,y,m)] = pid ;
 178       }
 179     }
 180   }
 181
 182   /* now divide over px in n dimension */
 183   n_px = n_py / px ;
 184   for ( p = 0 ; p < py ; p++ )
 185   {
 186     pid = -1 ;
 187     i = 0 ;
 188     for ( x = 0 ; x < n ; x++ )
 189     {
 190       for ( y = 0 ; y < m ; y++ )
 191       {
 192         if ( wk[INDEX_2(x,y,m)] == p )
 193         {
 194           if ( i % n_px == 0 ) pid++ ;
 195           i++ ;
 196           if ( pid > py-1 ) pid = py-1 ;
 197           wk[INDEX_2(x,y,m)] = pid*10000 + p ;
 198         }
 199       }
 200     }
 201   }
 202
 203   for ( x = 0 ; x < n ; x++ )
 204     for ( y = 0 ; y < m ; y++ )
 205   {
 206     if (( p = wk[ INDEX_2( x, y, m ) ] ) != -1 )
 207     {
 208       n_py = p % 10000 ;
 209       n_px = p / 10000 ;
 210       wrk[INDEX_2( x, y, m )] = n_py*px + n_px ;
 211     }
 212     else
 213     {
 214       wrk[INDEX_2( x, y, m )] = wk[INDEX_2( x, y, m )] ;
 215     }
 216   }
 217
 218   free(wk) ;
 219 }
 220
 221 decomp_region_2( wrk, m, n, py, px )
 222   int wrk[], m, n, py, px ;
 223 {
 224   int *wk ;
 225   int x, y, ncells, nprocs, n_p, n_py, n_px, i, pid, p ;
 226
 227
 228   printf("decomp_region_2( wrk, %d, %d, %d, %d )\n",m, n, py, px ) ;
 229   wk = (int*)malloc(m*n*sizeof(int)) ;
 230
 231   nprocs = px * py ;
 232
 233   for ( x = 0 ; x < m*n ; x++ )
 234   {
 235     wk[x] = -1 ;
 236     if ( wrk[x] != 0 ) ncells++ ;
 237   }
 238
 239   n_p = ncells / nprocs ;
 240
 241   dc2( 0, nprocs, wk, wrk, m, n, py, px, n_p ) ;
 242
 243   for ( x = 0 ; x < n ; x++ )
 244     for ( y = 0 ; y < m ; y++ )
 245   {
 246     if (( p = wk[ INDEX_2( x, y, m ) ] ) != -1 )
 247     {
 248       n_py = p % 10000 ;
 249       n_px = p / 10000 ;
 250       wrk[INDEX_2( x, y, m )] = n_py*px + n_px ;
 251     }
 252     else
 253     {
 254       wrk[INDEX_2( x, y, m )] = wk[INDEX_2( x, y, m )] ;
 255     }
 256   }
 257
 258   free(wk) ;
 259 }
 260
 261 dc2( p, nprocs, wk, wrk, m, n, py, px, n_p )
 262   int p, nprocs, wk[], wrk[], m, n, py, px, n_p ;
 263 {
 264   int x, y, i, v, flg, reach, oldi ;
 265
 266   if ( p >= nprocs ) return ;
 267
 268   printf("dc2(%d, %d, wk, wk, %d, %d, %d, %d, %d)\n",p, nprocs,m, n, py, px, n_p );
 269
 270
 271   for ( x = 0 ; x < n ; x++ )
 272   {
 273     flg = 1 ;
 274     for ( y = 0 ; y < m && flg ; y++ )
 275     {
 276       if ( wrk[INDEX_2(x,y,m)] != 0 )
 277       {
 278         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 279         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 280       }
 281     }
 282   }
 283
 284   for ( y = 0 ; y < m ; y++ )
 285   {
 286     flg = 1 ;
 287     for ( x = n-1 ; x >= 0 && flg ; x-- )
 288     {
 289       if ( wrk[INDEX_2(x,y,m)] != 0 )
 290       {
 291         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 292         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 293       }
 294     }
 295   }
 296
 297   for ( x = n-1 ; x >=0 ; x-- )
 298   {
 299     flg = 1 ;
 300     for ( y = m-1 ; y >=0 && flg ; y-- )
 301     {
 302       if ( wrk[INDEX_2(x,y,m)] != 0 )
 303       {
 304         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 305         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 306       }
 307     }
 308   }
 309
 310   for ( y = m-1 ; y >= 0 ; y-- )
 311   {
 312     flg = 1 ;
 313     for ( x = 0 ; x < n && flg ; x++ )
 314     {
 315       if ( wrk[INDEX_2(x,y,m)] != 0 )
 316       {
 317         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 318         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 319       }
 320     }
 321   }
 322
 323   for ( x = 0 ; x < n ; x++ )
 324   {
 325     for ( y = 0 ; y < m ; y++ )
 326     {
 327       v = wk[INDEX_2(x,y,m)] ;
 328       if ( v == -1 && wrk[INDEX_2(x,y,m)] != 0 )
 329       {
 330         goto breakout ;
 331       }
 332     }
 333   }
 334
 335 breakout:
 336   if ( x == n ) return ;                /* done, none found */
 337
 338   printf("dc2: p %d v %d x %d y %d  wrk %d\n", p,v,x,y,wrk[INDEX_2(x,y,m)]) ;
 339
 340   reach = 0 ;
 341   i = 0 ;
 342   /* start acreting outward until we're stopped dead or we get enough */
 343   while ( (i < n_p) )
 344   {
 345     reach++ ;
 346     oldi = i ;
 347     acrete( wk, wrk, p, &i, n_p, reach, x, y, m, n ) ;
 348     if ( i == oldi ) break ;
 349   }
 350
 351   printf("\n") ;
 352   print_region( wk, m, n ) ;
 353
 354   dc2( p+1, nprocs, wk, wrk, m, n, py, px, n_p ) ;
 355 }
 356
 357 #define T(X,Y) \
 358    {                                                                    \
 359      if( (X) >= 0 && (X) < n && (Y) >= 0 && (Y) < m && *i < n_p )       \
 360        if (wk[INDEX_2((X),(Y),m)]==-1&&wrk[INDEX_2((X),(Y),m)]!=0)      \
 361        {                                                                \
 362            setone = 1 ;                                                 \
 363            wk[INDEX_2((X),(Y),m)] = p ;                                 \
 364            *i = *i + 1 ;                                                \
 365        }                                                                \
 366    }
 367
 368 acrete( wk, wrk, p, i, n_p, reach, x, y, m, n )
 369   int wk[], wrk[], p, *i, n_p ;
 370 {
 371   int setone ;
 372
 373   if ( reach == 0 ) return ;
 374   if ( *i >= n_p ) return ;
 375   if ( !(x >= 0 && x < n && y >= 0 && y < m) ) return ;
 376   if ( wrk[INDEX_2(x,y,m)]==0 ) return ;
 377   if ( wk[INDEX_2(x,y,m)]!=-1 && wk[INDEX_2(x,y,m)]!= p) return ;
 378
 379   /*
 380   if ( p == 5 && reach == 1 )
 381   printf("acrete(wk,wrk, p %d, *i %d, n_p %d, reach %d, x %d, y %d, m %d, n %d)\n",p, *i, n_p, reach, x, y, m, n ) ;
 382   */
 383
 384   T( x   , y   ) ;
 385   T( x-1 , y+1 ) ;
 386   T( x   , y+1 ) ;
 387   T( x+1 , y+1 ) ;
 388   T( x-1 , y   ) ;
 389   T( x+1 , y   ) ;
 390   T( x-1 , y-1 ) ;
 391   T( x   , y-1 ) ;
 392   T( x+1 , y-1 ) ;
 393
 394   acrete( wk, wrk, p, i, n_p, reach-1, x-1 , y-1 , m, n ) ;
 395   acrete( wk, wrk, p, i, n_p, reach-1, x-1 , y   , m, n ) ;
 396   acrete( wk, wrk, p, i, n_p, reach-1, x-1 , y+1 , m, n ) ;
 397   acrete( wk, wrk, p, i, n_p, reach-1, x   , y+1 , m, n ) ;
 398   acrete( wk, wrk, p, i, n_p, reach-1, x+1 , y+1 , m, n ) ;
 399   acrete( wk, wrk, p, i, n_p, reach-1, x+1 , y   , m, n ) ;
 400   acrete( wk, wrk, p, i, n_p, reach-1, x+1 , y-1 , m, n ) ;
 401   acrete( wk, wrk, p, i, n_p, reach-1, x   , y-1 , m, n ) ;
 402
 403 }
 404
 405
 406 flood( x, y, v, v2, wrk, m, n )
 407   int x, y, v, v2, wrk[], m, n ;
 408 {
 409   if ( x < 0 || x >= n || y < 0 || y >= m )
 410     return ;
 411   if ( wrk[INDEX_2(x,y,m)] == v )
 412   {
 413     wrk[INDEX_2(x,y,m)] = v2 ;
 414     flood( x+1, y  , v, v2, wrk, m, n ) ;
 415     flood( x-1, y  , v, v2, wrk, m, n ) ;
 416     flood( x  , y+1, v, v2, wrk, m, n ) ;
 417     flood( x  , y-1, v, v2, wrk, m, n ) ;
 418   }
 419 }
 420
 421 print_region( wrk, m, n )
 422   int wrk[], m, n ;
 423 {
 424   int i, j ;
 425   for ( i = m-1 ; i >= 0 ; i-- )
 426   {
 427     for ( j = 0 ; j < n ; j++ )
 428     {
 429       if ( wrk[ INDEX_2( j, i, m ) ] == -1 )
 430         printf("  ." ) ;
 431       else
 432         printf("%3d", wrk[ INDEX_2( j, i, m ) ] ) ;
 433     }
 434     printf("\n") ;
 435   }
 436 }
 437
 438 /****************/
 439
 440 main()
 441 {
 442   /*      n   m */
 443   int wrk[100 * 100] ;
 444   int i, j, m, n, x1, y1, x2, y2, py, px, opt ;
 445
 446   for ( i = 0 ; i < 100*100 ; i++ ) wrk[i] = 0 ;
 447
 448   scanf("%d %d", &m, &n) ;
 449   scanf("%d %d %d", &py, &px, &opt) ;
 450   scanf("%d %d", &x1, &y1) ;
 451   while ( scanf("%d %d", &x2, &y2) != EOF )
 452   {
 453     mark_line( wrk, m, n, x1, y1, x2, y2, 1 ) ;
 454     x1 = x2 ;
 455     y1 = y2 ;
 456   }
 457
 458   print_region( wrk, m, n ) ;
 459
 460   fill_region( wrk, m, n, 0, 2 ) ;
 461   for ( j = 0 ; j < n ; j++ )
 462     for ( i = 0 ; i < m ; i++ )
 463     {
 464       if ( wrk[ INDEX_2( j, i, m ) ] == 2 )
 465         wrk[ INDEX_2( j, i, m ) ] = 0 ;
 466       else
 467         wrk[ INDEX_2( j, i, m ) ] = 1 ;
 468     }
 469
 470   print_region( wrk, m, n ) ;
 471
 472   switch ( opt )
 473   {
 474     case 1 : decomp_region_1( wrk, m, n, py, px ) ;
 475              print_region( wrk, m, n ) ;
 476              break ;
 477     case 2 : decomp_region_2( wrk, m, n, py, px ) ;
 478              break ;
 479     default : break ;
 480   }
 481
 482   printf("\n") ;
 483
 484
 485 }
 486
 487