wrfv2_fire/external/RSL/RSL/rsl_hemiforce.c

   1 #define KLUDGE_20000821
   2
   3 /***********************************************************************
   4
   5                               COPYRIGHT
   6
   7      The following is a notice of limited availability of the code and
   8      Government license and disclaimer which must be included in the
   9      prologue of the code and in all source listings of the code.
  10
  11      Copyright notice
  12        (c) 1977  University of Chicago
  13
  14      Permission is hereby granted to use, reproduce, prepare
  15      derivative works, and to redistribute to others at no charge.  If
  16      you distribute a copy or copies of the Software, or you modify a
  17      copy or copies of the Software or any portion of it, thus forming
  18      a work based on the Software and make and/or distribute copies of
  19      such work, you must meet the following conditions:
  20
  21           a) If you make a copy of the Software (modified or verbatim)
  22              it must include the copyright notice and Government
  23              license and disclaimer.
  24
  25           b) You must cause the modified Software to carry prominent
  26              notices stating that you changed specified portions of
  27              the Software.
  28
  29      This software was authored by:
  30
  31      Argonne National Laboratory
  32      J. Michalakes: (630) 252-6646; email: michalak@mcs.anl.gov
  33      Mathematics and Computer Science Division
  34      Argonne National Laboratory, Argonne, IL  60439
  35
  36      ARGONNE NATIONAL LABORATORY (ANL), WITH FACILITIES IN THE STATES
  37      OF ILLINOIS AND IDAHO, IS OWNED BY THE UNITED STATES GOVERNMENT,
  38      AND OPERATED BY THE UNIVERSITY OF CHICAGO UNDER PROVISION OF A
  39      CONTRACT WITH THE DEPARTMENT OF ENERGY.
  40
  41                       GOVERNMENT LICENSE AND DISCLAIMER
  42
  43      This computer code material was prepared, in part, as an account
  44      of work sponsored by an agency of the United States Government.
  45      The Government is granted for itself and others acting on its
  46      behalf a paid-up, nonexclusive, irrevocable worldwide license in
  47      this data to reproduce, prepare derivative works, distribute
  48      copies to the public, perform publicly and display publicly, and
  49      to permit others to do so.  NEITHER THE UNITED STATES GOVERNMENT
  50      NOR ANY AGENCY THEREOF, NOR THE UNIVERSITY OF CHICAGO, NOR ANY OF
  51      THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR
  52      ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY FOR THE ACCURACY,
  53      COMPLETENESS, OR USEFULNESS OF ANY INFORMATION, APPARATUS,
  54      PRODUCT, OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS USE WOULD
  55      NOT INFRINGE PRIVATELY OWNED RIGHTS.
  56
  57 ***************************************************************************/
  58
  59 #include <stdio.h>
  60 #include <stdlib.h>
  61 #include "rsl.h"
  62
  63 /*@
  64   RSL_TO_OH_INFO -- Get the next cell in a packing sequence for forcing.
  65
  66   Notes:
  67
  68   See also:
  69
  70 @*/
  71
  72 static rsl_domain_info_t *s_tinfo, *s_oinfo ;
  73 static int s_oig, s_ojg ;
  74 static int s_p, s_t, s_o ;
  75 static int s_msize ;
  76 static struct rsl_hemi_rec * s_q, * s_q1 ;
  77 static int s_nlen_o ;
  78 static int s_mlen_o ;
  79 static rsl_point_t *s_tdomain, *s_odomain ;
  80 static char * s_pointbuf = NULL ;
  81
  82 RSL_TO_OH_INFO ( t_p, o_p, msize_p, seed_p,
  83                  oig_p, ojg_p, retval_p )
  84   int_p
  85      t_p            /* (I) RSL domain descriptor of this hemi. */
  86     ,o_p            /* (I) RSL domain descriptor of other hemi. */
  87     ,msize_p        /* (I) Message size in bytes. */
  88     ,seed_p         /* (I) =1 start the traversal; =0 (zero) continue traversal */
  89     ,oig_p          /* (O) Global M index of other domain point. */
  90     ,ojg_p          /* (O) Global N index of other domain point. */
  91     ,retval_p ;     /* (O) =1 if a valid point returned; =0 (zero) otherwise. */
  92 {
  93   int kiddex ;
  94   int P ;
  95   rsl_hemi_rec_t * q ;
  96 #ifdef KLUDGE_20000821
  97   rsl_hemi_rec_t * qnuke ;
  98   rsl_hemi_rec_t * prev ;
  99 #endif
 100   int p, p1 ;
 101   int globalhemiPlist[RSL_MAXPROC][RSL_MAXPROC], work[RSL_MAXPROC][RSL_MAXPROC] ;
 102
 103 #ifndef STUBS
 104   s_msize = *msize_p ;
 105   s_t = *t_p ;
 106   s_o = *o_p ;
 107   RSL_TEST_ERR( s_t < 0 || s_t > RSL_MAXDOMAINS,
 108     "rsl_ready_bcast: bad 'this hemi' descriptor" ) ;
 109   RSL_TEST_ERR( s_o < 0 || s_o > RSL_MAXDOMAINS,
 110     "rsl_ready_bcast: bad 'other hemi' descriptor" ) ;
 111   RSL_TEST_ERR( s_t == s_o,
 112     "rsl_ready_bcast: hemispere cannot force itself" ) ;
 113
 114   s_tinfo = &( domain_info[s_t]) ;
 115   s_oinfo = &( domain_info[s_o]) ;
 116   s_mlen_o = s_oinfo->len_m ;
 117   s_nlen_o = s_oinfo->len_n ;
 118   s_odomain = s_oinfo->domain ;
 119
 120   if ( ! s_tinfo->other_hemi_proclist_built )
 121   {
 122     if ( *seed_p )
 123     {
 124       s_oig = 0 ; s_ojg = 0 ;
 125       for ( p = 0 ; p < RSL_MAXPROC ; p++ )
 126       {
 127 #ifdef KLUDGE_20000821
 128         for ( q = s_tinfo->other_hemi_procbufs[p], prev = NULL ; q ; )
 129         {
 130            if ( q->data ) RSL_FREE( q->data ) ;
 131            qnuke = q ;
 132            q = q->next ;
 133            RSL_FREE( qnuke ) ;
 134         }
 135 #endif
 136         s_tinfo->other_hemi_procbufs[p] = NULL ;
 137         s_tinfo->hemi_sendPlist[p] = 0 ;
 138         for ( p1 = 0 ; p1 < RSL_MAXPROC ; p1++ )
 139         {
 140           globalhemiPlist[p][p1] = 0 ;
 141         }
 142       }
 143     }
 144     else
 145     {
 146       s_oig++ ;
 147       if ( s_oig >= s_oinfo->len_m )
 148       {
 149         s_oig = 0 ;
 150         s_ojg++ ;
 151         if ( s_ojg >= s_oinfo->len_n )
 152         {
 153           *retval_p = 0 ;
 154 #ifndef KLUDGE_20000821
 155           s_tinfo->other_hemi_proclist_built = 1 ; /* FIX 20000818 JM */
 156  #endif
 157
 158 /* collapse the list and keep only entries that have data associated */
 159 /* also fill entries for processors I must send to, indicating the number
 160    of columns that go to each processors, or zero for processors I don't
 161    send to */
 162           for ( p = 0 ; p < RSL_MAXPROC ; p++ )
 163           {
 164             rsl_hemi_rec_t * prev ;
 165             for ( q = s_tinfo->other_hemi_procbufs[p], prev = NULL ; q ; )
 166             {
 167               if ( q->data == NULL )
 168               {
 169                 if ( prev == NULL )
 170                 {
 171                   s_tinfo->other_hemi_procbufs[p] = q->next ;
 172                   RSL_FREE(q) ;
 173                   q = s_tinfo->other_hemi_procbufs[p] ;
 174                 }
 175                 else if ( prev->next == q )
 176                 {
 177                   prev->next = q->next ;
 178                   RSL_FREE(q) ;
 179                   q = prev->next ;
 180                 }
 181                 else
 182                   RSL_TEST_ERR(1,"internal error") ;
 183               }
 184               else
 185               {
 186                 s_tinfo->hemi_sendPlist[p]++ ;
 187                 prev = q ;
 188                 q = q->next ;
 189               }
 190             }
 191           }
 192
 193           /* mpi all to all communication to share matrix of senders/receivers */
 194           MPI_Gather( s_tinfo->hemi_sendPlist, RSL_MAXPROC, MPI_INT,
 195                       globalhemiPlist, RSL_MAXPROC, MPI_INT,
 196                       0, rsl_mpi_communicator ) ;
 197           /* transpose */
 198           for ( p = 0 ; p < RSL_MAXPROC ; p++ )
 199             for ( p1 = 0 ; p1 < RSL_MAXPROC ; p1++ )
 200               work[p][p1] = globalhemiPlist[p1][p] ;
 201           for ( p = 0 ; p < RSL_MAXPROC ; p++ )
 202             for ( p1 = 0 ; p1 < RSL_MAXPROC ; p1++ )
 203               globalhemiPlist[p][p1] = work[p][p1] ;
 204           MPI_Scatter( globalhemiPlist, RSL_MAXPROC, MPI_INT,
 205                        s_tinfo->hemi_recvPlist, RSL_MAXPROC, MPI_INT,
 206                        0, rsl_mpi_communicator ) ;
 207
 208           return ;       /* EARLY RETURN */
 209         }
 210       }
 211     }
 212     kiddex = INDEX_2(s_ojg,s_oig,s_mlen_o) ;
 213     P = s_odomain[ kiddex ].P ;
 214     if ( s_tinfo->other_hemi_procbufs[P] == NULL )
 215     {
 216       q = RSL_MALLOC( rsl_hemi_rec_t, 1 ) ;
 217     }
 218     else
 219     {
 220       q = RSL_MALLOC( rsl_hemi_rec_t, 1 ) ;
 221       q->next = s_tinfo->other_hemi_procbufs[P] ;
 222     }
 223     q->oig = s_oig ;
 224     q->ojg = s_ojg ;
 225     q->data = NULL ;
 226     s_tinfo->other_hemi_procbufs[P] = q ;
 227     s_q1 = q ;
 228   }
 229   else
 230   {
 231     int * x ;
 232     if ( *seed_p )
 233     {
 234       s_p = -1 ;
 235       s_q = NULL ;
 236     }
 237     if ( s_q == NULL )
 238     {
 239       s_p++ ;
 240       while ( s_tinfo->other_hemi_procbufs[s_p] == NULL ) s_p++ ;
 241       if ( s_p >= rsl_nproc_all )
 242       {
 243         *retval_p = 0 ;
 244         return ;       /* EARLY RETURN */
 245       }
 246       s_q = s_tinfo->other_hemi_procbufs[s_p] ;
 247     }
 248     s_oig = s_q->oig ;
 249     s_ojg = s_q->ojg ;
 250     s_q1 = s_q ;
 251     s_q = s_q->next ;
 252   }
 253   *oig_p = s_oig + 1 ;  /* C to Fortran */
 254   *ojg_p = s_ojg + 1 ;  /* C to Fortran */
 255   *retval_p = 1 ;
 256 #else
 257    RSL_TEST_ERR( 1, "RSL_TO_OH_INFO STUBBED" ) ;
 258 #endif
 259   return ;
 260 }
 261
 262 /*@
 263   RSL_TO_OH_MSG -- Pack force data into a message for a nest point.
 264
 265   Notes:
 266   See also:
 267
 268 @*/
 269
 270 RSL_TO_OH_MSG ( nbuf_p, buf )
 271   int_p
 272     nbuf_p ;     /* (I) Number of bytes to be packed. */
 273   char *
 274     buf ;        /* (I) Buffer containing the data to be packed. */
 275 {
 276   int kiddex ;
 277   int nbuf ;
 278   int P ;
 279
 280   RSL_TEST_ERR(buf==NULL,"2nd argument is NULL.  Field allocated?") ;
 281   nbuf = *nbuf_p ;
 282   if ( s_q1->data == NULL )
 283   {
 284     s_q1->data = RSL_MALLOC( char, s_msize ) ;
 285     s_q1->curs = 0 ;
 286   }
 287   if ( s_q1->curs+nbuf >= s_msize )
 288   {
 289     sprintf(mess,"RSL_TO_OH_MSG: store of %d bytes would overflow %d sized buffer.\n",nbuf,s_msize ) ;
 290     RSL_TEST_ERR(1,mess) ;
 291   }
 292   bcopy( buf, &(s_q1->data[s_q1->curs]), nbuf ) ;
 293   s_q1->curs += nbuf ;
 294 }
 295
 296 /*@
 297   RSL_FORCE_HEMI -- Convey forcing data from this hemi to other hemi
 298
 299   Notes:
 300
 301   See also:
 302 @*/
 303
 304 RSL_FORCE_HEMI ()
 305 {
 306   int P ;
 307   int msglen, mdest, mtag ;
 308   int ii ;
 309   int ig, jg ;
 310   char * recvbuf, * sendbuf ;
 311   rsl_hemi_rec_t * q ;
 312
 313   /* post receives */
 314
 315   for ( P = 0 ; P < rsl_nproc_all ; P++ )
 316   {
 317     if ( s_tinfo->hemi_recvPlist[P] > 0 )
 318     {
 319       msglen = s_msize * s_tinfo->hemi_recvPlist[P] + 3*sizeof(int) ;
 320       recvbuf = buffer_for_proc( P, msglen,  RSL_RECVBUF ) ;
 321       mtag = MTYPE_FROMTO( MSG_FROM_PARENT, P, rsl_myproc ) ;
 322 #ifdef DEBUGGAL
 323 fprintf(stderr,"Posting receive on tag %d\n",mtag ) ;
 324 #endif
 325       RSL_RECVBEGIN( recvbuf, msglen, mtag ) ;
 326       s_tinfo->hemi_recv_tags[P] = mtag ; /* store tag */
 327     }
 328   }
 329
 330   /* do sends */
 331
 332   for ( P = 0 ; P < rsl_nproc_all ; P++ )
 333   {
 334 #ifdef DEBUGGAL
 335 fprintf(stderr,"s_tinfo->hemi_sendPlist[P] %d\n",s_tinfo->hemi_sendPlist[P]) ;
 336 #endif
 337     if ( s_tinfo->hemi_sendPlist[P] > 0 )
 338     {
 339       int curs ;
 340       int endofdata ;
 341
 342           /* oig,ojg,nbytes,buffer * # of points + end of data */
 343       msglen = (3*sizeof(int)+s_msize)*s_tinfo->hemi_sendPlist[P]+1*sizeof(int) ;
 344       sendbuf = buffer_for_proc( P, msglen,  RSL_SENDBUF ) ;
 345       curs = 0 ;
 346       for ( q = s_tinfo->other_hemi_procbufs[P] ; q ; q = q->next )
 347       {
 348 #ifdef DEBUGGAL
 349 {
 350 int *dp ;
 351 dp = (int *) q->data ;
 352 fprintf(stderr,"> curs %d, msglen %d msize %d (%d %d) data %d\n", curs, msglen, s_msize, q->oig, q->ojg, *dp ) ;
 353 }
 354 #endif
 355         bcopy( &(q->oig), &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
 356         bcopy( &(q->ojg), &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
 357         bcopy( &(q->curs), &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
 358         bcopy( q->data, &(sendbuf[curs]), q->curs) ; curs += q->curs ;
 359       }
 360       endofdata = RSL_INVALID ;
 361       bcopy( &endofdata, &(sendbuf[curs]), sizeof(int)) ; curs += sizeof(int) ;
 362       mtag = MTYPE_FROMTO( MSG_FROM_PARENT, rsl_myproc, P ) ;
 363 #ifdef DEBUGGAL
 364 fprintf(stderr,"sending sendbuf to %d, curs = %d\n",P,curs) ;
 365 #endif
 366       RSL_SEND( sendbuf, curs, mtag, P ) ;
 367     }
 368   }
 369 }
 370
 371 /*@
 372   RSL_FROM_TH_INFO -- Get the next cell in a unpacking sequence for forcing.
 373
 374   Notes:
 375
 376   See also:
 377 @*/
 378
 379 static int s_endofdata, s_remaining, s_ndata, s_curs ;
 380 static char * s_recvbuf ;
 381
 382 RSL_FROM_TH_INFO ( seed_p, oig_p, ojg_p, retval_p )
 383   int_p
 384     seed_p      /* (I) =1 if first call; =0 otherwise */
 385    ,oig_p        /* (O) Global index in M dimension of nest. */
 386    ,ojg_p        /* (O) Global index in N dimension of nest. */
 387    ,retval_p ;  /* (O) Return value; =1 valid point, =0 done. */
 388 {
 389   int mtag ;
 390
 391 #ifdef DEBUGGAL
 392 fprintf(stderr,"RSL_FROM_TH_INFO seed = %d, s_endofdata %d\n",*seed_p,s_endofdata) ;
 393 #endif
 394
 395   if ( *seed_p == 1 )
 396   {
 397     if ( s_pointbuf != NULL ) RSL_FREE(s_pointbuf) ;
 398     s_pointbuf = RSL_MALLOC( char, 2*s_msize ) ;  /* 2 times for safety */
 399     s_p = 0 ;
 400     s_endofdata = 1 ;
 401   }
 402
 403 nextproc:
 404   if ( s_endofdata )
 405   {
 406     while ( s_tinfo->hemi_recvPlist[s_p] <= 0 ) s_p++ ;
 407     if ( s_p >= rsl_nproc_all )
 408     {
 409       *retval_p = 0 ;
 410 #ifdef DEBUGGAL
 411 fprintf(stderr,"EARLY RETURN retval = 0\n") ;
 412 #endif
 413       return ;          /* EARLY RETURN */
 414     }
 415     mtag = s_tinfo->hemi_recv_tags[s_p] ;
 416 #ifdef DEBUGGAL
 417 fprintf(stderr,"Waiting for receive on tag %d\n",mtag ) ;
 418 #endif
 419     RSL_RECVEND ( mtag ) ;
 420 #ifdef DEBUGGAL
 421 fprintf(stderr,"got receive\n") ;
 422 #endif
 423     s_recvbuf = buffer_for_proc( s_p, 0, RSL_RECVBUF ) ;
 424     s_p++ ;
 425     s_curs = 0 ;
 426     s_endofdata = 0 ;
 427   }
 428
 429 #ifdef DEBUGGAL
 430 fprintf(stderr,"before bcopy  %d, s_recvbuf %08x\n",s_curs, s_recvbuf) ;
 431 #endif
 432
 433   bcopy ( &(s_recvbuf[s_curs]), oig_p, sizeof(int) ) ; s_curs += sizeof(int) ;
 434   if ( *oig_p == RSL_INVALID )
 435   {
 436 #ifdef DEBUGGAL
 437 fprintf(stderr,"hit end of data for s_p %d, %d\n", s_p, *oig_p ) ;
 438 #endif
 439     s_endofdata = 1 ;
 440     goto nextproc ;
 441   }
 442   bcopy ( &(s_recvbuf[s_curs]), ojg_p, sizeof(int) ) ; s_curs += sizeof(int) ;
 443   bcopy ( &(s_recvbuf[s_curs]), &s_ndata, sizeof(int) ) ; s_curs += sizeof(int) ;
 444   bcopy ( &(s_recvbuf[s_curs]), s_pointbuf, s_ndata ) ; s_curs += s_ndata ;
 445   s_remaining = s_ndata ;
 446 #ifdef DEBUGGAL
 447 fprintf(stderr,"s_remaining = %d\n",s_remaining) ;
 448 #endif
 449
 450   (*oig_p) ++ ;
 451   (*ojg_p) ++ ;
 452 #ifdef DEBUGGAL
 453 fprintf(stderr,"RETURN oig ojg %d %d\n", *oig_p, *ojg_p ) ;
 454 #endif
 455   *retval_p = 1 ;
 456   return ;
 457 }
 458
 459 /*@
 460   RSL_FROM_TH_MSG -- Unpack feedback data into a nest point.
 461
 462   Notes:
 463
 464 @*/
 465 RSL_FROM_TH_MSG ( len_p, buf )
 466   int_p
 467     len_p ;          /* (I) Number of bytes to unpack. */
 468   char *
 469     buf ;            /* (O) Destination buffer. */
 470 {
 471   if ( *len_p <= 0 ) return ;
 472   if ( *len_p > s_remaining )
 473   {
 474     sprintf(mess,
 475 "RSL_FROM_TH_MSG:\n   Requested number of bytes (%d) exceeds %d, the number remaining for this point.\n", *len_p, s_remaining) ;
 476     RSL_TEST_WRN(1,mess) ;
 477   }
 478   bcopy( &(s_pointbuf[s_ndata-s_remaining]),
 479          buf,
 480          *len_p ) ;
 481
 482   s_remaining -= *len_p ;
 483 }
 484
 485 /* retval =1 if point is local, =0 otherwise */
 486 RSL_POINT_ON_PROC ( d_p, ig_p, jg_p, retval_p )
 487   int_p d_p, ig_p, jg_p, retval_p ;
 488 {
 489   int d ;
 490   int kiddex ;
 491   int P ;
 492   int ig, jg ;
 493
 494   rsl_domain_info_t * info ;
 495   rsl_point_t       * domain ;
 496   ig = *ig_p - 1 ;
 497   jg = *jg_p - 1 ;
 498
 499   d = *d_p ;
 500   RSL_TEST_ERR( d < 0 || d > RSL_MAXDOMAINS,
 501     "rsl_ready_bcast: bad 'this hemi' descriptor" ) ;
 502   info = &( domain_info[d]) ;
 503 /* added 12/27/01 -- JM */
 504   if ( ig < 0 || ig >= info->len_m ||
 505        jg < 0 || jg >= info->len_n ) { *retval_p = 0 ; return ; }
 506   domain = info->domain ;
 507   kiddex = INDEX_2(jg,ig,info->len_m ) ;
 508   P = domain[ kiddex ].P ;
 509   *retval_p = 0 ;
 510   if ( P == rsl_myproc ) *retval_p = 1 ;
 511   return ;
 512 }
 513
 514 /* given a global point, return the processor number */
 515 RSL_PROC_FOR_POINT ( d_p, ig_p, jg_p, retval_p )
 516   int_p d_p, ig_p, jg_p, retval_p ;
 517 {
 518   int d ;
 519   int kiddex ;
 520   int P ;
 521   int ig, jg ;
 522
 523   rsl_domain_info_t * info ;
 524   rsl_point_t       * domain ;
 525   ig = *ig_p - 1 ;
 526   jg = *jg_p - 1 ;
 527
 528   d = *d_p ;
 529   RSL_TEST_ERR( d < 0 || d > RSL_MAXDOMAINS,
 530     "rsl_point_on_proc: bad descriptor" ) ;
 531   info = &( domain_info[d]) ;
 532   domain = info->domain ;
 533   kiddex = INDEX_2(jg,ig,info->len_m ) ;
 534   *retval_p = domain[ kiddex ].P ;
 535   return ;
 536 }
 537
 538