/* --------------------------------------------------------------------- * * -- 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_Cptrm( PBTYP_T * TYPE, PBTYP_T * UTYP, char * SIDE, char * UPLO, char * TRANS, char * DIAG, int N, int K, char * ALPHA, char * A, int IA, int JA, int * DESCA, char * X, int LDX, char * Y, int LDY, TZTRM_T TRM ) #else void PB_Cptrm( TYPE, UTYP, SIDE, UPLO, TRANS, DIAG, N, K, ALPHA, A, IA, JA, DESCA, X, LDX, Y, LDY, TRM ) /* * .. Scalar Arguments .. */ char * DIAG, * SIDE, * TRANS, * UPLO; int IA, JA, K, LDX, LDY, N; char * ALPHA; PBTYP_T * TYPE, * UTYP; TZTRM_T TRM; /* * .. Array Arguments .. */ int * DESCA; char * A, * X, * Y; #endif { /* * Purpose * ======= * * PB_Cptrm performs a triangular matrix-matrix or matrix-vector multi- * plication. In the following, sub( A ) denotes the triangular subma- * trix 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 submatrix sub( A ) is * an upper or lower triangular submatrix as follows: * * UPLO = 'U' or 'u' sub( A ) is an upper triangular * submatrix, * * UPLO = 'L' or 'l' sub( A ) is a lower triangular * submatrix. * * TRANS (global input) pointer to CHAR * On entry, TRANS specifies the operation to be performed as * follows: * * TRANS = 'N' or 'n' Y := alpha * sub( A ) * X + Y, * * TRANS = 'T' or 't' Y := alpha * sub( A )' * X + Y, * * TRANS = 'C' or 'c' Y := alpha * sub( A )' * X + Y, or * Y := alpha * conjg(sub( A )') * X + Y. * * DIAG (global input) pointer to CHAR * On entry, DIAG specifies whether or not sub( A ) is unit * triangular as follows: * * DIAG = 'U' or 'u' sub( A ) is assumed to be unit trian- * gular, * * DIAG = 'N' or 'n' sub( A ) is not assumed to be unit tri- * angular. * * 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 X and the local number of rows of the local array * Y when SIDE is 'L' or 'l' and TRANS is 'N' or 'n', or SIDE is * 'R' or 'r' and TRANS is 'T', 't', 'C' or 'c'. Otherwise, K * specifies the local number of rows of the local array X and * the local number of columns of the local array Y. 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 triangular 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 triangular submatrix sub( A ), and the local entries * corresponding to the strictly upper triangular of sub( A ) * are not referenced. * Note that when DIAG = 'U' or 'u', the local entries corres- * ponding to the diagonal elements of the submatrix sub( A ) * are not referenced either, but are assumed to be unity. * * 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. * * X (local input) pointer to CHAR * On entry, X is an array of dimension (LDX,Kx), where Kx is * at least Lc( JA, N ) when SIDE is 'L' or 'l' and TRANS is 'N' * or 'n', or SIDE is 'R' or 'r' and TRANS is 'T', 't', 'C' or * 'c', and K otherwise. Before entry, this array contains the * local entries of the matrix X. * * LDX (local input) INTEGER * On entry, LDX specifies the leading dimension of the array X. * LDX must be at least K when SIDE is 'L' or 'l' and TRANS is * 'N' or 'n', or SIDE is 'R' or 'r' and TRANS is 'T', 't', 'C' * or 'c', and max( 1, Lp( IA, N ) ) otherwise. * * Y (local input/local output) pointer to CHAR * On entry, Y is an array of dimension ( LDY, Ky ), where Ky is * at least max( 1, K ) when SIDE is 'L' or 'l' and TRANS is 'N' * or 'n', or SIDE is 'R' or 'r' and TRANS is 'T', 't', 'C' or * 'c', and max( 1, Lc( JA, N ) ) otherwise. Before entry, this * array contains the local entries of the matrix Y. On exit, * this array contains the updated vector Y. * * LDY (local input) INTEGER * On entry, LDY specifies the leading dimension of the array Y. * LDY must be at least max( 1, Lr( IA, N ) ) when SIDE is 'L' * or 'l' and TRANS is 'N' or 'n', or SIDE is 'R' or 'r' and * TRANS is 'T', 't', 'C' or 'c', and max( 1, K ) otherwise. * * TRM (local input) pointer to function of type TZTRM_T * On entry, TRM specifies the function performing the triangu- * lar multiplication 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, Xinc, Yinc, XYii=0, XYjj=0, XYoffi=-1, XYoffj=-1, XisRow, iimax, ilow, imbloc, inbloc, ioffd, ioffxy, iupp, jjmax, joffd, joffxy, 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 ) { TRM( TYPE, SIDE, UPLO, TRANS, DIAG, Amp, Anq, K, 0, ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald, X, LDX, Y, LDY ); return; } if( Mupcase( SIDE[0] ) == CLEFT ) { if( Mupcase( TRANS[0] ) == CNOTRAN ) { XisRow = 1; Xinc = LDX * size; Yinc = UTYP->size; } else { XisRow = 0; Xinc = size; Yinc = LDY * UTYP->size; } } else { if( Mupcase( TRANS[0] ) == CNOTRAN ) { XisRow = 0; Xinc = size; Yinc = LDY * UTYP->size; } else { XisRow = 1; Xinc = LDX * size; Yinc = 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 ); if( XisRow ) { /* * Go through the table looking for blocks owning diagonal entries. */ if( ( !( GoSouth ) ) && ( !( GoEast ) ) ) { /* * The upper left block owns diagonal entries lcmt00 >= ilow && lcmt00 <= iupp */ TRM( TYPE, SIDE, UPLO, TRANS, DIAG, imbloc, inbloc, K, lcmt00, ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald, X+XYjj*Xinc, LDX, Y+XYii*Yinc, LDY ); /* * 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; TRM( TYPE, SIDE, ALL, TRANS, DIAG, imbloc, tmp1, K, 0, ALPHA, Mptr( A, Aii, Ajj+inbloc, Ald, size ), Ald, X+(XYjj+inbloc)*Xinc, LDX, Y+XYii*Yinc, LDY ); } Aii += imbloc; XYii += 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; TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, inbloc, K, 0, ALPHA, Mptr( A, Aii+imbloc, Ajj, Ald, size ), Ald, X+XYjj*Xinc, LDX, Y+(XYii+imbloc)*Yinc, LDY ); } Ajj += inbloc; XYjj += 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; XYoffi += 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; XYoffi += 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 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+XYii*Yinc, LDY ); Aii += tmp1; XYii += 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; ioffxy = XYoffi; mbloc = Amb; while( ( mblkd > 0 ) && ( lcmt >= ilow ) ) { /* * A block owning diagonals lcmt00 >= ilow && lcmt00 <= upp has been found. */ if( mblkd == 1 ) mbloc = lmbloc; TRM( TYPE, SIDE, UPLO, TRANS, DIAG, mbloc, inbloc, K, lcmt, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+(ioffxy+1)*Yinc, LDY ); lcmt00 = lcmt; lcmt -= pmb; mblks = mblkd; mblkd--; Aoffi = ioffd; XYoffi = ioffxy; ioffd += mbloc; ioffxy += mbloc; } /* * Operate with the lower triangular part of sub( A ). */ tmp1 = m1 - ioffd + Aii - 1; if( lower && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, inbloc, K, 0, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+(ioffxy+1)*Yinc, LDY ); } tmp1 = Aoffi - Aii + 1; m1 -= tmp1; n1 -= inbloc; lcmt00 += low - ilow + qnb; nblks--; Aoffj += inbloc; XYoffj += inbloc; /* * Operate with the upper triangular part of sub( A ). */ if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+XYii*Yinc, LDY ); } Aii = Aoffi + 1; Ajj = Aoffj + 1; XYii = XYoffi + 1; XYjj = XYoffj + 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; XYoffj += 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; XYoffj += Anb; } /* * Operate with the lower triangular part of sub( A ). */ tmp1 = MIN( Aoffj, jjmax ) - Ajj + 1; if( lower && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, m1, tmp1, K, 0, ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald, X+XYjj*Xinc, LDX, Y+XYii*Yinc, LDY ); Ajj += tmp1; XYjj += 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; joffxy = XYoffj; nbloc = Anb; while( ( nblkd > 0 ) && ( lcmt <= iupp ) ) { /* * A block owning diagonals lcmt00 >= low && lcmt00 <= iupp has been found. */ if( nblkd == 1 ) nbloc = lnbloc; TRM( TYPE, SIDE, UPLO, TRANS, DIAG, imbloc, nbloc, K, lcmt, ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald, X+(joffxy+1)*Xinc, LDX, Y+XYii*Yinc, LDY ); lcmt00 = lcmt; lcmt += qnb; nblks = nblkd; nblkd--; Aoffj = joffd; XYoffj = joffxy; joffd += nbloc; joffxy += nbloc; } /* * Operate with the upper triangular part of sub( A ). */ tmp1 = n1 - joffd + Ajj - 1; if( upper && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, imbloc, tmp1, K, 0, ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald, X+(joffxy+1)*Xinc, LDX, Y+XYii*Yinc, LDY ); } tmp1 = Aoffj - Ajj + 1; m1 -= imbloc; n1 -= tmp1; lcmt00 -= ( iupp - upp + pmb ); mblks--; Aoffi += imbloc; XYoffi += imbloc; /* * Operate with the lower triangular part of sub( A ). */ if( lower && ( m1 > 0 ) && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, m1, tmp1, K, 0, ALPHA, Mptr( A, Aoffi+1, Ajj, Ald, size ), Ald, X+XYjj*Xinc, LDX, Y+(XYoffi+1)*Yinc, LDY ); } Aii = Aoffi + 1; Ajj = Aoffj + 1; XYii = XYoffi + 1; XYjj = XYoffj + 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; XYoffi += Amb; } /* * Operate with the upper triangular part of sub( A ). */ tmp1 = MIN( Aoffi, iimax ) - Aii + 1; if( upper && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+XYii*Yinc, LDY ); Aii += tmp1; XYii += 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; ioffxy = XYoffi; mbloc = Amb; while( ( mblkd > 0 ) && ( lcmt >= low ) ) { /* * A block owning diagonals lcmt00 >= low && lcmt00 <= upp has been found. */ if( mblkd == 1 ) mbloc = lmbloc; TRM( TYPE, SIDE, UPLO, TRANS, DIAG, mbloc, nbloc, K, lcmt, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+(ioffxy+1)*Yinc, LDY ); lcmt00 = lcmt; lcmt -= pmb; mblks = mblkd; mblkd--; Aoffi = ioffd; XYoffi = ioffxy; ioffd += mbloc; ioffxy += mbloc; } /* * Operate with the lower triangular part of sub( A ). */ tmp1 = m1 - ioffd + Aii - 1; if( lower && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, nbloc, K, 0, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+(ioffxy+1)*Yinc, LDY ); } tmp1 = MIN( Aoffi, iimax ) - Aii + 1; m1 -= tmp1; n1 -= nbloc; lcmt00 += qnb; nblks--; Aoffj += nbloc; XYoffj += nbloc; /* * Operate with the upper triangular part of sub( A ). */ if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+(XYoffj+1)*Xinc, LDX, Y+XYii*Yinc, LDY ); } Aii = Aoffi + 1; Ajj = Aoffj + 1; XYii = XYoffi + 1; XYjj = XYoffj + 1; } } else { /* * Go through the table looking for blocks owning diagonal entries. */ if( ( !( GoSouth ) ) && ( !( GoEast ) ) ) { /* * The upper left block owns diagonal entries lcmt00 >= ilow && lcmt00 <= iupp */ TRM( TYPE, SIDE, UPLO, TRANS, DIAG, imbloc, inbloc, K, lcmt00, ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+XYjj*Yinc, LDY ); /* * 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; TRM( TYPE, SIDE, ALL, TRANS, DIAG, imbloc, tmp1, K, 0, ALPHA, Mptr( A, Aii, Ajj+inbloc, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(XYjj+inbloc)*Yinc, LDY ); } Aii += imbloc; XYii += 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; TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, inbloc, K, 0, ALPHA, Mptr( A, Aii+imbloc, Ajj, Ald, size ), Ald, X+(XYii+imbloc)*Xinc, LDX, Y+XYjj*Yinc, LDY ); } Ajj += inbloc; XYjj += 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; XYoffi += 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; XYoffi += 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 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); Aii += tmp1; XYii += 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; ioffxy = XYoffi; mbloc = Amb; while( ( mblkd > 0 ) && ( lcmt >= ilow ) ) { /* * A block owning diagonals lcmt00 >= ilow && lcmt00 <= upp has been found. */ if( mblkd == 1 ) mbloc = lmbloc; TRM( TYPE, SIDE, UPLO, TRANS, DIAG, mbloc, inbloc, K, lcmt, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(ioffxy+1)*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); lcmt00 = lcmt; lcmt -= pmb; mblks = mblkd; mblkd--; Aoffi = ioffd; XYoffi = ioffxy; ioffd += mbloc; ioffxy += mbloc; } /* * Operate with the lower triangular part of sub( A ). */ tmp1 = m1 - ioffd + Aii - 1; if( lower && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, inbloc, K, 0, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(ioffxy+1)*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); } tmp1 = Aoffi - Aii + 1; m1 -= tmp1; n1 -= inbloc; lcmt00 += low - ilow + qnb; nblks--; Aoffj += inbloc; XYoffj += inbloc; /* * Operate with the upper triangular part of sub( A ). */ if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); } Aii = Aoffi + 1; Ajj = Aoffj + 1; XYii = XYoffi + 1; XYjj = XYoffj + 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; XYoffj += 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; XYoffj += Anb; } /* * Operate with the lower triangular part of sub( A ). */ tmp1 = MIN( Aoffj, jjmax ) - Ajj + 1; if( lower && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, m1, tmp1, K, 0, ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+XYjj*Yinc, LDY ); Ajj += tmp1; XYjj += 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; joffxy = XYoffj; nbloc = Anb; while( ( nblkd > 0 ) && ( lcmt <= iupp ) ) { /* * A block owning diagonals lcmt00 >= low && lcmt00 <= iupp has been found. */ if( nblkd == 1 ) nbloc = lnbloc; TRM( TYPE, SIDE, UPLO, TRANS, DIAG, imbloc, nbloc, K, lcmt, ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(joffxy+1)*Yinc, LDY ); lcmt00 = lcmt; lcmt += qnb; nblks = nblkd; nblkd--; Aoffj = joffd; XYoffj = joffxy; joffd += nbloc; joffxy += nbloc; } /* * Operate with the upper triangular part of sub( A ). */ tmp1 = n1 - joffd + Ajj - 1; if( upper && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, imbloc, tmp1, K, 0, ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(joffxy+1)*Yinc, LDY ); } tmp1 = Aoffj - Ajj + 1; m1 -= imbloc; n1 -= tmp1; lcmt00 -= ( iupp - upp + pmb ); mblks--; Aoffi += imbloc; XYoffi += imbloc; /* * Operate with the lower triangular part of sub( A ). */ if( lower && ( m1 > 0 ) && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, m1, tmp1, K, 0, ALPHA, Mptr( A, Aoffi+1, Ajj, Ald, size ), Ald, X+(XYoffi+1)*Xinc, LDX, Y+XYjj*Yinc, LDY ); } Aii = Aoffi + 1; Ajj = Aoffj + 1; XYii = XYoffi + 1; XYjj = XYoffj + 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; XYoffi += Amb; } /* * Operate with the upper triangular part of sub( A ). */ tmp1 = MIN( Aoffi, iimax ) - Aii + 1; if( upper && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); Aii += tmp1; XYii += 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; ioffxy = XYoffi; mbloc = Amb; while( ( mblkd > 0 ) && ( lcmt >= low ) ) { /* * A block owning diagonals lcmt00 >= low && lcmt00 <= upp has been found. */ if( mblkd == 1 ) mbloc = lmbloc; TRM( TYPE, SIDE, UPLO, TRANS, DIAG, mbloc, nbloc, K, lcmt, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(ioffxy+1)*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); lcmt00 = lcmt; lcmt -= pmb; mblks = mblkd; mblkd--; Aoffi = ioffd; XYoffi = ioffxy; ioffd += mbloc; ioffxy += mbloc; } /* * Operate with the lower triangular part of sub( A ). */ tmp1 = m1 - ioffd + Aii - 1; if( lower && ( tmp1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, nbloc, K, 0, ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald, X+(ioffxy+1)*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); } tmp1 = MIN( Aoffi, iimax ) - Aii + 1; m1 -= tmp1; n1 -= nbloc; lcmt00 += qnb; nblks--; Aoffj += nbloc; XYoffj += nbloc; /* * Operate with the upper triangular part of sub( A ). */ if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) ) { TRM( TYPE, SIDE, ALL, TRANS, DIAG, tmp1, n1, K, 0, ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald, X+XYii*Xinc, LDX, Y+(XYoffj+1)*Yinc, LDY ); } Aii = Aoffi + 1; Ajj = Aoffj + 1; XYii = XYoffi + 1; XYjj = XYoffj + 1; } } /* * End of PB_Cptrm */ }