/* ---------------------------------------------------------------------
*
*  -- 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 pstran_( m, n, alpha, A, ia, ja, desc_A, beta, C, ic, jc, desc_C )
/*
*  .. Scalar Arguments ..
*/
   int         * ia, * ic, * ja, * jc, * m, * n;
   float       * alpha, * beta;
/* ..
*  .. Array Arguments ..
*/
   int         desc_A[], desc_C[];
   float       A[], C[];
{
/*
*  Purpose
*  =======
*
*  PSTRAN  transposes a distributed matrix
*
*     sub( C ) = beta * sub( C ) + alpha * op( sub( A ) )
*
*  where sub( C ) denotes C(IC:IC+M-1,JC:JC+N-1),
*        sub( A ) denotes A(IA:IA+N-1,JA:JA+M-1),
*        op( A )  denotes A**T.
*
*  Beta is a scalar, sub( C ) is an M-by-N distributed matrix, sub( A )
*  is an N-by-M distributed matrix.
*
*  Notes
*  =====
*
*  Each global data object is described by an associated description
*  vector.  This vector stores the information required to establish
*  the mapping between an object element and its corresponding process
*  and memory location.
*
*  Let A be a generic term for any 2D block cyclicly distributed array.
*  Such a global array has an associated description vector descA.
*  In the following comments, the character _ should be read as
*  "of the global array".
*
*  NOTATION        STORED IN      EXPLANATION
*  --------------- -------------- --------------------------------------
*  DT_A   (global) descA[ DT_ ]   The descriptor type.  In this case,
*                                 DT_A = 1.
*  CTXT_A (global) descA[ CTXT_ ] The BLACS context handle, indicating
*                                 the BLACS process grid A is distribu-
*                                 ted over. The context itself is glo-
*                                 bal, but the handle (the integer
*                                 value) may vary.
*  M_A    (global) descA[ M_ ]    The number of rows in the global
*                                 array A.
*  N_A    (global) descA[ N_ ]    The number of columns in the global
*                                 array A.
*  MB_A   (global) descA[ MB_ ]   The blocking factor used to distribu-
*                                 te the rows of the array.
*  NB_A   (global) descA[ NB_ ]   The blocking factor used to distribu-
*                                 te the columns of the array.
*  RSRC_A (global) descA[ RSRC_ ] The process row over which the first
*                                 row of the array A is distributed.
*  CSRC_A (global) descA[ CSRC_ ] The process column over which the
*                                 first column of the array A is
*                                 distributed.
*  LLD_A  (local)  descA[ LLD_ ]  The leading dimension of the local
*                                 array.  LLD_A >= MAX(1,LOCr(M_A)).
*
*  Let K be the number of rows or columns of a distributed matrix,
*  and assume that its process grid has dimension p x q.
*  LOCr( K ) denotes the number of elements of K that a process
*  would receive if K were distributed over the p processes of its
*  process column.
*  Similarly, LOCc( K ) denotes the number of elements of K that a
*  process would receive if K were distributed over the q processes of
*  its process row.
*  The values of LOCr() and LOCc() may be determined via a call to the
*  ScaLAPACK tool function, NUMROC:
*          LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
*          LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
*  An upper bound for these quantities may be computed by:
*          LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
*          LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
*
*  If M+MOD(JA-1,NB_A) <= NB_A, sub( A ) is considered to be distributed
*  columnwise, IA-1 must be a multiple of MB_A, JC-1 must be a multiple
*  of NB_C and the active block sizes must be equal: MB_A = NB_C,
*  else if N+MOD(IA-1,MB_A) <= MB_A, sub( A ) is considered to be
*  distributed rowwise, JA-1 must be a multiple of NB_A, IC-1 must be a
*  multiple of MB_C and the active block sizes must be equal:
*  NB_A = MB_C.
*  else if M < N, sub( A ) is considered to be distributed columnwise
*  for looping purposes, IA-1 must be a multiple of MB_A, JC-1 must be a
*  multiple of NB_C, MOD(JA-1,NB_A) must be equal to MOD(IC-1,MB_C) and
*  the active block sizes must be equal: MB_A = NB_C and NB_A = MB_C,
*  otherwise sub( A ) is considered to be distributed rowwise for
*  looping purposes, JA-1 must be a multiple of NB_A, IC-1 must be a
*  multiple of MB_C, MOD(IA-1,MB_A) must be equal to MOD(JC-1,NB_C) and
*  the active block sizes must be equal: NB_A = MB_C and MB_A = NB_C.
*
*  Arguments
*  =========
*
*  M       (global input) pointer to INTEGER
*          The number of rows to be operated on i.e the number of rows
*          of the distributed submatrix sub( C ). M >= 0.
*
*  N       (global input) pointer to INTEGER
*          The number of columns to be operated on i.e the number of
*          columns of the distributed submatrix sub( C ). N >= 0.
*
*  ALPHA   (global input) pointer to REAL
*          On entry,  ALPHA  specifies the scalar  alpha.
*
*  A       (local input) REAL pointer into the local memory
*          to an array of dimension (LLD_A,LOCc(JA+M-1)) containing the
*          local pieces of the distributed matrix sub( A ).
*
*  IA      (global input) pointer to INTEGER
*          The global row index of the submatrix of the distributed
*          matrix A to operate on.
*
*  JA      (global input) pointer to INTEGER
*          The global column index of the submatrix of the distributed
*          matrix A to operate on.
*
*  DESCA   (global and local input) INTEGER array of dimension 8.
*          The array descriptor of the distributed matrix A.
*
*  BETA    (global input) pointer to REAL
*          On entry, BETA  specifies the scalar beta. When  BETA
*          is supplied as zero then sub( C ) need not be set on input.
*
*  C       (local input/local output) REAL pointer into the
*          local memory to an array of dimension (LLD_C,LOCc(JC+N-1))
*          containing the local pieces of the distributed matrix
*          sub( C ). On exit, the distributed matrix sub( C ) is over-
*          written by the updated matrix.
*
*  IC      (global input) pointer to INTEGER
*          The global row index of the submatrix of the distributed
*          matrix C to operate on.
*
*  JC      (global input) pointer to INTEGER
*          The global column index of the submatrix of the distributed
*          matrix C to operate on.
*
*  DESCC   (global and local input) INTEGER array of dimension 8.
*          The array descriptor of the distributed matrix C.
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   char        adist;
   int         cablock, ccblock, i, iacol, iarow, iblk, iccol, icoffa,
               icoffc, icrow, ictxt, iia, iic, in, info, iroffa, iroffc,
               j, jblk, jja, jjc, jn, lcm, lcmp, lcmq, mp0, mycol,
               myrow, nca, ncc, nq0, nprow, npcol, nra, nrc, rablock,
               rcblock, tmp0, wksz;
/* ..
*  .. PBLAS Buffer ..
*/
   float       * buff;
/* ..
*  .. External Functions ..
*/
   void        blacs_gridinfo_();
   void        pbchkmat();
   void        pberror_();
   char        * getpbbuf();
   void        pstra2();
   F_INTG_FCT  ilcm_();
/* ..
*  .. Executable Statements ..
*
*  Get grid parameters
*/
   ictxt = desc_A[CTXT_];
   blacs_gridinfo_( &ictxt, &nprow, &npcol, &myrow, &mycol );
/*
*  Test the input parameters
*/
   info = 0;
   if( nprow == -1 )
      info = -(700+CTXT_+1);
   else
   {
      pbchkmat( *n, 2, *m, 1, *ia, *ja, desc_A, 7, &iia, &jja,
                &iarow, &iacol, nprow, npcol, myrow, mycol,
                &nra, &nca, &info );
      pbchkmat( *m, 1, *n, 2, *ic, *jc, desc_C, 12, &iic, &jjc,
                &icrow, &iccol, nprow, npcol, myrow, mycol,
                &nrc, &ncc, &info );
      if( info == 0 )
      {
         iroffa = (*ia-1) % desc_A[MB_];
         icoffa = (*ja-1) % desc_A[NB_];
         iroffc = (*ic-1) % desc_C[MB_];
         icoffc = (*jc-1) % desc_C[NB_];
         cablock = ( ( (*m) + icoffa ) <= desc_A[NB_] );
         rablock = ( ( (*n) + iroffa ) <= desc_A[MB_] );
         if( cablock )
            adist = 'C';
         else if( rablock )
            adist = 'R';
         else if( *m <= *n )
            adist = 'C';
         else
            adist = 'R';
         if( adist == 'C' )
         {
            rcblock = ( ( (*m) + iroffc ) <= desc_C[MB_] );
            if( iroffa != 0 )
               info = -5;
            else if( icoffc != 0 )
               info = -11;
            else if( desc_A[MB_] != desc_C[NB_] )
               info = -(1200+NB_+1);
            if( !cablock || !rcblock )
            {
               if( icoffa != iroffc )
                  info = -10;
               else if( desc_A[NB_] != desc_C[MB_] )
                  info = -(1200+MB_+1);
            }
         }
         else
         {
            ccblock = ( ( (*n) + icoffc ) <= desc_C[NB_] );
            if( icoffa != 0 )
               info = -6;
            else if( iroffc != 0 )
               info = -10;
            else if( desc_A[NB_] != desc_C[MB_] )
               info = -(1200+MB_+1);
            if( !rablock || !ccblock )
            {
               if( iroffa != icoffc )
                  info = -11;
               else if( desc_A[MB_] != desc_C[NB_] )
               info = -(1200+NB_+1);
            }
         }
         if( ictxt != desc_C[CTXT_] )
            info = -(1200+CTXT_+1);
      }
   }
   if( info )
   {
      pberror_( &ictxt, "PSTRAN", &info );
      return;
   }
/*
*  Quick return if possible.
*/
   if( ( *m == 0 ) || ( *n == 0 ) ||
       ( ( *alpha == ZERO ) && ( *beta == ONE ) ) )
      return;
/*
*  Figure out the arguments to be passed to pstra2
*/
   lcm = ilcm_( &nprow, &npcol );
   if( adist == 'C' )
   {
      lcmq = lcm / npcol;
      tmp0 = (*n) / desc_A[MB_];
      nq0 = MYROC0( tmp0, *n, desc_A[MB_], npcol );
      tmp0 = nq0 / desc_A[MB_];
      wksz = desc_A[NB_] * MYROC0( tmp0, nq0, desc_A[MB_], lcmq );
   }
   else
   {
      lcmp = lcm / nprow;
      tmp0 = (*m) / desc_A[NB_];
      mp0 = MYROC0( tmp0, *m, desc_A[NB_], nprow );
      tmp0 = mp0 / desc_A[NB_];
      wksz = desc_A[MB_] * MYROC0( tmp0, mp0, desc_A[NB_], lcmp );
   }
   buff = (float *)getpbbuf( "PSTRAN", wksz*sizeof(float) );
/*
*  Call PB-BLAS routine
*/
   if( adist == 'R' )
   {
      i = CEIL( (*ia), desc_A[MB_] ) * desc_A[MB_];
      in = (*ia)+(*n)-1;
      in = MIN( i, in );
                                     /* Handle first block separately */
      iblk = in-(*ia)+1;
      pstra2( &ictxt, &adist, &iblk, m, &desc_A[NB_], alpha,
              &A[iia-1+(jja-1)*desc_A[LLD_]], &desc_A[LLD_],
              beta, &C[iic-1+(jjc-1)*desc_C[LLD_]], &desc_C[LLD_],
              &iarow, &iacol, &icrow, &iccol, buff );
      if( myrow == iarow )
      {
         iia += iblk;
         iia = MIN( iia, nra );
      }
      if( mycol == iccol )
      {
         jjc += iblk;
         jjc = MIN( jjc, ncc );
      }
      iarow = (iarow+1) % nprow;
      iccol = (iccol+1) % npcol;
                                 /* loop over remaining block of rows */
      tmp0 = (*ia)+(*n)-1;
      for( i=in+1; i <= tmp0; i+=desc_A[MB_] )
      {
          iblk = (*n)-i+(*ia);
          iblk = MIN( desc_A[MB_], iblk );
          pstra2( &ictxt, &adist, &iblk, m, &desc_A[NB_], alpha,
                  &A[iia-1+(jja-1)*desc_A[LLD_]], &desc_A[LLD_],
                  beta, &C[iic-1+(jjc-1)*desc_C[LLD_]], &desc_C[LLD_],
                  &iarow, &iacol, &icrow, &iccol, buff );

          if( myrow == iarow )
          {
             iia += iblk;
             iia = MIN( iia, nra );
          }
          if( mycol == iccol )
          {
             jjc += iblk;
             jjc = MIN( jjc, ncc );
          }
          iarow = (iarow+1) % nprow;
          iccol = (iccol+1) % npcol;
      }
   }
   else
   {
      j = CEIL( (*ja), desc_A[NB_] ) * desc_A[NB_];
      jn = (*ja)+(*m)-1;
      jn = MIN( j, jn );
                                     /* Handle first block separately */
      jblk = jn-(*ja)+1;
      pstra2( &ictxt, &adist, n, &jblk, &desc_A[MB_], alpha,
              &A[iia-1+(jja-1)*desc_A[LLD_]], &desc_A[LLD_], beta,
              &C[iic-1+(jjc-1)*desc_C[LLD_]], &desc_C[LLD_], &iarow,
              &iacol, &icrow, &iccol, buff );

      if( myrow == icrow )
      {
         iic += jblk;
         iic = MIN( iic, nrc );
      }
      if( mycol == iacol )
      {
         jja += jblk;
         jja = MIN( jja, nca );
      }
      icrow = (icrow+1) % nprow;
      iacol = (iacol+1) % npcol;
                              /* loop over remaining block of columns */
      tmp0 = (*ja)+(*m)-1;
      for( j=jn+1; j <= tmp0; j+=desc_A[NB_] )
      {
          jblk = (*m)-j+(*ja);
          jblk = MIN( desc_A[NB_], jblk );
          pstra2( &ictxt, &adist, n, &jblk, &desc_A[MB_], alpha,
                  &A[iia-1+(jja-1)*desc_A[LLD_]], &desc_A[LLD_], beta,
                  &C[iic-1+(jjc-1)*desc_C[LLD_]], &desc_C[LLD_], &iarow,
                  &iacol, &icrow, &iccol, buff );
          if( myrow == icrow )
          {
             iic += jblk;
             iic = MIN( iic, nrc );
          }
          if( mycol == iacol )
          {
             jja += jblk;
             jja = MIN( jja, nca );
          }
          icrow = (icrow+1) % nprow;
          iacol = (iacol+1) % npcol;
      }
   }
}

