/* ---------------------------------------------------------------------
*
*  -- ScaLAPACK routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*  ---------------------------------------------------------------------
*/
/*
*  Include files
*/
#include "pblas.h"

void pbchkvect( n, npos0, ix, jx, desc_X, incx, dpos0, iix, jjx, ixrow,
                ixcol, nprow, npcol, myrow, mycol, info )
/*
*  .. Scalar Arguments ..
*/
   int         dpos0, * iix, incx, * info, ix, * ixcol, * ixrow, * jjx,
               jx, myrow, mycol, npcol, nprow, n, npos0;
/*
*  .. Array Arguments ..
*/
   int         desc_X[];
{
/*
*
*  Purpose
*  =======
*
*  pbchkvect checks the validity of a descriptor vector DESCX, the
*  related global indexes IX, JX and the global increment INCX. It also
*  computes the starting local indexes (IIX,JJX) corresponding to the
*  submatrix starting globally at the entry pointed by (IX,JX). 
*  Moreover, this routine returns the coordinates in the grid of the
*  process owning the global matrix entry of indexes (IX,JX), namely
*  (IXROW,IXCOL). The routine prevents out-of-bound memory access
*  by performing the appropriate MIN operation on iix and JJX.  Finally,
*  if an inconsistency is found among its parameters IX, JX, DESCX and
*  INCX, the routine returns an error code in info.
*
*  Arguments
*  =========
*
*  N       (global input) INTEGER
*          The length of the vector X being operated on.
*
*  NPOS0   (global input) INTEGER
*          Where in the calling routine's parameter list N appears.
*
*  IX      (global input) INTEGER
*          The row index in the global array X indicating the first
*          row of sub( X ).
*
*  JX      (global input) INTEGER
*          The column index in the global array X indicating the
*          first column of sub( X ).
*
*  DESCX   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix X.
*
*  INCX    (global input) INTEGER
*          The global increment for the elements of X. Only two values
*          of INCX are supported in this version, namely 1 and M_X. 
*          INCX must not be zero.
*
*  DPOS0   (global input) INTEGER
*          Where in the calling routine's parameter list DESCX
*          appears.  Note that we assume IX and JX are respectively 2
*          and 1 entries behind DESCX, and INCX is 1 entry after DESCX.
*
*  IIX     (local output) pointer to INTEGER
*          The local rows starting index of the submatrix.
*
*  JJX     (local output) pointer to INTEGER
*          The local columns starting index of the submatrix.
*
*  IXROW   (global output) pointer to INTEGER
*          The row coordinate of the process that possesses the first
*          row and column of the submatrix.
*
*  IXCOL   (global output) pointer to INTEGER
*          The column coordinate of the process that possesses the
*          first row and column of the submatrix.
*
*  NPROW   (global input) INTEGER
*          The total number of process rows over which the distributed
*          matrix is distributed.
*
*  NPCOL   (global input) INTEGER
*          The total number of process columns over which the
*          distributed matrix is distributed.
*
*  MYROW   (local input) INTEGER
*          The row coordinate of the process calling this routine.
*
*  MYCOL   (local input) INTEGER
*          The column coordinate of the process calling this routine.
*
*  INFO    (local input/local output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*
*  =====================================================================
*
*  .. Parameters ..
*/
#define DESCMULT      100
#define BIGNUM      10000
/* ..
*  .. Local Scalars ..
*/
   int         descpos, ExtraColBlock, ExtraRowBlock, icpos, ixpos,
               jxpos, MyColBlock, MyColDist, MyRowBlock, MyRowDist,
               NColBlock, np, npos, nq, NRowBlock;
/* ..
*  .. External Functions ..
*/
   F_INTG_FCT  numroc_();
/*
*  .. Executable Statements ..
*/
   if( *info >= 0 )
      *info = BIGNUM;
   else if( *info < -DESCMULT )
      *info = -(*info);
   else
      *info = -(*info) * DESCMULT;
/*
*  Figure where in parameter list each parameter was, factoring in
*  descriptor multiplier
*/
   npos = npos0 * DESCMULT;
   ixpos = ( dpos0 - 2 ) * DESCMULT; 
   jxpos = ( dpos0 - 1 ) * DESCMULT; 
   icpos = ( dpos0 + 1 ) * DESCMULT; 
   descpos = dpos0 * DESCMULT + 1;
/*
 * Check that we have a legal descriptor type
 */
   if(desc_X[DT_] != BLOCK_CYCLIC_2D) *info = MIN( *info, descpos + DT_ );
/*
*  Check that matrix values make sense from local viewpoint
*/
   if( n < 0 )
      *info = MIN( *info, npos );
   else if( ix < 1 )
      *info = MIN( *info, ixpos );
   else if( jx < 1 )
      *info = MIN( *info, jxpos );
   else if( desc_X[MB_] < 1 )
      *info = MIN( *info, descpos + MB_ );
   else if( desc_X[NB_] < 1 )
      *info = MIN( *info, descpos + NB_ );
   else if( ( desc_X[RSRC_] < 0 ) || ( desc_X[RSRC_] >= nprow ) )
      *info = MIN( *info, descpos + RSRC_ );
   else if( ( desc_X[CSRC_] < 0 ) || ( desc_X[CSRC_] >= npcol ) )
      *info = MIN( *info, descpos + CSRC_ );
   else if( incx != 1 && incx != desc_X[M_] )
      *info = MIN( *info, icpos );
   else if( desc_X[LLD_] < 1 )
      *info = MIN( *info, descpos + LLD_ );

   if( n == 0 )
   {
/*
*     NULL matrix, relax some checks
*/
      if( desc_X[M_] < 0 )
         *info = MIN( *info, descpos + M_ );
      if( desc_X[N_] < 0 )
         *info = MIN( *info, descpos + N_ );
   }
   else
   {
/*
*     more rigorous checks for non-degenerate matrices
*/
      if( desc_X[M_] < 1 )
         *info = MIN( *info, descpos + M_ );
      else if( desc_X[N_] < 1 )
         *info = MIN( *info, descpos + N_ );
      else if( ( incx == desc_X[M_] ) && ( jx+n-1 > desc_X[N_] ) )
         *info = MIN( *info, jxpos );
      else if( ( incx == 1 ) && ( incx != desc_X[M_] ) &&
               ( ix+n-1 > desc_X[M_] ) )
         *info = MIN( *info, ixpos );
      else
      {
         if( ix > desc_X[M_] )
            *info = MIN( *info, ixpos );
         else if( jx > desc_X[N_] )
            *info = MIN( *info, jxpos );
      }
   }
/*
*  Retrieve local information for vector X, and prepare output:
*  set info = 0 if no error, and divide by DESCMULT if error is not
*  in a descriptor entry.
*/
   if( *info == BIGNUM )
   {
      MyRowDist = ( myrow + nprow - desc_X[RSRC_] ) % nprow;
      MyColDist = ( mycol + npcol - desc_X[CSRC_] ) % npcol;
      NRowBlock = desc_X[M_] / desc_X[MB_];
      NColBlock = desc_X[N_] / desc_X[NB_];
      np = ( NRowBlock / nprow ) * desc_X[MB_];
      nq = ( NColBlock / npcol ) * desc_X[NB_];
      ExtraRowBlock = NRowBlock % nprow;
      ExtraColBlock = NColBlock % npcol;

      ix--;
      jx--;
      MyRowBlock = ix / desc_X[MB_];
      MyColBlock = jx / desc_X[NB_];
      *ixrow = ( MyRowBlock + desc_X[RSRC_] ) % nprow;
      *ixcol = ( MyColBlock + desc_X[CSRC_] ) % npcol;

      *iix = ( MyRowBlock / nprow + 1 ) * desc_X[MB_] + 1;
      *jjx = ( MyColBlock / npcol + 1 ) * desc_X[NB_] + 1;

      if( MyRowDist >= ( MyRowBlock % nprow ) )
      {
         if( myrow == *ixrow )
            *iix += ix % desc_X[MB_];
         *iix -= desc_X[MB_];
      }
      if( MyRowDist  < ExtraRowBlock )
         np += desc_X[MB_];
      else if( MyRowDist == ExtraRowBlock )
         np += ( desc_X[M_] % desc_X[MB_] );
      np = MAX( 1, np );

      if( MyColDist >= ( MyColBlock % npcol ) )
      {
         if( mycol == *ixcol )
            *jjx += jx % desc_X[NB_];
         *jjx -= desc_X[NB_];
      }
      if( MyColDist < ExtraColBlock )
         nq += desc_X[NB_];
      else if( MyColDist == ExtraColBlock )
         nq += ( desc_X[N_] % desc_X[NB_] );
      nq = MAX( 1, nq );

      *iix = MIN( *iix, np );
      *jjx = MIN( *jjx, nq );

      if( desc_X[LLD_] < np )
      {
         if( numroc_(&desc_X[N_], &desc_X[NB_], &mycol, &desc_X[CSRC_], &npcol) )
            *info = -( descpos + LLD_ );
         else *info = 0;
      }
      else *info = 0;
   }
   else if( *info % DESCMULT == 0 )
   {
      *info = -(*info) / DESCMULT;
   }
   else
   {
      *info = -(*info);
   }
}

