wrfv2_fire/external/RSL/RSL/cd1.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 )
  65   int wrk[] ;
  66   int m, n, x1, y1, x2, y2 ;
  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 )] = 1 ;
  84     else
  85       for ( i = y2 ; i <= y1 ; i++ )
  86         wrk[ INDEX_2( x1, i, m )] = 1 ;
  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 )] = 1 ;
 104       else
 105         for ( i = Y ; i >= (k=Y+dY) ; i-- )
 106           wrk[ INDEX_2( j, i, m )] = 1 ;
 107       Y = Y + dY ;
 108     }
 109     wrk[ INDEX_2( x2, y2, m )] = 1 ;
 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 )] = 1 ;
 120       else
 121         for ( i = Y ; i >= (k=Y+dY) ; i-- )
 122           wrk[ INDEX_2( j, i, m )] = 1 ;
 123       Y = Y + dY ;
 124     }
 125     wrk[ INDEX_2( x1, y1, m )] = 1 ;
 126   }
 127 }
 128
 129 fill_region( wrk, m, n )
 130   int wrk[], m, n ;
 131 {
 132   int x, y ;
 133
 134   for ( y = 0 ; y < m ; y++ )
 135   {
 136     flood( 0, y, 0, 2, wrk, m, n ) ;
 137     flood( n-1, y, 0, 2, wrk, m, n ) ;
 138   }
 139   for ( x = 0 ; x < n ; x++ )
 140   {
 141     flood( x, 0, 0, 2, wrk, m, n ) ;
 142     flood( x, m-1, 0, 2, wrk, m, n ) ;
 143   }
 144
 145   for ( y = 0 ; y < n ; y++ )
 146   {
 147     for ( x = 0 ; x < m ; x++ )
 148     {
 149       if      ( wrk[ INDEX_2(x,y,m) ] == 0 )
 150         wrk[ INDEX_2(x,y,m) ] = 1 ;
 151       else if ( wrk[ INDEX_2(x,y,m) ] == 2 )
 152         wrk[ INDEX_2(x,y,m) ] = 0 ;
 153     }
 154   }
 155 }
 156
 157 decomp_region_1( wrk, m, n, py, px )
 158   int wrk[], m, n, py, px ;
 159 {
 160   int *wk ;
 161   int x, y, ncells, nprocs, n_p, n_py, n_px, i, pid, p ;
 162
 163   wk = (int*)malloc(m*n*sizeof(int)) ;
 164
 165   nprocs = px * py ;
 166
 167   for ( x = 0 ; x < m*n ; x++ )
 168   {
 169     wk[x] = -1 ;
 170     if ( wrk[x] != 0 ) ncells++ ;
 171   }
 172
 173   n_p = ncells / nprocs ;
 174
 175   /* divide over py in m dimension first */
 176   pid = 0 ;
 177   i = 0 ;
 178   n_py = ncells / py ;
 179   for ( y = 0 ; y < m ; y++ )
 180   {
 181     for ( x = 0 ; x < n ; x++ )
 182     {
 183       if ( wrk[INDEX_2(x,y,m)] != 0 )   /* only do cells in partition */
 184       {
 185         i++ ;
 186         if ( i % n_py == 0 ) pid++ ;
 187         wk[INDEX_2(x,y,m)] = pid ;
 188       }
 189     }
 190   }
 191
 192   /* now divide over px in n dimension */
 193   n_px = n_py / px ;
 194   for ( p = 0 ; p < py ; p++ )
 195   {
 196     pid = 0 ;
 197     i = 0 ;
 198     for ( x = 0 ; x < n ; x++ )
 199     {
 200       for ( y = 0 ; y < m ; y++ )
 201       {
 202         if ( wk[INDEX_2(x,y,m)] == p )
 203         {
 204           i++ ;
 205           if ( i % n_px == 0 ) pid++ ;
 206           if ( pid > px-1 ) pid = px-1 ;
 207           wk[INDEX_2(x,y,m)] = pid*10000 + p ;
 208         }
 209       }
 210     }
 211   }
 212
 213   for ( x = 0 ; x < n ; x++ )
 214     for ( y = 0 ; y < m ; y++ )
 215   {
 216     if (( p = wk[ INDEX_2( x, y, m ) ] ) != -1 )
 217     {
 218       n_py = p % 10000 ;
 219       n_px = p / 10000 ;
 220       wrk[INDEX_2( x, y, m )] = n_py*px + n_px ;
 221     }
 222     else
 223     {
 224       wrk[INDEX_2( x, y, m )] = wk[INDEX_2( x, y, m )] ;
 225     }
 226   }
 227
 228   free(wk) ;
 229 }
 230
 231 decomp_region_2( wrk, m, n, py, px )
 232   int wrk[], m, n, py, px ;
 233 {
 234   int *wk ;
 235   int x, y, ncells, nprocs, n_p, n_py, n_px, i, pid, p ;
 236
 237
 238   printf("decomp_region_2( wrk, %d, %d, %d, %d )\n",m, n, py, px ) ;
 239   wk = (int*)malloc(m*n*sizeof(int)) ;
 240
 241   nprocs = px * py ;
 242
 243   for ( x = 0 ; x < m*n ; x++ )
 244   {
 245     wk[x] = -1 ;
 246     if ( wrk[x] != 0 ) ncells++ ;
 247   }
 248
 249   n_p = ncells / nprocs ;
 250
 251   dc2( 0, nprocs, wk, wrk, m, n, py, px, n_p ) ;
 252
 253   for ( x = 0 ; x < n ; x++ )
 254     for ( y = 0 ; y < m ; y++ )
 255   {
 256     if (( p = wk[ INDEX_2( x, y, m ) ] ) != -1 )
 257     {
 258       n_py = p % 10000 ;
 259       n_px = p / 10000 ;
 260       wrk[INDEX_2( x, y, m )] = n_py*px + n_px ;
 261     }
 262     else
 263     {
 264       wrk[INDEX_2( x, y, m )] = wk[INDEX_2( x, y, m )] ;
 265     }
 266   }
 267
 268   free(wk) ;
 269 }
 270
 271 dc2( p, nprocs, wk, wrk, m, n, py, px, n_p )
 272   int p, nprocs, wk[], wrk[], m, n, py, px, n_p ;
 273 {
 274   int x, y, i, v, flg, reach, oldi ;
 275
 276   if ( p >= nprocs ) return ;
 277
 278   printf("dc2(%d, %d, wk, wk, %d, %d, %d, %d, %d)\n",p, nprocs,m, n, py, px, n_p );
 279
 280
 281   for ( x = 0 ; x < n ; x++ )
 282   {
 283     flg = 1 ;
 284     for ( y = 0 ; y < m && flg ; y++ )
 285     {
 286       if ( wrk[INDEX_2(x,y,m)] != 0 )
 287       {
 288         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 289         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 290       }
 291     }
 292   }
 293
 294   for ( y = 0 ; y < m ; y++ )
 295   {
 296     flg = 1 ;
 297     for ( x = n-1 ; x >= 0 && flg ; x-- )
 298     {
 299       if ( wrk[INDEX_2(x,y,m)] != 0 )
 300       {
 301         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 302         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 303       }
 304     }
 305   }
 306
 307   for ( x = n-1 ; x >=0 ; x-- )
 308   {
 309     flg = 1 ;
 310     for ( y = m-1 ; y >=0 && flg ; y-- )
 311     {
 312       if ( wrk[INDEX_2(x,y,m)] != 0 )
 313       {
 314         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 315         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 316       }
 317     }
 318   }
 319
 320   for ( y = m-1 ; y >= 0 ; y-- )
 321   {
 322     flg = 1 ;
 323     for ( x = 0 ; x < n && flg ; x++ )
 324     {
 325       if ( wrk[INDEX_2(x,y,m)] != 0 )
 326       {
 327         if ( wk[INDEX_2(x,y,m)] == -1 ) goto breakout ;
 328         if ( wk[INDEX_2(x,y,m)] != -1 ) flg = 0 ;
 329       }
 330     }
 331   }
 332
 333   for ( x = 0 ; x < n ; x++ )
 334   {
 335     for ( y = 0 ; y < m ; y++ )
 336     {
 337       v = wk[INDEX_2(x,y,m)] ;
 338       if ( v == -1 && wrk[INDEX_2(x,y,m)] != 0 )
 339       {
 340         goto breakout ;
 341       }
 342     }
 343   }
 344
 345 breakout:
 346   if ( x == n ) return ;                /* done, none found */
 347
 348   printf("dc2: p %d v %d x %d y %d  wrk %d\n", p,v,x,y,wrk[INDEX_2(x,y,m)]) ;
 349
 350   reach = 0 ;
 351   i = 0 ;
 352   /* start acreting outward until we're stopped dead or we get enough */
 353   while ( (i < n_p) )
 354   {
 355     reach++ ;
 356     oldi = i ;
 357     acrete( wk, wrk, p, &i, n_p, reach, x, y, m, n ) ;
 358     if ( i == oldi ) break ;
 359   }
 360
 361   printf("\n") ;
 362   print_region( wk, m, n ) ;
 363
 364   dc2( p+1, nprocs, wk, wrk, m, n, py, px, n_p ) ;
 365 }
 366
 367 #define T(X,Y) \
 368    {                                                                    \
 369      if( (X) >= 0 && (X) < n && (Y) >= 0 && (Y) < m && *i < n_p )       \
 370        if (wk[INDEX_2((X),(Y),m)]==-1&&wrk[INDEX_2((X),(Y),m)]!=0)      \
 371        {                                                                \
 372            setone = 1 ;                                                 \
 373            wk[INDEX_2((X),(Y),m)] = p ;                                 \
 374            *i = *i + 1 ;                                                \
 375        }                                                                \
 376    }
 377
 378 acrete( wk, wrk, p, i, n_p, reach, x, y, m, n )
 379   int wk[], wrk[], p, *i, n_p ;
 380 {
 381   int setone ;
 382
 383   if ( reach == 0 ) return ;
 384   if ( *i >= n_p ) return ;
 385   if ( !(x >= 0 && x < n && y >= 0 && y < m) ) return ;
 386   if ( wrk[INDEX_2(x,y,m)]==0 ) return ;
 387   if ( wk[INDEX_2(x,y,m)]!=-1 && wk[INDEX_2(x,y,m)]!= p) return ;
 388
 389   /*
 390   if ( p == 5 && reach == 1 )
 391   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 ) ;
 392   */
 393
 394   T( x   , y   ) ;
 395   T( x-1 , y+1 ) ;
 396   T( x   , y+1 ) ;
 397   T( x+1 , y+1 ) ;
 398   T( x-1 , y   ) ;
 399   T( x+1 , y   ) ;
 400   T( x-1 , y-1 ) ;
 401   T( x   , y-1 ) ;
 402   T( x+1 , y-1 ) ;
 403
 404   acrete( wk, wrk, p, i, n_p, reach-1, x-1 , y-1 , m, n ) ;
 405   acrete( wk, wrk, p, i, n_p, reach-1, x-1 , y   , m, n ) ;
 406   acrete( wk, wrk, p, i, n_p, reach-1, x-1 , y+1 , m, n ) ;
 407   acrete( wk, wrk, p, i, n_p, reach-1, x   , y+1 , m, n ) ;
 408   acrete( wk, wrk, p, i, n_p, reach-1, x+1 , y+1 , m, n ) ;
 409   acrete( wk, wrk, p, i, n_p, reach-1, x+1 , y   , m, n ) ;
 410   acrete( wk, wrk, p, i, n_p, reach-1, x+1 , y-1 , m, n ) ;
 411   acrete( wk, wrk, p, i, n_p, reach-1, x   , y-1 , m, n ) ;
 412
 413 }
 414
 415
 416 flood( x, y, v, v2, wrk, m, n )
 417   int x, y, v, v2, wrk[], m, n ;
 418 {
 419   if ( x < 0 || x >= n || y < 0 || y >= m )
 420     return ;
 421   if ( wrk[INDEX_2(x,y,m)] == v )
 422   {
 423     wrk[INDEX_2(x,y,m)] = v2 ;
 424     flood( x+1, y  , v, v2, wrk, m, n ) ;
 425     flood( x-1, y  , v, v2, wrk, m, n ) ;
 426     flood( x  , y+1, v, v2, wrk, m, n ) ;
 427     flood( x  , y-1, v, v2, wrk, m, n ) ;
 428   }
 429 }
 430
 431 print_region( wrk, m, n )
 432   int wrk[], m, n ;
 433 {
 434   int i, j ;
 435   for ( i = m-1 ; i >= 0 ; i-- )
 436   {
 437     for ( j = 0 ; j < n ; j++ )
 438     {
 439       if ( wrk[ INDEX_2( j, i, m ) ] == -1 )
 440         printf("  ." ) ;
 441       else
 442         printf("%3d", wrk[ INDEX_2( j, i, m ) ] ) ;
 443     }
 444     printf("\n") ;
 445   }
 446 }
 447
 448 /****************/
 449
 450 main()
 451 {
 452   /*      n   m */
 453   int wrk[100 * 100] ;
 454   int i, j, m, n, x1, y1, x2, y2, py, px, opt ;
 455
 456   for ( i = 0 ; i < 100*100 ; i++ ) wrk[i] = 0 ;
 457
 458   scanf("%d %d", &m, &n) ;
 459   scanf("%d %d %d", &py, &px, &opt) ;
 460   scanf("%d %d", &x1, &y1) ;
 461   while ( scanf("%d %d", &x2, &y2) != EOF )
 462   {
 463     mark_line( wrk, m, n, x1, y1, x2, y2 ) ;
 464     x1 = x2 ;
 465     y1 = y2 ;
 466   }
 467
 468   /* print_region( wrk, m, n ) ; */
 469
 470   fill_region( wrk, m, n ) ;
 471
 472   /* print_region( wrk, m, n ) ; */
 473
 474   switch ( opt )
 475   {
 476     case 1 : decomp_region_1( wrk, m, n, py, px ) ;
 477              print_region( wrk, m, n ) ;
 478              break ;
 479     case 2 : decomp_region_2( wrk, m, n, py, px ) ;
 480              break ;
 481     default : break ;
 482   }
 483
 484   printf("\n") ;
 485
 486
 487 }
 488
 489