void pstra2( ictxt, adist, m, n, nb, alpha, A, lda, beta, C, ldc,
             iarow, iacol, icrow, iccol, WORK )
/*
*  .. Scalar Arguments ..
*/
   char        * adist;
   int         * iacol, * iarow, * iccol, * icrow, * ictxt, * lda,
               * ldc, * m, * n, * nb;
   float       * alpha, * beta;
/* ..
*  .. Array Arguments ..
*/
   float       A[], C[], WORK[];
{
/*
*
*  Purpose
*  =======
*
*  PSTRA2 transposes a column block to row block, or a row block to
*  column block by reallocating data distribution:
*
*                         C := alpha*A^T + beta*C,
*
*  where A is an M-by-N matrix  and C is an N-by-M matrix, and the size
*  of M or N is limited to its block size NB.
*
*  The first elements  of the matrices A, and C  should  be  located  at
*  the beginnings of their first blocks. (not the middle of the blocks.)
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   int         i, idex, iend, intv, izero=0, jdex, lcm, lcmp, lcmq, ma,
               mccol, mcrow, mp, mq, mq0, na, np, np0, nq, mycol, myrow,
               npcol, nprow;
   float       one, zero;
/* ..
*  .. External Functions ..
*/
   void        blacs_gridinfo_();
   void        sgerv2d_();
   void        sgesd2d_();
   void        pstr2af();
   void        pstr2at();
   void        pstr2bt();

   F_INTG_FCT  ilcm_();
   F_INTG_FCT  numroc_();
/* ..
*  .. Executable Statements ..
*/
   blacs_gridinfo_( ictxt, &nprow, &npcol, &myrow, &mycol );

   lcm = ilcm_( &nprow, &npcol );
   lcmp = lcm / nprow;
   lcmq = lcm / npcol;
   one = ONE;
   zero = ZERO;

   if( (*adist) == 'C' )
   {
      mp  = numroc_( m, nb, &myrow, iarow, &nprow );
      mq  = numroc_( m, nb, &mycol, iccol, &npcol );
      mq0 = numroc_( m, nb, &izero, &izero, &npcol );
      mq0 = numroc_( &mq0, nb, &izero, &izero, &lcmq );

      intv = lcmp*(*nb);
      iend = CEIL( *m, *nb );
      iend = MIN( lcm, iend ) - 1;
      for( i = 0; i <= iend; i++ )
      {
         mcrow = ( (i % nprow) + (*iarow) ) % nprow;
         mccol = ( (i % npcol) + (*iccol) ) % npcol;
         if( lcmq == 1 )
            mq0 = numroc_( m, nb, &i, &izero, &npcol );
         jdex = ( i / npcol ) * (*nb);
/*
*        A source node copies the blocks to WORK, and send it
*/
         if( ( myrow == mcrow ) && ( mycol == (*iacol) ) )
         {                   /* The source node is a destination node */
            idex = ( i / nprow ) * (*nb);
            if( ( myrow == (*icrow) ) && ( mycol == mccol ) )
            {
               ma = mp - idex;
               pstr2at( ictxt, COLUMN, &ma, n, nb, alpha, &A[idex],
                        lda, beta, &C[jdex*(*ldc)], ldc, &lcmp, &lcmq );
            }
            else
            {       /* Source node sends blocks to a destination node */
               ma = mp - idex;
               pstr2bt( ictxt, COLUMN, &ma, n, nb, alpha, &A[idex], lda,
                        &zero, WORK, n, &intv );
               sgesd2d_( ictxt, n, &mq0, WORK, n, icrow, &mccol );
            }
         }
         else if( ( myrow == (*icrow) ) && ( mycol == mccol ) )
         {           /* A destination node receives the copied blocks */
            if( ( lcmq == 1 ) && ( *beta == ZERO ) )
            {
               sgerv2d_( ictxt, n, &mq0, C, ldc, &mcrow, iacol );
            }
            else
            {
               sgerv2d_( ictxt, n, &mq0, WORK, n, &mcrow, iacol );
               ma = mq - jdex;
               pstr2af( ictxt, ROW, n, &ma, nb, &one, WORK, n, beta,
                        &C[jdex*(*ldc)], ldc, &lcmp, &lcmq, &mq0 );
            }
         }
      }
   }
   else
   {
      np  = numroc_( n, nb, &myrow, icrow, &nprow );
      nq  = numroc_( n, nb, &mycol, iacol, &npcol );
      np0 = numroc_( n, nb, &izero, &izero, &nprow );
      np0 = numroc_( &np0, nb, &izero, &izero, &lcmp );

      intv = lcmq * (*nb);
      iend = CEIL( *n, *nb );
      iend = MIN( lcm, iend ) - 1;
      for( i = 0; i <= iend; i++ )
      {
         mcrow = ( ( i % nprow ) + (*icrow) ) % nprow;
         mccol = ( ( i % npcol ) + (*iacol) ) % npcol;
         if( lcmp == 1 )
            np0 = numroc_( n, nb, &i, &izero, &nprow );
         idex = ( i / nprow ) * (*nb);
/*
*        A source node copies the blocks to WORK, and send it
*/
         if( ( myrow == (*iarow) ) && ( mycol == mccol ) )
         {                   /* The source node is a destination node */
            jdex = (i / npcol) * (*nb);
            if( ( myrow == mcrow ) && ( mycol == (*iccol) ) )
            {
               na = nq - jdex;
               pstr2at( ictxt, ROW, m, &na, nb, alpha, &A[jdex*(*lda)],
                             lda, beta, &C[idex], ldc, &lcmp, &lcmq );
            }
            else
            {       /* Source node sends blocks to a destination node */
               na = nq - jdex;
               pstr2bt( ictxt, ROW, m, &na, nb, alpha, &A[jdex*(*lda)],
                         lda, &zero, WORK, &np0, &intv );
               sgesd2d_( ictxt, &np0, m, WORK, &np0, &mcrow, iccol );
            }
         }
         else if( ( myrow == mcrow ) && ( mycol == (*iccol) ) )
         {           /* A destination node receives the copied blocks */
            if( ( lcmp == 1 ) && ( *beta == ZERO ) )
            {
               sgerv2d_( ictxt, &np0, m, C, ldc, iarow, &mccol );
            }
            else
            {
               sgerv2d_( ictxt, &np0, m, WORK, &np0, iarow, &mccol );
               na = np - idex;
               pstr2af( ictxt, COLUMN, &na, m, nb, &one, WORK, &np0,
                        beta, &C[idex], ldc, &lcmp, &lcmq, &np0 );
            }
         }
      }
   }
}