/* ---------------------------------------------------------------------
*
*  -- ScaLAPACK routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*  ---------------------------------------------------------------------
*/

void pbchkmat( m, mpos0, n, npos0, ia, ja, desc_A, dpos0, iia, jja,
               iarow, iacol, nprow, npcol, myrow, mycol, nra, nca,
               info )
/*
*  .. Scalar Arguments ..
*/
   int         dpos0, * iia, * info, ia, * iacol, * iarow, * jja,
               ja, m, mpos0, myrow, mycol, npcol, nprow, n, * nca,
               npos0, * nra;
/*
*  .. Array Arguments ..
*/
   int         desc_A[];
{
/*
*
*  Purpose
*  =======
*
*  pbmatvect checks the validity of a descriptor vector DESCA, the
*  related global indexes IA, JA. It also computes the starting local 
*  indexes (IIA,JJA) corresponding to the submatrix starting globally at
*  the entry pointed by (IA,JA). Moreover, this routine returns the 
*  coordinates in the grid of the process owning the global matrix entry
*  of indexes (IA,JA), namely (IAROW,IACOL). The routine prevents from
*  out-of-bound memory access, by performing the adequate MIN operation
*  on IIA and JJA.  Finally, if an inconsitency is found among its 
*  parameters ia, ja and desc_A, the routine returns an error code in
*  info.
*
*  Arguments
*  =========
*
*  M       (global input) INTEGER
*          The number or matrix rows of A being operated on.
*
*  MPOS0   (global input) INTEGER
*          Where in the calling routine's parameter list M appears.
*
*  N       (global input) INTEGER
*          The number or matrix columns of A being operated on.
*
*  NPOS0   (global input) INTEGER
*          Where in the calling routine's parameter list N appears.
*
*  IA      (global input) INTEGER
*          The row index in the global array A indicating the first
*          row of sub( A ).
*
*  JA      (global input) INTEGER
*          The column index in the global array A indicating the
*          first column of sub( A ).
*
*  DESCA   (global and local input) INTEGER array of dimension DLEN_.
*          The array descriptor for the distributed matrix A.
*
*  DPOS0   (global input) INTEGER
*          Where in the calling routine's parameter list DESCA
*          appears.  Note that we assume IA and JA are respectively 2
*          and 1 entries behind DESCA.
*
*  IIA     (local output) pointer to INTEGER
*          The local rows starting index of the submatrix.
*
*  JJA     (local output) pointer to INTEGER
*          The local columns starting index of the submatrix.
*
*  IAROW   (global output) pointer to INTEGER
*          The row coordinate of the process that possesses the first
*          row and column of the submatrix.
*
*  IACOL   (global output) pointer to INTEGER
*          The column coordinate of the process that possesses the
*          first row and column of the submatrix.
*
*  NPROW   (global input) INTEGER
*          The total number of process rows over which the distributed
*          matrix is distributed.
*
*  NPCOL   (global input) INTEGER
*          The total number of process columns over which the
*          distributed matrix is distributed.
*
*  MYROW   (local input) INTEGER
*          The row coordinate of the process calling this routine.
*
*  MYCOL   (local input) INTEGER
*          The column coordinate of the process calling this routine.
*
*  NRA     (local output) pointer to INTEGER
*          The maximum between the number of local rows owned by the
*          calling process and 1.
*
*  NCA     (local output) pointer to INTEGER
*          The maximum between the number of local columns owned by the
*          calling process and 1.
*
*  INFO    (local input/local output) INTEGER
*          = 0:  successful exit
*          < 0:  If the i-th argument is an array and the j-entry had
*                an illegal value, then INFO = -(i*100+j), if the i-th
*                argument is a scalar and had an illegal value, then
*                INFO = -i.
*
*  =====================================================================
*
*  .. Parameters ..
*/
#define DESCMULT      100
#define BIGNUM      10000
/* ..
*  .. Local Scalars ..
*/
   int         descpos, ExtraColBlock, ExtraRowBlock, iapos, japos,
               mpos, MyColBlock, MyColDist, MyRowBlock, MyRowDist,
               NColBlock, npos, NRowBlock;
/* ..
*  .. External Functions ..
*/
   F_INTG_FCT  numroc_();
/*
*  .. Executable Statements ..
*/
   if( *info >= 0 )
      *info = BIGNUM;
   else if( *info < -DESCMULT )
      *info = -(*info);
   else
      *info = -(*info) * DESCMULT;
/*
*  Figure where in parameter list each parameter was, factoring in
*  descriptor multiplier
*/
   mpos = mpos0 * DESCMULT;
   npos = npos0 * DESCMULT;
   iapos = ( dpos0 - 2 ) * DESCMULT; 
   japos = ( dpos0 - 1 ) * DESCMULT; 
   descpos = dpos0 * DESCMULT + 1;
/*
 * Check that we have a legal descriptor type
 */
   if(desc_A[DT_] != BLOCK_CYCLIC_2D) *info = MIN( *info, descpos + DT_ );
/*
*  Check that matrix values make sense from local viewpoint
*/
   if( m < 0 )
      *info = MIN( *info, mpos );
   if( n < 0 )
      *info = MIN( *info, npos );
   else if( ia < 1 )
      *info = MIN( *info, iapos );
   else if( ja < 1 )
      *info = MIN( *info, japos );
   else if( desc_A[MB_] < 1 )
      *info = MIN( *info, descpos + MB_ );
   else if( desc_A[NB_] < 1 )
      *info = MIN( *info, descpos + NB_ );
   else if( ( desc_A[RSRC_] < 0 ) || ( desc_A[RSRC_] >= nprow ) )
      *info = MIN( *info, descpos + RSRC_ );
   else if( ( desc_A[CSRC_] < 0 ) || ( desc_A[CSRC_] >= npcol ) )
      *info = MIN( *info, descpos + CSRC_ );
   else if( desc_A[LLD_] < 1 )
      *info = MIN( *info, descpos + LLD_ );

   if( m == 0 || n == 0 )
   {
/*
*     NULL matrix, relax some checks
*/
      if( desc_A[M_] < 0 )
         *info = MIN( *info, descpos + M_ );
      if( desc_A[N_] < 0 )
         *info = MIN( *info, descpos + N_ );
   }
   else
   {
/*
*     more rigorous checks for non-degenerate matrices
*/
      if( desc_A[M_] < 1 )
         *info = MIN( *info, descpos + M_ );
      else if( desc_A[N_] < 1 )
         *info = MIN( *info, descpos + N_ );
      else
      {
         if( ia > desc_A[M_] )
            *info = MIN( *info, iapos );
         else if( ja > desc_A[N_] )
            *info = MIN( *info, japos );
         else
         {
            if( ia+m-1 > desc_A[M_] )
               *info = MIN( *info, mpos );
            if( ja+n-1 > desc_A[N_] )
               *info = MIN( *info, npos );
         }
      }
   }
/*
*  Retrieve local information for matrix A, and prepare output:
*  set info = 0 if no error, and divide by DESCMULT if error is not
*  in a descriptor entry.
*/
   if( *info == BIGNUM )
   {
      MyRowDist = ( myrow + nprow - desc_A[RSRC_] ) % nprow;
      MyColDist = ( mycol + npcol - desc_A[CSRC_] ) % npcol;
      NRowBlock = desc_A[M_] / desc_A[MB_];
      NColBlock = desc_A[N_] / desc_A[NB_];
      *nra = ( NRowBlock / nprow ) * desc_A[MB_];
      *nca = ( NColBlock / npcol ) * desc_A[NB_];
      ExtraRowBlock = NRowBlock % nprow;
      ExtraColBlock = NColBlock % npcol;

      ia--;
      ja--;
      MyRowBlock = ia / desc_A[MB_];
      MyColBlock = ja / desc_A[NB_];
      *iarow = ( MyRowBlock + desc_A[RSRC_] ) % nprow;
      *iacol = ( MyColBlock + desc_A[CSRC_] ) % npcol;

      *iia = ( MyRowBlock / nprow + 1 ) * desc_A[MB_] + 1;
      *jja = ( MyColBlock / npcol + 1 ) * desc_A[NB_] + 1;

      if( MyRowDist >= ( MyRowBlock % nprow ) )
      {
         if( myrow == *iarow )
            *iia += ia % desc_A[MB_];
         *iia -= desc_A[MB_];
      }
      if( MyRowDist  < ExtraRowBlock )
         *nra += desc_A[MB_];
      else if( MyRowDist == ExtraRowBlock )
         *nra += ( desc_A[M_] % desc_A[MB_] );
      *nra = MAX( 1, *nra );

      if( MyColDist >= ( MyColBlock % npcol ) )
      {
         if( mycol == *iacol )
            *jja += ja % desc_A[NB_];
         *jja -= desc_A[NB_];
      }
      if( MyColDist < ExtraColBlock )
         *nca += desc_A[NB_];
      else if( MyColDist == ExtraColBlock )
         *nca += ( desc_A[N_] % desc_A[NB_] );
      *nca = MAX( 1, *nca );

      *iia = MIN( *iia, *nra );
      *jja = MIN( *jja, *nca );

      if( desc_A[LLD_] < *nra )
      {
         if( numroc_(&desc_A[N_], &desc_A[NB_], &mycol, &desc_A[CSRC_], &npcol) )
            *info = -( descpos + LLD_ );
         else *info = 0;
      }
      else
         *info = 0;
   }
   else if( *info % DESCMULT == 0 )
   {
      *info = -(*info) / DESCMULT;
   }
   else
   {
      *info = -(*info);
   }
}

