PB_Cpsym.c

Go to the documentation of this file.

/* ---------------------------------------------------------------------
*
*  -- PBLAS auxiliary routine (version 2.0) --
*     University of Tennessee, Knoxville, Oak Ridge National Laboratory,
*     and University of California, Berkeley.
*     April 1, 1998
*
*  ---------------------------------------------------------------------
*/
/*
*  Include files
*/
#include "../pblas.h"
#include "../PBpblas.h"
#include "../PBtools.h"
#include "../PBblacs.h"
#include "../PBblas.h"
 
#ifdef __STDC__
void PB_Cpsym( PBTYP_T * TYPE, PBTYP_T * UTYP, char * SIDE,
               char * UPLO, Int N, Int K, char * ALPHA, char * A,
               Int IA, Int JA, Int * DESCA, char * XC, Int LDXC,
               char * XR, Int LDXR, char * YC, Int LDYC, char * YR,
               Int LDYR, TZSYM_T SYM )
#else
void PB_Cpsym( TYPE, UTYP, SIDE, UPLO, N, K, ALPHA, A, IA, JA, DESCA,
               XC, LDXC, XR, LDXR, YC, LDYC, YR, LDYR, SYM )
/*
*  .. Scalar Arguments ..
*/
   char           * SIDE, * UPLO;
   Int            IA, JA, K, LDXC, LDXR, LDYC, LDYR, N;
   char           * ALPHA;
   PBTYP_T        * TYPE, * UTYP;
   TZSYM_T        SYM;
/*
*  .. Array Arguments ..
*/
   Int            * DESCA;
   char           * A, * XC, * XR, * YC, * YR;