void pstr2at( ictxt, adist, m, n, nb, alpha, A, lda, beta, B, ldb,
              lcmp, lcmq )
/*
*  .. Scalar Arguments ..
*/
   char        * adist;
   int         * ictxt, * lcmp, * lcmq, * lda, * ldb, * m, * n, * nb;
   float       * alpha, * beta;
/* ..
*  .. Array Arguments ..
*/
   float       A[], B[];
{
/*
*
*  Purpose
*  =======
*
*  B := alpha*A^T + beta*B, B is a transposed scattered block row
*  (or column), copied from a scattered block column (or row) of A.
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   char        trans='T';
   int         ia=1, iendk, inca, incb, ioffa=0, ioffb=0, jntv, k, ma,
               na;
/* ..
*  .. External Functions ..
*/
   F_VOID_FCT  pbsmatadd_();
/* ..
*  .. Executable Statements ..
*/
   if( (*lcmp) == (*lcmq) )
   {
      pbsmatadd_( ictxt, C2F_CHAR( &trans ), n, m, alpha, A, lda,
                  beta, B, ldb );
   }
   else
   {
      if( (*adist) == 'C' )
      {                                        /* A is a column block */
         inca = (*nb) * (*lcmp);
         incb = (*nb) * (*lcmq) * (*ldb);
         iendk = CEIL( (*m), inca );

         for( k = 1; k <= iendk; k++ )
         {
            ma = (*m) - ioffa;
            ma = MIN( ma, (*nb) );

            pbsmatadd_( ictxt, C2F_CHAR( &trans ), n, &ma, alpha,
                        &A[ioffa], lda, beta, &B[ioffb], ldb );

            ioffa += inca;
            ioffb += incb;
         }
      }
      else
      {                                           /* A is a row block */
         jntv = (*lcmq) * (*nb);
         iendk = CEIL( (*n), jntv );
         inca = jntv * (*lda);
         incb = (*lcmp) * (*nb);

         for( k = 1; k <= iendk; k++ )
         {
            na = (*n) - ia + 1;
            na = MIN( na, (*nb) );

            pbsmatadd_( ictxt, C2F_CHAR( &trans ), &na, m, alpha,
                        &A[ioffa], lda, beta, &B[ioffb], ldb );

            ia += jntv;
            ioffa += inca;
            ioffb += incb;
         }
      }
   }
}