/* ---------------------------------------------------------------------
*
*  -- ScaLAPACK routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*  ---------------------------------------------------------------------
*/

char * ptop( op, scope, top )
/*
*  .. Scalar Arguments ..
*/
   char        * op, * scope, * top;
{
/* 
*  Purpose
*  =======
*
*  ptop get or set the row, column or all broadcast or combine 
*  topologies.
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   static char rowbtop = ' ';         /* Default broadcast topologies */
   static char colbtop = ' ';
   static char allbtop = ' ';

   static char rowctop = ' ';           /* Default combine topologies */
   static char colctop = ' ';
   static char allctop = ' ';

   if( *op == 'B' )
   {
      if( *top == '!' )
      {
         if( *scope == 'R' )
         {
            return &rowbtop;
         }
         else if( *scope == 'C' ) 
         {
            return &colbtop;
         }
         else
         {
            return &allbtop;
         }
      }
      else
      {
         if( *scope == 'R' )
         {
            rowbtop = *top;
            return &rowbtop;
         }
         else if( *scope == 'C' )
         {
            colbtop = *top;
            return &colbtop;
         }
         else
         {
            allbtop = *top;
            return &allbtop;
         }
      }
   }
   else
   { 
      if( *top == '!' )
      {
         if( *scope == 'R' )
         {
            return &rowctop;
         }
         else if( *scope == 'C' )
         {
            return &colctop;
         }
         else
         {
            return &allctop;
         }
      }
      else
      {
         if( *scope == 'R' )
         {
            rowctop = *top;
            return &rowctop;
         }
         else if( *scope == 'C' )
         {
            colctop = *top;
            return &colctop;
         }
         else
         {
            allctop = *top;
            return &allctop;
         }
      }
   }
}

