Distributed-Run-Time-Resource-Manager/RCCE_V2.0/apps/XHPL/src/pgesv/HPL_pdlaswp00T.c

/* 
 * -- High Performance Computing Linpack Benchmark (HPL)                
 *    HPL - 2.0 - September 10, 2008                          
 *    Antoine P. Petitet                                                
 *    University of Tennessee, Knoxville                                
 *    Innovative Computing Laboratory                                 
 *    (C) Copyright 2000-2008 All Rights Reserved                       
 *                                                                      
 * -- Copyright notice and Licensing terms:                             
 *                                                                      
 * Redistribution  and  use in  source and binary forms, with or without
 * modification, are  permitted provided  that the following  conditions
 * are met:                                                             
 *                                                                      
 * 1. Redistributions  of  source  code  must retain the above copyright
 * notice, this list of conditions and the following disclaimer.        
 *                                                                      
 * 2. Redistributions in binary form must reproduce  the above copyright
 * notice, this list of conditions,  and the following disclaimer in the
 * documentation and/or other materials provided with the distribution. 
 *                                                                      
 * 3. All  advertising  materials  mentioning  features  or  use of this
 * software must display the following acknowledgement:                 
 * This  product  includes  software  developed  at  the  University  of
 * Tennessee, Knoxville, Innovative Computing Laboratory.             
 *                                                                      
 * 4. The name of the  University,  the name of the  Laboratory,  or the
 * names  of  its  contributors  may  not  be used to endorse or promote
 * products  derived   from   this  software  without  specific  written
 * permission.                                                          
 *                                                                      
 * -- Disclaimer:                                                       
 *                                                                      
 * THIS  SOFTWARE  IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES,  INCLUDING,  BUT NOT
 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
 * A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE UNIVERSITY
 * OR  CONTRIBUTORS  BE  LIABLE FOR ANY  DIRECT,  INDIRECT,  INCIDENTAL,
 * SPECIAL,  EXEMPLARY,  OR  CONSEQUENTIAL DAMAGES  (INCLUDING,  BUT NOT
 * LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
 * DATA OR PROFITS; OR BUSINESS INTERRUPTION)  HOWEVER CAUSED AND ON ANY
 * THEORY OF LIABILITY, WHETHER IN CONTRACT,  STRICT LIABILITY,  OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
 * OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
 * ---------------------------------------------------------------------
 */ 
/*
 * Include files
 */
#include "hpl.h"

#ifdef STDC_HEADERS
void HPL_pdlaswp00T
(
   HPL_T_panel *                    PBCST,
   int *                            IFLAG,
   HPL_T_panel *                    PANEL,
   const int                        NN
)
#else
void HPL_pdlaswp00T
( PBCST, IFLAG, PANEL, NN )
   HPL_T_panel *                    PBCST;
   int *                            IFLAG;
   HPL_T_panel *                    PANEL;
   const int                        NN;