void pstr2bt( ictxt, adist, m, n, nb, alpha, A, lda, beta, B, ldb, intv )
/*
*  .. Scalar Arguments ..
*/
   char        * adist;
   int         * ictxt, * intv, * lda, * ldb, * m, * n, * nb;
   float       * alpha, * beta;
/* ..
*  .. Array Arguments ..
*/
   float       A[], B[];
{
/*
*
*  Purpose
*  =======
*
*  B := alpha*A^T + beta*B, B is a transposed condensed block row
*  (or column), copied from a scattered block column (or row) of A.
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   char        trans='T';
   int         ia=1, iendk, inca, incb, ioffa=0, ioffb=0, k, ma, na;
/* ..
*  .. External Functions ..
*/
   F_VOID_FCT  pbsmatadd_();
/* ..
*  .. Executable Statements ..
*/
   if( (*intv) == (*nb) )
   {
      pbsmatadd_( ictxt, C2F_CHAR( &trans ), n, m, alpha, A, lda, beta,
               B, ldb );
   }
   else
   {
      if( (*adist) == 'C' )
      {                                        /* A is a column block */
         iendk = CEIL( (*m), (*intv) );
         inca = (*intv);
         incb = (*nb) * (*ldb);

         for( k = 1; k <= iendk; k++ )
         {
            ma = (*m) - ioffa;
            ma = MIN( ma, (*nb) );

            pbsmatadd_( ictxt, C2F_CHAR( &trans ), n, &ma, alpha,
                        &A[ioffa], lda, beta, &B[ioffb], ldb );

            ioffa += inca;
            ioffb += incb;
         }
      }
      else
      {                                           /* A is a row block */
         iendk = CEIL( (*n), (*intv) );
         inca = (*intv) * (*lda);
         incb = (*nb);

         for( k = 1; k <= iendk; k++ )
         {
            na = (*n) - ia + 1;
            na = MIN( na, (*nb) );

            pbsmatadd_( ictxt, C2F_CHAR( &trans ), &na, m, alpha,
                        &A[ioffa], lda, beta, &B[ioffb], ldb );

            ia += (*intv);
            ioffa += inca;
            ioffb += incb;
         }
      }
   }
}