/* ---------------------------------------------------------------------
*
*  -- ScaLAPACK routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*  ---------------------------------------------------------------------
*/

void ptopset_( ictxt, op, scope, top )
/*
*  .. Scalar Arguments ..
*/

   int         *ictxt;
   F_CHAR      op, scope, top;
{
/*
*  Purpose
*  =======
*
*  ptopset_ set the row, column or all broadcast and combine topologies.
*
* ======================================================================
*
*  .. External Functions ..
*/
   char        * ptop();
/*
*  .. Executable Statements ..
*/
   if( * F2C_CHAR( top ) != '!' )
      ptop( F2C_CHAR( op ), F2C_CHAR( scope ), F2C_CHAR( top ) );
}

/* ---------------------------------------------------------------------
*
*  -- ScaLAPACK routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*  ---------------------------------------------------------------------
*/

void ptopget_( ictxt, op, scope, top  )
/*
*  .. Scalar Arguments ..
*/
   int         *ictxt;
   F_CHAR      op, scope, top;
{
/*
*  Purpose
*  =======
*
*  ptopget_ get the row, column or all broadcast and combine topologies.
*
* ======================================================================
*
*  .. External Functions ..
*/
   char        * ptop(), * topc;
/*
*  .. Executable Statements ..
*/
   topc = F2C_CHAR( top );
   * topc = * ptop( F2C_CHAR( op ), F2C_CHAR( scope ), "!" );
   top = C2F_CHAR( topc );
}