#endif
{
/* 
 * Purpose
 * =======
 *
 * HPL_pdlaswp00T applies the  NB  row interchanges to  NN columns of the
 * trailing submatrix and broadcast a column panel.
 *  
 * Bi-directional  exchange  is used to perform the  swap :: broadcast of
 * the row  panel U at once, resulting in a lower number of messages than
 * usual as well as a lower communication volume. With P process rows and
 * assuming  bi-directional links,  the running time of this function can
 * be approximated by:
 *  
 *    log_2(P) * (lat + NB*LocQ(N) / bdwth)
 *  
 * where  NB  is the number of rows of the row panel U,  N is the global
 * number of columns being updated,  lat and bdwth  are the latency  and
 * bandwidth  of  the  network  for  double  precision real words.  Mono
 * directional links will double this communication cost.
 *
 * Arguments
 * =========
 *
 * PBCST   (local input/output)          HPL_T_panel *
 *         On entry,  PBCST  points to the data structure containing the
 *         panel (to be broadcast) information.
 *
 * IFLAG   (local input/output)          int *
 *         On entry, IFLAG  indicates  whether or not  the broadcast has
 *         already been completed.  If not,  probing will occur, and the
 *         outcome will be contained in IFLAG on exit.
 *
 * PANEL   (local input/output)          HPL_T_panel *
 *         On entry,  PANEL  points to the data structure containing the
 *         panel (to be broadcast and swapped) information.
 *
 * NN      (local input)                 const int
 *         On entry, NN specifies  the  local  number  of columns of the
 *         trailing  submatrix  to  be swapped and broadcast starting at
 *         the current position. NN must be at least zero.
 *
 * ---------------------------------------------------------------------
 */ 
/*
 * .. Local Variables ..
 */
   MPI_Comm                  comm;
   HPL_T_grid                * grid;
   double                    * A, * U, * W;
   void                       * vptr = NULL;
   int                       * ipID, * lindxA, * lindxAU, * llen,
                             * llen_sv;
   unsigned int              ip2, ip2_=1, ipdist, ipow=1, mask=1,
                             mydist, mydis_;
   int                       Cmsgid=MSGID_BEGIN_PFACT, Np2, align,
                             hdim, i, icurrow, *iflag, ipA, ipW, *ipl,
                             iprow, jb, k, lda, ldW, myrow, n, nprow,
                             partner, root, size_, usize;
#define LDU                  n
/* ..
 * .. Executable Statements ..
 */
   n = Mmin( NN, PANEL->n ); jb = PANEL->jb;
/*
 * Quick return if there is nothing to do
 */
   if( ( n <= 0 ) || ( jb <= 0 ) ) return;

#ifdef HPL_DETAILED_TIMING
   HPL_ptimer( HPL_TIMING_LASWP );
#endif
/*
 * Retrieve parameters from the PANEL data structure
 */
   grid  = PANEL->grid;    nprow   = grid->nprow; myrow = grid->myrow;
   comm  = grid->col_comm; ip2     = (unsigned int)grid->row_ip2;
   hdim  = grid->row_hdim; align   = PANEL->algo->align;
   A     = PANEL->A;       U       = PANEL->U;    iflag = PANEL->IWORK;
   lda   = PANEL->lda;     icurrow = PANEL->prow; usize = jb * n;
   ldW   = n + 1;
/*
 * Allocate space for temporary W (ldW * jb)
 */
   vptr = (void*)malloc( ( (size_t)(align) + 
                           ((size_t)(jb) * (size_t)(ldW))) * 
                           sizeof(double) );
   if( vptr == NULL )
   { HPL_pabort( __LINE__, "HPL_pdlaswp00T", "Memory allocation failed" ); }

   W = (double *)HPL_PTR( vptr, ((size_t)(align) * sizeof(double) ) );
/*
 * Construct ipID and its local counter parts lindxA, lindxAU -  llen is
 * the number of rows/columns that I have in workspace and that I should
 * send.  Compute  lindx_, ipA, llen if it has not already been done for
 * this panel;
 */
   k = (int)((unsigned int)(jb) << 1); ipl = iflag + 1; ipID = ipl + 1;
   lindxA  = ipID + ((unsigned int)(k) << 1); lindxAU = lindxA + k;
   llen    = lindxAU + k; llen_sv = llen + nprow;

   if( *iflag == -1 )    /* no index arrays have been computed so far */
   {
      HPL_pipid(   PANEL,  ipl, ipID );
      HPL_plindx0( PANEL, *ipl, ipID, lindxA, lindxAU, llen_sv );
      *iflag = 0;
   }
   else if( *iflag == 1 ) /* HPL_pdlaswp01T called before: reuse ipID */
   {
      HPL_plindx0( PANEL, *ipl, ipID, lindxA, lindxAU, llen_sv );
      *iflag = 0;
   }
/*
 * Copy the llen_sv into llen - Reset ipA to its correct value
 */
   ipA = llen_sv[myrow];
   for( i = 0; i < nprow; i++ ) { llen[i]  = llen_sv[i]; }
/*
 * For i in [0..2*jb),  lindxA[i] is the offset in A of a row that ulti-
 * mately goes to U( lindxAU[i], : ) or U( :, lindxAU[i] ).  In icurrow,
 * we directly pack into U, otherwise we pack into workspace. The  first
 * entry of each column packed in workspace is in fact the row or column
 * offset in U where it should go to.
 */
   if( myrow == icurrow ) 
   {
      HPL_dlaswp01T( ipA, n, A, lda, U, LDU, lindxA, lindxAU );
   }
   else
   {
      HPL_dlaswp02N( ipA, n, A, lda, W, W+1, ldW, lindxA, lindxAU );
   }
/*
 * Probe for column panel - forward it when available 
 */
   if( *IFLAG == HPL_KEEP_TESTING ) (void) HPL_bcast( PBCST, IFLAG );
/*
 * Algorithm for bi-directional data exchange:
 *
 * As long as I have not talked to a process that  already  had the data
 * from icurrow,  I will be sending the workspace,  otherwise  I will be
 * sending U. Note that the columns in workspace contain the local index
 * in U they should go to.
 *
 * If I am receiving from a process that  has the data from  icurrow,  I
 * will be receiving in  U, copy the data of  U  that stays into  A, and
 * then the columns I have in workspace into U; otherwise  I will be re-
 * ceiving in the remaining workspace.  If I am one  of  those processes 
 * that already has the data from icurrow, I will be immediately copying
 * the data I have in my workspace into U.
 *
 * When I receive U, some of U should be copied in my piece of A  before
 * I can copy the rows I have in my workspace into  U.  This information
 * is kept in the lists  lindx_:  the row lindxAU[i] should be copied in
 * the row  lindxA[i] of my piece of  A, just as in the reversed initial
 * packing operation. Those rows are thus the first ones in the work ar-
 * ray.  After  this  operation  has  been  performed,  I will not  need
 * those lindx arrays,  and  I  will  always be sending a buffer of size
 * jb x n, or n x jb, that is, U.
 *
 * At  every  step  of  the algorithm, it is necesary to update the list 
 * llen,  so that I can figure out how large the next messages I will be
 * sending/receiving are.  It is  obvious when I am sending U. It is not
 * otherwise.
 *
 * We  choose  icurrow  to be the source of the bi-directional exchange.
 * This allows the processes in the non-power 2 part to receive U at the
 * first exchange,  and  then  broadcast internally this U so that those 
 * processes can grab their piece of A.
 */
   if( myrow == icurrow ) { llen[myrow] = 0; ipA = 0; }
   ipW    = ipA;
   Np2    = ( ( size_ = nprow - ip2 ) != 0 );
   mydist = (unsigned int)MModSub( myrow, icurrow, nprow );
/*
 * bi-directional exchange:   If nprow is not a power of 2,  proc[i-ip2]
 * receives local data from proc[i] for all i in  [ip2..nprow);  icurrow
 * is the source, these last process indexes are relative to icurrow.
 */
   if( ( Np2 != 0 ) && ( ( partner = (int)(mydist ^ ip2) ) < nprow ) )
   {
      partner = MModAdd( icurrow, partner, nprow );

      if( mydist == 0 )  /* I am the current row: I send U and recv W */
      {
         (void) HPL_sdrv( U, usize, Cmsgid, W, llen[partner] * ldW,
                          Cmsgid, partner, comm );
         if( llen[partner] > 0 )
            HPL_dlaswp03T( llen[partner], n, U, LDU, W, W+1, ldW );
      }
      else if( mydist == ip2 )
      {                      /* I recv U for later Bcast, I send my W */
         (void) HPL_sdrv( W, llen[myrow]*ldW, Cmsgid, U, usize,
                          Cmsgid, partner, comm );
      }
      else               /* None of us is icurrow, we exchange our Ws */
      {
         if( ( mydist & ip2 ) != 0 ) 
         {
            (void) HPL_send( W, llen[myrow]*ldW, partner, Cmsgid, comm );
         }
         else
         {
            (void) HPL_recv( Mptr( W, 0, ipW, ldW ), llen[partner]*ldW,
                             partner, Cmsgid, comm );
            if( llen[partner] > 0 ) ipW += llen[partner];
         }
      }
   }
/*
 * Update llen
 */
   for( i = 1; i < size_; i++ )
   {
      iprow   = MModAdd( icurrow, i,          nprow );
      partner = MModAdd( iprow,   (int)(ip2), nprow );
      llen[ iprow ] += llen[ partner ];
   }
/*
 * Probe for column panel - forward it when available 
 */
   if( *IFLAG == HPL_KEEP_TESTING ) (void) HPL_bcast( PBCST, IFLAG );
/*
 * power of 2 part of the processes collection:  only processes [0..ip2)
 * are working;  some of them  (mydist >> (k+1) == 0) either send or re-
 * ceive U.  At every step k, k is in [0 .. hdim),  of the algorithm,  a
 * process pair that exchanges  U  is such that  (mydist >> (k+1) == 0).
 * Among  those  processes,  the  ones  that are sending U are such that 
 * mydist >> k == 0.
 */
   if( mydist < ip2 )
   {
      k = 0;

      while( k < hdim )
      {
         partner = (int)(mydist ^ ipow);
         partner = MModAdd( icurrow, partner, nprow );
/*
 * Exchange and combine the local results - If I receive U,  then I must
 * copy from U the rows that belong to my piece of A, and then update  U
 * by  copying in it the rows I have accumulated in W.  Otherwise, I re-
 * ceive W.  In this later case, and I have U, I shall update my copy of
 * U by copying in it the rows I have accumulated in  W.  If  I  did not
 * have U before, I simply need to update my pointer in W for later use.
 */
         if( ( mydist >> (unsigned int)( k + 1 ) ) == 0 )
         {
            if( ( mydist >> (unsigned int)(k) ) == 0 )
            {
               (void) HPL_sdrv( U, usize, Cmsgid, Mptr( W, 0, ipW,
                                ldW ), llen[partner]*ldW, Cmsgid,
                                partner, comm );
               HPL_dlaswp03T( llen[partner], n, U, LDU, Mptr( W, 0, ipW,
                              ldW ), Mptr( W, 1, ipW, ldW ), ldW );
               ipW += llen[partner];
            }
            else
            {
               (void) HPL_sdrv( W, llen[myrow]*ldW, Cmsgid, U, usize,
                                Cmsgid, partner, comm );
               HPL_dlaswp04T( ipA, llen[myrow], n, U, LDU, A, lda, W,
                              W+1, ldW, lindxA, lindxAU );
            }
         }
         else
         {
            (void) HPL_sdrv( W, llen[myrow]*ldW, Cmsgid, Mptr( W, 0,
                             ipW, ldW ), llen[partner]*ldW, Cmsgid,
                             partner, comm );
            ipW += llen[partner];
         }
/*
 * Update llen - Go to next process pairs
 */
         iprow = icurrow; ipdist = 0;
         do
         {
            if( (unsigned int)( partner = (int)(ipdist ^ ipow) ) > ipdist )
            {
               partner = MModAdd( icurrow, partner, nprow );
               llen[iprow]  += llen[partner];
               llen[partner] = llen[iprow];
            }
            iprow = MModAdd( iprow, 1, nprow ); ipdist++;

         } while( ipdist < ip2 );

         ipow <<= 1; k++;
/*
 * Probe for column panel - forward it when available 
 */
         if( *IFLAG == HPL_KEEP_TESTING ) (void) HPL_bcast( PBCST, IFLAG );
      }
   }
   else
   {
/*
 * non power of 2 part of the process collection:  proc[ip2] broadcast U
 * to procs[ip2..nprow) (relatively to icurrow).
 */
      if( size_ > 1 )
      {
         k = size_ - 1;
         while( k > 1 ) { k >>= 1; ip2_ <<= 1; mask <<= 1; mask++; }
         root   = MModAdd( icurrow, (int)(ip2), nprow );
         mydis_ = (unsigned int)MModSub( myrow,  root, nprow );

         do
         {
            mask ^= ip2_;
            if( ( mydis_ & mask ) == 0 )
            {
               partner = (int)(mydis_ ^ ip2_);
               if( ( mydis_ & ip2_ ) != 0 )
               {
                  (void) HPL_recv( U, usize, MModAdd( root, partner,
                                   nprow ), Cmsgid, comm );

               }
               else if( partner < size_ )
               {
                  (void) HPL_send( U, usize, MModAdd( root, partner,
                                   nprow ), Cmsgid, comm );
               }
            }
            ip2_ >>= 1;
/*
 * Probe for column panel - forward it when available 
 */
            if( *IFLAG == HPL_KEEP_TESTING ) (void) HPL_bcast( PBCST, IFLAG );

         } while( ip2_ > 0 );
      }
/*
 * Every process in [ip2..nprow) (relatively to icurrow) grabs its piece
 * of A.
 */
      HPL_dlaswp05T( ipA, n, A, lda, U, LDU, lindxA, lindxAU );
   }
/*
 * If  nprow  is not a power of 2,  proc[i-ip2]  sends  global result to
 * proc[i] for all i in [ip2..nprow);
 */
   if( ( Np2 != 0 ) && ( ( partner = (int)(mydist ^ ip2) ) < nprow ) )
   {
      partner = MModAdd( icurrow, partner, nprow );
      if( ( mydist & ip2 ) != 0 )
      { (void) HPL_recv( U, usize, partner, Cmsgid, comm ); }
      else
      { (void) HPL_send( U, usize, partner, Cmsgid, comm ); }
   }

   if( vptr ) free( vptr );
/*
 * Probe for column panel - forward it when available 
 */
   if( *IFLAG == HPL_KEEP_TESTING ) (void) HPL_bcast( PBCST, IFLAG );

#ifdef HPL_DETAILED_TIMING
   HPL_ptimer( HPL_TIMING_LASWP );
#endif
/*
 * End of HPL_pdlaswp00T
 */
}