|
ScaLAPACK
2.0.2
ScaLAPACK: Scalable Linear Algebra PACKage
|
#include "../pblas.h"#include "../PBpblas.h"#include "../PBtools.h"#include "../PBblacs.h"#include "../PBblas.h"Go to the source code of this file.
Functions/Subroutines | |
| 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) |
| 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 | ||
| ) |
Definition at line 26 of file PB_Cpsym.c.
{
/*
* 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
*/
}