char * getpbbuf( mess, length )
/*
*  .. Scalar Arguments ..
*/
   int         length;
/*
*  .. Array Arguments ..
*/
   char        * mess;
{
/*
*  Purpose
*  =======
*
*  getpbbuf returns a pointer to a working buffer of size length alloca-
*  ted for the PBLAS routines.
*
* ======================================================================
*
*  .. Local Scalars ..
*/
   static char * pblasbuf = NULL;
   static int  pbbuflen = 0, mone = -1;
/* ..
*  .. External Functions ..
*/
   void        blacs_abort_();
/* ..
*  .. Executable Statements ..
*/
   if( length >= 0 )
   {
      if( length > pbbuflen )
      {
         if( pblasbuf ) 
            free( pblasbuf );
         pblasbuf = (char *) malloc((unsigned)length);
         if( !pblasbuf )
         {
            fprintf( stderr, 
                     "PBLAS %s ERROR: Memory allocation failed\n",
                     mess );
            blacs_abort_( &mone, &mone );
         }
         pbbuflen = length;
      }
   }
   else if( pblasbuf )
   {
      free( pblasbuf );
      pblasbuf = NULL;
      pbbuflen = 0;
   }
   return( pblasbuf );
}

/* ---------------------------------------------------------------------
*
*  -- ScaLAPACK routine (version 1.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     November 17, 1996
*
*  ---------------------------------------------------------------------
*/

void pbfreebuf_()
{
/*
*  Purpose
*  =======
*
*  pbfreebuf_ disposes the buffer allocated for the PBLAS routines.
*
*  ======================================================================
*
*  .. External Functions ..
*/
   char        * getpbbuf();
/* ..
*  .. Executable Statements ..
*/
   getpbbuf( " ", -1 );                          /* free PBLAS buffer */
}