#endif
{
/*
*  Purpose
*  =======
*
*  PB_Cpsym  performs  a symmetric or Hermitian matrix-matrix or matrix-
*  vector multiplication.  In the following, sub( A ) denotes the symme-
*  tric or Hermitian submatrix operand A( IA:IA+N-1, JA:JA+N-1 ).
*
*  Notes
*  =====
*
*  A description  vector  is associated with each 2D block-cyclicly dis-
*  tributed matrix.  This  vector  stores  the  information  required to
*  establish the  mapping  between a  matrix entry and its corresponding
*  process and memory location.
*
*  In  the  following  comments,   the character _  should  be  read  as
*  "of  the  distributed  matrix".  Let  A  be a generic term for any 2D
*  block cyclicly distributed matrix.  Its description vector is DESC_A:
*
*  NOTATION         STORED IN       EXPLANATION
*  ---------------- --------------- ------------------------------------
*  DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
*  CTXT_A  (global) DESCA[ CTXT_  ] The BLACS context handle, indicating
*                                   the NPROW x NPCOL BLACS process grid
*                                   A  is  distributed over. The context
*                                   itself  is  global,  but  the handle
*                                   (the integer value) may vary.
*  M_A     (global) DESCA[ M_     ] The  number of rows in the distribu-
*                                   ted matrix A, M_A >= 0.
*  N_A     (global) DESCA[ N_     ] The number of columns in the distri-
*                                   buted matrix A, N_A >= 0.
*  IMB_A   (global) DESCA[ IMB_   ] The number of rows of the upper left
*                                   block of the matrix A, IMB_A > 0.
*  INB_A   (global) DESCA[ INB_   ] The  number  of columns of the upper
*                                   left   block   of   the  matrix   A,
*                                   INB_A > 0.
*  MB_A    (global) DESCA[ MB_    ] The blocking factor used to  distri-
*                                   bute the last  M_A-IMB_A  rows of A,
*                                   MB_A > 0.
*  NB_A    (global) DESCA[ NB_    ] The blocking factor used to  distri-
*                                   bute the last  N_A-INB_A  columns of
*                                   A, NB_A > 0.
*  RSRC_A  (global) DESCA[ RSRC_  ] The process row over which the first
*                                   row of the matrix  A is distributed,
*                                   NPROW > RSRC_A >= 0.
*  CSRC_A  (global) DESCA[ CSRC_  ] The  process column  over  which the
*                                   first column of  A  is  distributed.
*                                   NPCOL > CSRC_A >= 0.
*  LLD_A   (local)  DESCA[ LLD_   ] The  leading dimension  of the local
*                                   array  storing  the  local blocks of
*                                   the distributed matrix A,
*                                   IF( Lc( 1, N_A ) > 0 )
*                                      LLD_A >= MAX( 1, Lr( 1, M_A ) )
*                                   ELSE
*                                      LLD_A >= 1.
*
*  Let K be the number of  rows of a matrix A starting at the global in-
*  dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
*  that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
*  receive if these K rows were distributed over NPROW processes.  If  K
*  is the number of columns of a matrix  A  starting at the global index
*  JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number  of co-
*  lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would  receive if
*  these K columns were distributed over NPCOL processes.
*
*  The values of Lr() and Lc() may be determined via a call to the func-
*  tion PB_Cnumroc:
*  Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
*  Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
*
*  Arguments
*  =========
*
*  TYPE    (local input) pointer to a PBTYP_T structure
*          On entry,  TYPE  is a pointer to a structure of type PBTYP_T,
*          that contains type information (See pblas.h).
*
*  UTYP    (local input) pointer to a PBTYP_T structure
*          On entry,  UTYP  is a pointer to a structure of type PBTYP_T,
*          that contains type information for the Y's (See pblas.h).
*
*  SIDE    (global input) pointer to CHAR
*          On entry,  SIDE  specifies whether  op( sub( A ) ) multiplies
*          its operand X from the left or right as follows:
*
*          SIDE = 'L' or 'l'  Y := alpha*op( sub( A ) )*X + Y,
*
*          SIDE = 'R' or 'r'  Y := alpha*X*op( sub( A ) ) + Y.
*
*  UPLO    (global input) pointer to CHAR
*          On  entry,   UPLO  specifies  whether  the  local  pieces  of
*          the array  A  containing the  upper or lower triangular  part
*          of the symmetric or Hermitian submatrix  sub( A )  are to  be
*          referenced as follows:
*
*             UPLO = 'U' or 'u'   Only the local pieces corresponding to
*                                 the upper triangular part of  the sym-
*                                 metric or Hermitian submatrix sub( A )
*                                 are to be referenced,
*
*             UPLO = 'L' or 'l'   Only the local pieces corresponding to
*                                 the  lower triangular part of the sym-
*                                 metric or Hermitian submatrix sub( A )
*                                 are to be referenced.
*
*  N       (global input) INTEGER
*          On entry,  N specifies the order of the  submatrix  sub( A ).
*          N must be at least zero.
*
*  K       (global input) INTEGER
*          On entry, K  specifies the local number of columns of the lo-
*          cal array XC  and the local number of rows of the local array
*          XR. K mut be at least zero.
*
*  ALPHA   (global input) pointer to CHAR
*          On entry, ALPHA specifies the scalar alpha.
*
*  A       (local input) pointer to CHAR
*          On entry, A is an array of dimension (LLD_A, Ka), where Ka is
*          at least Lc( 1, JA+N-1 ).  Before  entry, this array contains
*          the local entries of the matrix A.
*          Before  entry  with  UPLO = 'U' or 'u', this  array  contains
*          the local entries corresponding to the upper triangular  part
*          of the  @(syhec)  submatrix  sub( A ), and the local entries
*          corresponding to the  strictly lower triangular  of  sub( A )
*          are not  referenced.
*          Before  entry  with  UPLO = 'L' or 'l', this  array  contains
*          the local entries corresponding to the lower triangular  part
*          of the  @(syhec)  submatrix  sub( A ), and the local entries
*          corresponding to the  strictly upper triangular  of  sub( A )
*          are not  referenced.
*
*  IA      (global input) INTEGER
*          On entry, IA  specifies A's global row index, which points to
*          the beginning of the submatrix sub( A ).
*
*  JA      (global input) INTEGER
*          On entry, JA  specifies A's global column index, which points
*          to the beginning of the submatrix sub( A ).
*
*  DESCA   (global and local input) INTEGER array
*          On entry, DESCA  is an integer array of dimension DLEN_. This
*          is the array descriptor for the matrix A.
*
*  XC      (local input) pointer to CHAR
*          On entry, XC is an array of dimension (LDXC,K). Before entry,
*          this array contains the local entries of the matrix XC.
*
*  LDXC    (local input) INTEGER
*          On entry, LDXC  specifies  the leading dimension of the array
*          XC. LDXC must be at least max( 1, Lp( IA, N ) ).
*
*  XR      (local input) pointer to CHAR
*          On entry, XR is an array of dimension (LDXR,Kx),  where Kx is
*          at least Lc( JA, N ). Before  entry, this array contains  the
*          local entries of the matrix XR.
*
*  LDXR    (local input) INTEGER
*          On entry, LDXR  specifies  the leading dimension of the array
*          XR. LDXR must be at least max( 1, K ).
*
*  YC      (local input/local output) pointer to CHAR
*          On entry, YC is an array of dimension (LDYC,K). Before entry,
*          this  array  contains  the local entries of the matrix YC. On
*          exit, this array contains the updated vector YC.
*
*  LDYC    (local input) INTEGER
*          On entry, LDYC  specifies  the leading dimension of the array
*          YC. LDYC must be at least max( 1, Lp( IA, N ) ).
*
*  YR      (local input/local output) pointer to CHAR
*          On entry, YR is an array of dimension (LDYR,Ky),  where Ky is
*          at least Lc( JA, N ).  Before  entry, this array contains the
*          local entries of the matrix YR. On  exit, this array contains
*          the updated vector YR.
*
*  LDYR    (local input) INTEGER
*          On entry, LDYR  specifies  the leading dimension of the array
*          YR. LDYR must be at least max( 1, K ).
*
*  SYM     (local input) pointer to function of type TZSYM_T
*          On entry, SYM specifies the function performing the symmetric
*          or  Hermitian  matrix-vector  or  matrix-matrix multiply of a
*          single block.
*
*  -- Written on April 1, 1998 by
*     Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
*
*  ---------------------------------------------------------------------
*/
/*
*  .. Local Scalars ..
*/
   Int            Acol, Arow, Aii, Aimb1, Ainb1, Ajj, Ald, Amp, Amb, Anb, Anq,
                  Aoffi, Aoffj, Arcol, Arrow, GoEast, GoSouth, IsColRepl,
                  IsRowRepl, XCinc, XRinc, Xii=0, Xjj=0, Xoffi=-1, Xoffj=-1,
                  YCinc, YRinc, iimax, ilow, imbloc, inbloc, ioffd, ioffx, iupp,
                  jjmax, joffd, joffx, lcmt, lcmt00, lmbloc, lnbloc, low, lower,
                  m1, mbloc, mblkd, mblks, mycol, myrow, n1, nbloc, nblkd,
                  nblks, npcol, nprow, pmb, qnb, size, tmp1, upp, upper;
/* ..
*  .. Executable Statements ..
*
*/
/*
*  Quick return if possible
*/
   if( N <= 0 ) return;
/*
*  Retrieve process grid information
*/
   Cblacs_gridinfo( DESCA[CTXT_], &nprow, &npcol, &myrow, &mycol );
/*
*  Retrieve sub( A )'s local information: Aii, Ajj, Arow, Acol ...
*/
   PB_Cainfog2l( N, N, IA, JA, DESCA, nprow, npcol, myrow, mycol, &Aimb1,
                 &Ainb1, &Amp, &Anq, &Aii, &Ajj, &Arow, &Acol, &Arrow, &Arcol );
/*
*  Quick return if I don't own any of sub( A ) or if sub( A ) is replicated in
*  all processes.
*/
   if( ( Amp <= 0 ) || ( Anq <= 0 ) ) return;
 
   IsRowRepl = ( ( Arow < 0 ) || ( nprow == 1 ) );
   IsColRepl = ( ( Acol < 0 ) || ( npcol == 1 ) );
   Amb  = DESCA[ MB_ ]; Anb = DESCA[ NB_ ]; Ald = DESCA[ LLD_ ];
   size = TYPE->size;
 
   if( IsRowRepl && IsColRepl )
   {
      SYM( TYPE, SIDE, UPLO, Amp, Anq, K, 0, ALPHA,  Mptr( A, Aii, Ajj, Ald,
           size ), Ald, XC, LDXC, XR, LDXR, YC, LDYC, YR, LDYR );
      return;
   }
 
   XCinc = size;         XRinc = LDXR * size;
   YCinc = UTYP->size;   YRinc = LDYR * UTYP->size;
   upper = ( Mupcase( UPLO[0] ) == CUPPER );
   lower = ( Mupcase( UPLO[0] ) == CLOWER );
/*
*  Initialize lcmt00, mblks, nblks, imbloc, inbloc, lmbloc, lnbloc, ilow, low,
*  iupp, and upp.
*/
   PB_Cbinfo( 0, Amp, Anq, Aimb1, Ainb1, Amb, Anb, Arrow, Arcol, &lcmt00,
              &mblks, &nblks, &imbloc, &inbloc, &lmbloc, &lnbloc, &ilow, &low,
              &iupp, &upp );
 
   iimax = ( Aoffi = Aii - 1 ) + ( m1 = Amp );
   jjmax = ( Aoffj = Ajj - 1 ) + ( n1 = Anq );
   pmb   = ( IsRowRepl ? Amb : nprow * Amb );
   qnb   = ( IsColRepl ? Anb : npcol * Anb );
/*
*  Handle separately the first row and/or column of the LCM table. Update the
*  LCM value of the curent block lcmt00, as well as the number of rows and
*  columns mblks and nblks remaining in the LCM table.
*/
   GoSouth = ( lcmt00 > iupp );
   GoEast  = ( lcmt00 < ilow );
/*
*  Go through the table looking for blocks owning diagonal entries.
*/
   if( ( !( GoSouth ) ) && ( !( GoEast ) ) )
   {
/*
*  The upper left block owns diagonal entries lcmt00 >= ilow && lcmt00 <= iupp
*/
      SYM( TYPE, SIDE, UPLO, imbloc, inbloc, K, lcmt00, ALPHA, Mptr( A, Aii,
           Ajj, Ald, size ), Ald, XC+Xii*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
           YC+Xii*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
/*
*  Decide whether one should go south or east in the table: Go east if the
*  block below the current one only owns lower entries. If this block, however,
*  owns diagonals, then go south.
*/
      GoSouth = !( GoEast = ( ( lcmt00 - ( iupp - upp + pmb ) ) < ilow ) );
 
      if( GoSouth )
      {
/*
*  When the upper triangular part of sub( A ) should be operated with and
*  one is planning to go south in the table, it is neccessary to take care
*  of the remaining columns of these imbloc rows immediately.
*/
         if( upper && ( Anq > inbloc ) )
         {
            tmp1 = Anq - inbloc;
            SYM( TYPE, SIDE, ALL, imbloc, tmp1, K, 0,
                 ALPHA, Mptr( A, Aii, Ajj+inbloc, Ald, size ), Ald,
                 XC+Xii*XCinc, LDXC, XR+(Xjj+inbloc)*XRinc, LDXR,
                 YC+Xii*YCinc, LDYC, YR+(Xjj+inbloc)*YRinc, LDYR );
         }
         Aii += imbloc; Xii += imbloc; m1  -= imbloc;
      }
      else
      {
/*
*  When the lower triangular part of sub( A ) should be operated with and
*  one is planning to go east in the table, it is neccessary to take care
*  of the remaining rows of these inbloc columns immediately.
*/
         if( lower && ( Amp > imbloc ) )
         {
            tmp1 = Amp - imbloc;
            SYM( TYPE, SIDE, ALL, tmp1, inbloc, K, 0,
                 ALPHA, Mptr( A, Aii+imbloc, Ajj, Ald, size ), Ald,
                 XC+(Xii+imbloc)*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
                 YC+(Xii+imbloc)*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
         }
         Ajj += inbloc; Xjj += inbloc; n1  -= inbloc;
      }
   }
 
   if( GoSouth )
   {
/*
*  Go one step south in the LCM table. Adjust the current LCM value as well as
*  the local row indexes in A and XC.
*/
      lcmt00 -= ( iupp - upp + pmb ); mblks--;
      Aoffi  += imbloc; Xoffi  += imbloc;
/*
*  While there are blocks remaining that own upper entries, keep going south.
*  Adjust the current LCM value as well as the local row indexes in A and XC.
*/
      while( ( mblks > 0 ) && ( lcmt00 > upp ) )
      { lcmt00 -= pmb; mblks--; Aoffi += Amb; Xoffi += Amb; }
/*
*  Operate with the upper triangular part of sub( A ) we just skipped when
*  necessary.
*/
      tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
      if( upper && ( tmp1 > 0 ) )
      {
         SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
              ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
         Aii += tmp1; Xii += tmp1; m1  -= tmp1;
      }
/*
*  Return if no more row in the LCM table.
*/
      if( mblks <= 0 ) return;
/*
*  lcmt00 <= upp. The current block owns either diagonals or lower entries.
*  Save the current position in the LCM table. After this column has been
*  completely taken care of, re-start from this row and the next column of
*  the LCM table.
*/
      lcmt = lcmt00; mblkd = mblks; ioffd = Aoffi; ioffx = Xoffi;
 
      mbloc = Amb;
      while( ( mblkd > 0 ) && ( lcmt >= ilow ) )
      {
/*
*  A block owning diagonals lcmt00 >= ilow && lcmt00 <= upp has been found.
*/
         if( mblkd == 1 ) mbloc = lmbloc;
         SYM( TYPE, SIDE, UPLO, mbloc, inbloc, K, lcmt,
              ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
              XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
         lcmt00 = lcmt;  lcmt  -= pmb;
         mblks  = mblkd; mblkd--;
         Aoffi  = ioffd; ioffd += mbloc;
         Xoffi  = ioffx; ioffx += mbloc;
      }
/*
*  Operate with the lower triangular part of sub( A ).
*/
      tmp1 = m1 - ioffd + Aii - 1;
      if( lower && ( tmp1 > 0 ) )
         SYM( TYPE, SIDE, ALL, tmp1, inbloc, K, 0,
              ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
              XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
 
      tmp1    = Aoffi - Aii + 1;
      m1     -= tmp1;
      n1     -= inbloc;
      lcmt00 += low - ilow + qnb;
      nblks--;
      Aoffj  += inbloc;
      Xoffj  += inbloc;
/*
*  Operate with the upper triangular part of sub( A ).
*/
      if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
         SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
              ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
      Aii = Aoffi + 1; Ajj = Aoffj + 1;
      Xii = Xoffi + 1; Xjj = Xoffj + 1;
   }
   else if( GoEast )
   {
/*
*  Go one step east in the LCM table. Adjust the current LCM value as well as
*  the local column index in A and XR.
*/
      lcmt00 += low - ilow + qnb; nblks--;
      Aoffj  += inbloc; Xoffj  += inbloc;
/*
*  While there are blocks remaining that own lower entries, keep going east.
*  Adjust the current LCM value as well as the local column index in A and XR.
*/
      while( ( nblks > 0 ) && ( lcmt00 < low ) )
      { lcmt00 += qnb; nblks--; Aoffj += Anb; Xoffj += Anb; }
/*
*  Operate with the lower triangular part of sub( A ).
*/
      tmp1 = MIN( Aoffj, jjmax ) - Ajj + 1;
      if( lower && ( tmp1 > 0 ) )
      {
         SYM( TYPE, SIDE, ALL, m1, tmp1, K, 0,
              ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
         Ajj += tmp1; Xjj += tmp1; n1  -= tmp1;
      }
/*
*  Return if no more column in the LCM table.
*/
      if( nblks <= 0 ) return;
/*
*  lcmt00 >= low. The current block owns either diagonals or upper entries.
*  Save the current position in the LCM table. After this row has been
*  completely taken care of, re-start from this column and the next row of
*  the LCM table.
*/
      lcmt = lcmt00; nblkd = nblks; joffd = Aoffj; joffx = Xoffj;
 
      nbloc = Anb;
      while( ( nblkd > 0 ) && ( lcmt <= iupp ) )
      {
/*
*  A block owning diagonals lcmt00 >= low && lcmt00 <= iupp has been found.
*/
         if( nblkd == 1 ) nbloc = lnbloc;
         SYM( TYPE, SIDE, UPLO, imbloc, nbloc, K, lcmt,
              ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+(joffx+1)*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+(joffx+1)*YRinc, LDYR );
         lcmt00 = lcmt;  lcmt  += qnb;
         nblks  = nblkd; nblkd--;
         Aoffj  = joffd; joffd += nbloc;
         Xoffj  = joffx; joffx += nbloc;
      }
/*
*  Operate with the upper triangular part of sub( A ).
*/
      tmp1 = n1 - joffd + Ajj - 1;
      if( upper && ( tmp1 > 0 ) )
         SYM( TYPE, SIDE, ALL, imbloc, tmp1, K, 0,
              ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+(joffx+1)*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+(joffx+1)*YRinc, LDYR );
 
      tmp1    = Aoffj - Ajj + 1;
      m1     -= imbloc;
      n1     -= tmp1;
      lcmt00 -= ( iupp - upp + pmb );
      mblks--;
      Aoffi  += imbloc;
      Xoffi  += imbloc;
/*
*  Operate with the lower triangular part of sub( A ).
*/
      if( lower && ( m1 > 0 ) && ( tmp1 > 0 ) )
         SYM( TYPE, SIDE, ALL, m1, tmp1, K, 0,
              ALPHA, Mptr( A, Aoffi+1, Ajj, Ald, size ), Ald,
              XC+(Xoffi+1)*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
              YC+(Xoffi+1)*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
      Aii = Aoffi + 1; Ajj = Aoffj + 1;
      Xii = Xoffi + 1; Xjj = Xoffj + 1;
   }
/*
*  Loop over the remaining columns of the LCM table.
*/
   nbloc = Anb;
   while( nblks > 0 )
   {
      if( nblks == 1 ) nbloc = lnbloc;
/*
*  While there are blocks remaining that own upper entries, keep going south.
*  Adjust the current LCM value as well as the local row index in A and XC.
*/
      while( ( mblks > 0 ) && ( lcmt00 > upp ) )
      { lcmt00 -= pmb; mblks--; Aoffi += Amb; Xoffi += Amb; }
/*
*  Operate with the upper triangular part of sub( A ).
*/
      tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
      if( upper && ( tmp1 > 0 ) )
      {
         SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
              ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
         Aii += tmp1;
         Xii += tmp1;
         m1  -= tmp1;
      }
/*
*  Return if no more row in the LCM table.
*/
      if( mblks <= 0 ) return;
/*
*  lcmt00 <= upp. The current block owns either diagonals or lower entries.
*  Save the current position in the LCM table. After this column has been
*  completely taken care of, re-start from this row and the next column of
*  the LCM table.
*/
      lcmt  = lcmt00; mblkd = mblks; ioffd = Aoffi; ioffx = Xoffi;
 
      mbloc = Amb;
      while( ( mblkd > 0 ) && ( lcmt >= low ) )
      {
/*
*  A block owning diagonals lcmt00 >= low && lcmt00 <= upp has been found.
*/
         if( mblkd == 1 ) mbloc = lmbloc;
         SYM( TYPE, SIDE, UPLO, mbloc, nbloc, K, lcmt,
              ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
              XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
         lcmt00 = lcmt;  lcmt  -= pmb;
         mblks  = mblkd; mblkd--;
         Aoffi  = ioffd; Xoffi  = ioffx;
         ioffd += mbloc; ioffx += mbloc;
      }
/*
*  Operate with the lower triangular part of sub( A ).
*/
      tmp1 = m1 - ioffd + Aii - 1;
      if( lower && ( tmp1 > 0 ) )
         SYM( TYPE, SIDE, ALL, tmp1, nbloc, K, 0,
              ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
              XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
 
      tmp1    = MIN( Aoffi, iimax ) - Aii + 1;
      m1     -= tmp1;
      n1     -= nbloc;
      lcmt00 += qnb;
      nblks--;
      Aoffj  += nbloc;
      Xoffj  += nbloc;
/*
*  Operate with the upper triangular part of sub( A ).
*/
      if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
         SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
              ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
              XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
              YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
      Aii = Aoffi + 1; Ajj = Aoffj + 1;
      Xii = Xoffi + 1; Xjj = Xoffj + 1;
   }
/*
*  End of PB_Cpsym
*/
}