void pstr2af( ictxt, adist, m, n, nb, alpha, A, lda, beta, B, ldb,
              lcmp, lcmq, nint )
/*
*  .. Scalar Arguments ..
*/
   char        * adist;
   int         * ictxt, * lcmp, * lcmq, * lda, * ldb, * m, * n, * nb,
               * nint;
   float       * alpha, * beta;
/* ..
*  .. Array Arguments ..
*/
   float       A[], B[];
{
/*
*
*  Purpose
*  =======
*
*     B <== alpha * A + beta * B
*     B is a scattered block row (or column)
*       copied from a (condensed) block column (or row) of A
*
*  =====================================================================
*
*  .. Local Scalars ..
*/
   char        full='G';
   int         iendk, inca, incb, intv, ioffa=0, ioffb=0, jb=1, k, ma,
               na;
/* ..
*  .. External Functions ..
*/
   F_VOID_FCT  pbsmatadd_();
/* ..
*  .. Executable Statements ..
*/
   iendk = CEIL( (*nint), (*nb) );

   if( (*adist) == 'R' )
   {
      intv = (*nb) * (*lcmq);
      inca = (*nb) * (*lda);
      incb = intv * (*ldb);

      for( k = 1; k <= iendk; k++ )
      {
         na = (*n) - jb + 1;
         na = MIN( na, (*nb) );

         pbsmatadd_( ictxt, C2F_CHAR( &full ), m, &na, alpha, &A[ioffa],
                     lda, beta, &B[ioffb], ldb );
         jb += intv;
         ioffa += inca;
         ioffb += incb;
      }
   }
   else
   {
      inca = (*nb);
      incb = (*nb) * (*lcmp);

      for( k = 1; k <= iendk; k++ )
      {
         ma = (*m) - ioffb;
         ma = MIN( ma, (*nb) );

         pbsmatadd_( ictxt, C2F_CHAR( &full ), &ma, n, alpha, &A[ioffa],
                     lda, beta, &B[ioffb], ldb );

         ioffa += inca;
         ioffb += incb;
      }
   }
}