ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
PB_Cpsym.c
Go to the documentation of this file.
1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS auxiliary routine (version 2.0) --
4 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
5 * and University of California, Berkeley.
6 * April 1, 1998
7 *
8 * ---------------------------------------------------------------------
9 */
10 /*
11 * Include files
12 */
13 #include "../pblas.h"
14 #include "../PBpblas.h"
15 #include "../PBtools.h"
16 #include "../PBblacs.h"
17 #include "../PBblas.h"
18 
19 #ifdef __STDC__
20 void PB_Cpsym( PBTYP_T * TYPE, PBTYP_T * UTYP, char * SIDE,
21  char * UPLO, int N, int K, char * ALPHA, char * A,
22  int IA, int JA, int * DESCA, char * XC, int LDXC,
23  char * XR, int LDXR, char * YC, int LDYC, char * YR,
24  int LDYR, TZSYM_T SYM )
25 #else
26 void PB_Cpsym( TYPE, UTYP, SIDE, UPLO, N, K, ALPHA, A, IA, JA, DESCA,
27  XC, LDXC, XR, LDXR, YC, LDYC, YR, LDYR, SYM )
28 /*
29 * .. Scalar Arguments ..
30 */
31  char * SIDE, * UPLO;
32  int IA, JA, K, LDXC, LDXR, LDYC, LDYR, N;
33  char * ALPHA;
34  PBTYP_T * TYPE, * UTYP;
35  TZSYM_T SYM;
36 /*
37 * .. Array Arguments ..
38 */
39  int * DESCA;
40  char * A, * XC, * XR, * YC, * YR;
41 #endif
42 {
43 /*
44 * Purpose
45 * =======
46 *
47 * PB_Cpsym performs a symmetric or Hermitian matrix-matrix or matrix-
48 * vector multiplication. In the following, sub( A ) denotes the symme-
49 * tric or Hermitian submatrix operand A( IA:IA+N-1, JA:JA+N-1 ).
50 *
51 * Notes
52 * =====
53 *
54 * A description vector is associated with each 2D block-cyclicly dis-
55 * tributed matrix. This vector stores the information required to
56 * establish the mapping between a matrix entry and its corresponding
57 * process and memory location.
58 *
59 * In the following comments, the character _ should be read as
60 * "of the distributed matrix". Let A be a generic term for any 2D
61 * block cyclicly distributed matrix. Its description vector is DESC_A:
62 *
63 * NOTATION STORED IN EXPLANATION
64 * ---------------- --------------- ------------------------------------
65 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
66 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
67 * the NPROW x NPCOL BLACS process grid
68 * A is distributed over. The context
69 * itself is global, but the handle
70 * (the integer value) may vary.
71 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
72 * ted matrix A, M_A >= 0.
73 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
74 * buted matrix A, N_A >= 0.
75 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
76 * block of the matrix A, IMB_A > 0.
77 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
78 * left block of the matrix A,
79 * INB_A > 0.
80 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
81 * bute the last M_A-IMB_A rows of A,
82 * MB_A > 0.
83 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
84 * bute the last N_A-INB_A columns of
85 * A, NB_A > 0.
86 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
87 * row of the matrix A is distributed,
88 * NPROW > RSRC_A >= 0.
89 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
90 * first column of A is distributed.
91 * NPCOL > CSRC_A >= 0.
92 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
93 * array storing the local blocks of
94 * the distributed matrix A,
95 * IF( Lc( 1, N_A ) > 0 )
96 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
97 * ELSE
98 * LLD_A >= 1.
99 *
100 * Let K be the number of rows of a matrix A starting at the global in-
101 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
102 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
103 * receive if these K rows were distributed over NPROW processes. If K
104 * is the number of columns of a matrix A starting at the global index
105 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
106 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
107 * these K columns were distributed over NPCOL processes.
108 *
109 * The values of Lr() and Lc() may be determined via a call to the func-
110 * tion PB_Cnumroc:
111 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
112 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
113 *
114 * Arguments
115 * =========
116 *
117 * TYPE (local input) pointer to a PBTYP_T structure
118 * On entry, TYPE is a pointer to a structure of type PBTYP_T,
119 * that contains type information (See pblas.h).
120 *
121 * UTYP (local input) pointer to a PBTYP_T structure
122 * On entry, UTYP is a pointer to a structure of type PBTYP_T,
123 * that contains type information for the Y's (See pblas.h).
124 *
125 * SIDE (global input) pointer to CHAR
126 * On entry, SIDE specifies whether op( sub( A ) ) multiplies
127 * its operand X from the left or right as follows:
128 *
129 * SIDE = 'L' or 'l' Y := alpha*op( sub( A ) )*X + Y,
130 *
131 * SIDE = 'R' or 'r' Y := alpha*X*op( sub( A ) ) + Y.
132 *
133 * UPLO (global input) pointer to CHAR
134 * On entry, UPLO specifies whether the local pieces of
135 * the array A containing the upper or lower triangular part
136 * of the symmetric or Hermitian submatrix sub( A ) are to be
137 * referenced as follows:
138 *
139 * UPLO = 'U' or 'u' Only the local pieces corresponding to
140 * the upper triangular part of the sym-
141 * metric or Hermitian submatrix sub( A )
142 * are to be referenced,
143 *
144 * UPLO = 'L' or 'l' Only the local pieces corresponding to
145 * the lower triangular part of the sym-
146 * metric or Hermitian submatrix sub( A )
147 * are to be referenced.
148 *
149 * N (global input) INTEGER
150 * On entry, N specifies the order of the submatrix sub( A ).
151 * N must be at least zero.
152 *
153 * K (global input) INTEGER
154 * On entry, K specifies the local number of columns of the lo-
155 * cal array XC and the local number of rows of the local array
156 * XR. K mut be at least zero.
157 *
158 * ALPHA (global input) pointer to CHAR
159 * On entry, ALPHA specifies the scalar alpha.
160 *
161 * A (local input) pointer to CHAR
162 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
163 * at least Lc( 1, JA+N-1 ). Before entry, this array contains
164 * the local entries of the matrix A.
165 * Before entry with UPLO = 'U' or 'u', this array contains
166 * the local entries corresponding to the upper triangular part
167 * of the @(syhec) submatrix sub( A ), and the local entries
168 * corresponding to the strictly lower triangular of sub( A )
169 * are not referenced.
170 * Before entry with UPLO = 'L' or 'l', this array contains
171 * the local entries corresponding to the lower triangular part
172 * of the @(syhec) submatrix sub( A ), and the local entries
173 * corresponding to the strictly upper triangular of sub( A )
174 * are not referenced.
175 *
176 * IA (global input) INTEGER
177 * On entry, IA specifies A's global row index, which points to
178 * the beginning of the submatrix sub( A ).
179 *
180 * JA (global input) INTEGER
181 * On entry, JA specifies A's global column index, which points
182 * to the beginning of the submatrix sub( A ).
183 *
184 * DESCA (global and local input) INTEGER array
185 * On entry, DESCA is an integer array of dimension DLEN_. This
186 * is the array descriptor for the matrix A.
187 *
188 * XC (local input) pointer to CHAR
189 * On entry, XC is an array of dimension (LDXC,K). Before entry,
190 * this array contains the local entries of the matrix XC.
191 *
192 * LDXC (local input) INTEGER
193 * On entry, LDXC specifies the leading dimension of the array
194 * XC. LDXC must be at least max( 1, Lp( IA, N ) ).
195 *
196 * XR (local input) pointer to CHAR
197 * On entry, XR is an array of dimension (LDXR,Kx), where Kx is
198 * at least Lc( JA, N ). Before entry, this array contains the
199 * local entries of the matrix XR.
200 *
201 * LDXR (local input) INTEGER
202 * On entry, LDXR specifies the leading dimension of the array
203 * XR. LDXR must be at least max( 1, K ).
204 *
205 * YC (local input/local output) pointer to CHAR
206 * On entry, YC is an array of dimension (LDYC,K). Before entry,
207 * this array contains the local entries of the matrix YC. On
208 * exit, this array contains the updated vector YC.
209 *
210 * LDYC (local input) INTEGER
211 * On entry, LDYC specifies the leading dimension of the array
212 * YC. LDYC must be at least max( 1, Lp( IA, N ) ).
213 *
214 * YR (local input/local output) pointer to CHAR
215 * On entry, YR is an array of dimension (LDYR,Ky), where Ky is
216 * at least Lc( JA, N ). Before entry, this array contains the
217 * local entries of the matrix YR. On exit, this array contains
218 * the updated vector YR.
219 *
220 * LDYR (local input) INTEGER
221 * On entry, LDYR specifies the leading dimension of the array
222 * YR. LDYR must be at least max( 1, K ).
223 *
224 * SYM (local input) pointer to function of type TZSYM_T
225 * On entry, SYM specifies the function performing the symmetric
226 * or Hermitian matrix-vector or matrix-matrix multiply of a
227 * single block.
228 *
229 * -- Written on April 1, 1998 by
230 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
231 *
232 * ---------------------------------------------------------------------
233 */
234 /*
235 * .. Local Scalars ..
236 */
237  int Acol, Arow, Aii, Aimb1, Ainb1, Ajj, Ald, Amp, Amb, Anb, Anq,
238  Aoffi, Aoffj, Arcol, Arrow, GoEast, GoSouth, IsColRepl,
239  IsRowRepl, XCinc, XRinc, Xii=0, Xjj=0, Xoffi=-1, Xoffj=-1,
240  YCinc, YRinc, iimax, ilow, imbloc, inbloc, ioffd, ioffx, iupp,
241  jjmax, joffd, joffx, lcmt, lcmt00, lmbloc, lnbloc, low, lower,
242  m1, mbloc, mblkd, mblks, mycol, myrow, n1, nbloc, nblkd,
243  nblks, npcol, nprow, pmb, qnb, size, tmp1, upp, upper;
244 /* ..
245 * .. Executable Statements ..
246 *
247 */
248 /*
249 * Quick return if possible
250 */
251  if( N <= 0 ) return;
252 /*
253 * Retrieve process grid information
254 */
255  Cblacs_gridinfo( DESCA[CTXT_], &nprow, &npcol, &myrow, &mycol );
256 /*
257 * Retrieve sub( A )'s local information: Aii, Ajj, Arow, Acol ...
258 */
259  PB_Cainfog2l( N, N, IA, JA, DESCA, nprow, npcol, myrow, mycol, &Aimb1,
260  &Ainb1, &Amp, &Anq, &Aii, &Ajj, &Arow, &Acol, &Arrow, &Arcol );
261 /*
262 * Quick return if I don't own any of sub( A ) or if sub( A ) is replicated in
263 * all processes.
264 */
265  if( ( Amp <= 0 ) || ( Anq <= 0 ) ) return;
266 
267  IsRowRepl = ( ( Arow < 0 ) || ( nprow == 1 ) );
268  IsColRepl = ( ( Acol < 0 ) || ( npcol == 1 ) );
269  Amb = DESCA[ MB_ ]; Anb = DESCA[ NB_ ]; Ald = DESCA[ LLD_ ];
270  size = TYPE->size;
271 
272  if( IsRowRepl && IsColRepl )
273  {
274  SYM( TYPE, SIDE, UPLO, Amp, Anq, K, 0, ALPHA, Mptr( A, Aii, Ajj, Ald,
275  size ), Ald, XC, LDXC, XR, LDXR, YC, LDYC, YR, LDYR );
276  return;
277  }
278 
279  XCinc = size; XRinc = LDXR * size;
280  YCinc = UTYP->size; YRinc = LDYR * UTYP->size;
281  upper = ( Mupcase( UPLO[0] ) == CUPPER );
282  lower = ( Mupcase( UPLO[0] ) == CLOWER );
283 /*
284 * Initialize lcmt00, mblks, nblks, imbloc, inbloc, lmbloc, lnbloc, ilow, low,
285 * iupp, and upp.
286 */
287  PB_Cbinfo( 0, Amp, Anq, Aimb1, Ainb1, Amb, Anb, Arrow, Arcol, &lcmt00,
288  &mblks, &nblks, &imbloc, &inbloc, &lmbloc, &lnbloc, &ilow, &low,
289  &iupp, &upp );
290 
291  iimax = ( Aoffi = Aii - 1 ) + ( m1 = Amp );
292  jjmax = ( Aoffj = Ajj - 1 ) + ( n1 = Anq );
293  pmb = ( IsRowRepl ? Amb : nprow * Amb );
294  qnb = ( IsColRepl ? Anb : npcol * Anb );
295 /*
296 * Handle separately the first row and/or column of the LCM table. Update the
297 * LCM value of the curent block lcmt00, as well as the number of rows and
298 * columns mblks and nblks remaining in the LCM table.
299 */
300  GoSouth = ( lcmt00 > iupp );
301  GoEast = ( lcmt00 < ilow );
302 /*
303 * Go through the table looking for blocks owning diagonal entries.
304 */
305  if( ( !( GoSouth ) ) && ( !( GoEast ) ) )
306  {
307 /*
308 * The upper left block owns diagonal entries lcmt00 >= ilow && lcmt00 <= iupp
309 */
310  SYM( TYPE, SIDE, UPLO, imbloc, inbloc, K, lcmt00, ALPHA, Mptr( A, Aii,
311  Ajj, Ald, size ), Ald, XC+Xii*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
312  YC+Xii*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
313 /*
314 * Decide whether one should go south or east in the table: Go east if the
315 * block below the current one only owns lower entries. If this block, however,
316 * owns diagonals, then go south.
317 */
318  GoSouth = !( GoEast = ( ( lcmt00 - ( iupp - upp + pmb ) ) < ilow ) );
319 
320  if( GoSouth )
321  {
322 /*
323 * When the upper triangular part of sub( A ) should be operated with and
324 * one is planning to go south in the table, it is neccessary to take care
325 * of the remaining columns of these imbloc rows immediately.
326 */
327  if( upper && ( Anq > inbloc ) )
328  {
329  tmp1 = Anq - inbloc;
330  SYM( TYPE, SIDE, ALL, imbloc, tmp1, K, 0,
331  ALPHA, Mptr( A, Aii, Ajj+inbloc, Ald, size ), Ald,
332  XC+Xii*XCinc, LDXC, XR+(Xjj+inbloc)*XRinc, LDXR,
333  YC+Xii*YCinc, LDYC, YR+(Xjj+inbloc)*YRinc, LDYR );
334  }
335  Aii += imbloc; Xii += imbloc; m1 -= imbloc;
336  }
337  else
338  {
339 /*
340 * When the lower triangular part of sub( A ) should be operated with and
341 * one is planning to go east in the table, it is neccessary to take care
342 * of the remaining rows of these inbloc columns immediately.
343 */
344  if( lower && ( Amp > imbloc ) )
345  {
346  tmp1 = Amp - imbloc;
347  SYM( TYPE, SIDE, ALL, tmp1, inbloc, K, 0,
348  ALPHA, Mptr( A, Aii+imbloc, Ajj, Ald, size ), Ald,
349  XC+(Xii+imbloc)*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
350  YC+(Xii+imbloc)*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
351  }
352  Ajj += inbloc; Xjj += inbloc; n1 -= inbloc;
353  }
354  }
355 
356  if( GoSouth )
357  {
358 /*
359 * Go one step south in the LCM table. Adjust the current LCM value as well as
360 * the local row indexes in A and XC.
361 */
362  lcmt00 -= ( iupp - upp + pmb ); mblks--;
363  Aoffi += imbloc; Xoffi += imbloc;
364 /*
365 * While there are blocks remaining that own upper entries, keep going south.
366 * Adjust the current LCM value as well as the local row indexes in A and XC.
367 */
368  while( ( mblks > 0 ) && ( lcmt00 > upp ) )
369  { lcmt00 -= pmb; mblks--; Aoffi += Amb; Xoffi += Amb; }
370 /*
371 * Operate with the upper triangular part of sub( A ) we just skipped when
372 * necessary.
373 */
374  tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
375  if( upper && ( tmp1 > 0 ) )
376  {
377  SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
378  ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
379  XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
380  YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
381  Aii += tmp1; Xii += tmp1; m1 -= tmp1;
382  }
383 /*
384 * Return if no more row in the LCM table.
385 */
386  if( mblks <= 0 ) return;
387 /*
388 * lcmt00 <= upp. The current block owns either diagonals or lower entries.
389 * Save the current position in the LCM table. After this column has been
390 * completely taken care of, re-start from this row and the next column of
391 * the LCM table.
392 */
393  lcmt = lcmt00; mblkd = mblks; ioffd = Aoffi; ioffx = Xoffi;
394 
395  mbloc = Amb;
396  while( ( mblkd > 0 ) && ( lcmt >= ilow ) )
397  {
398 /*
399 * A block owning diagonals lcmt00 >= ilow && lcmt00 <= upp has been found.
400 */
401  if( mblkd == 1 ) mbloc = lmbloc;
402  SYM( TYPE, SIDE, UPLO, mbloc, inbloc, K, lcmt,
403  ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
404  XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
405  YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
406  lcmt00 = lcmt; lcmt -= pmb;
407  mblks = mblkd; mblkd--;
408  Aoffi = ioffd; ioffd += mbloc;
409  Xoffi = ioffx; ioffx += mbloc;
410  }
411 /*
412 * Operate with the lower triangular part of sub( A ).
413 */
414  tmp1 = m1 - ioffd + Aii - 1;
415  if( lower && ( tmp1 > 0 ) )
416  SYM( TYPE, SIDE, ALL, tmp1, inbloc, K, 0,
417  ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
418  XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
419  YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
420 
421  tmp1 = Aoffi - Aii + 1;
422  m1 -= tmp1;
423  n1 -= inbloc;
424  lcmt00 += low - ilow + qnb;
425  nblks--;
426  Aoffj += inbloc;
427  Xoffj += inbloc;
428 /*
429 * Operate with the upper triangular part of sub( A ).
430 */
431  if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
432  SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
433  ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
434  XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
435  YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
436  Aii = Aoffi + 1; Ajj = Aoffj + 1;
437  Xii = Xoffi + 1; Xjj = Xoffj + 1;
438  }
439  else if( GoEast )
440  {
441 /*
442 * Go one step east in the LCM table. Adjust the current LCM value as well as
443 * the local column index in A and XR.
444 */
445  lcmt00 += low - ilow + qnb; nblks--;
446  Aoffj += inbloc; Xoffj += inbloc;
447 /*
448 * While there are blocks remaining that own lower entries, keep going east.
449 * Adjust the current LCM value as well as the local column index in A and XR.
450 */
451  while( ( nblks > 0 ) && ( lcmt00 < low ) )
452  { lcmt00 += qnb; nblks--; Aoffj += Anb; Xoffj += Anb; }
453 /*
454 * Operate with the lower triangular part of sub( A ).
455 */
456  tmp1 = MIN( Aoffj, jjmax ) - Ajj + 1;
457  if( lower && ( tmp1 > 0 ) )
458  {
459  SYM( TYPE, SIDE, ALL, m1, tmp1, K, 0,
460  ALPHA, Mptr( A, Aii, Ajj, Ald, size ), Ald,
461  XC+Xii*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
462  YC+Xii*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
463  Ajj += tmp1; Xjj += tmp1; n1 -= tmp1;
464  }
465 /*
466 * Return if no more column in the LCM table.
467 */
468  if( nblks <= 0 ) return;
469 /*
470 * lcmt00 >= low. The current block owns either diagonals or upper entries.
471 * Save the current position in the LCM table. After this row has been
472 * completely taken care of, re-start from this column and the next row of
473 * the LCM table.
474 */
475  lcmt = lcmt00; nblkd = nblks; joffd = Aoffj; joffx = Xoffj;
476 
477  nbloc = Anb;
478  while( ( nblkd > 0 ) && ( lcmt <= iupp ) )
479  {
480 /*
481 * A block owning diagonals lcmt00 >= low && lcmt00 <= iupp has been found.
482 */
483  if( nblkd == 1 ) nbloc = lnbloc;
484  SYM( TYPE, SIDE, UPLO, imbloc, nbloc, K, lcmt,
485  ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald,
486  XC+Xii*XCinc, LDXC, XR+(joffx+1)*XRinc, LDXR,
487  YC+Xii*YCinc, LDYC, YR+(joffx+1)*YRinc, LDYR );
488  lcmt00 = lcmt; lcmt += qnb;
489  nblks = nblkd; nblkd--;
490  Aoffj = joffd; joffd += nbloc;
491  Xoffj = joffx; joffx += nbloc;
492  }
493 /*
494 * Operate with the upper triangular part of sub( A ).
495 */
496  tmp1 = n1 - joffd + Ajj - 1;
497  if( upper && ( tmp1 > 0 ) )
498  SYM( TYPE, SIDE, ALL, imbloc, tmp1, K, 0,
499  ALPHA, Mptr( A, Aii, joffd+1, Ald, size ), Ald,
500  XC+Xii*XCinc, LDXC, XR+(joffx+1)*XRinc, LDXR,
501  YC+Xii*YCinc, LDYC, YR+(joffx+1)*YRinc, LDYR );
502 
503  tmp1 = Aoffj - Ajj + 1;
504  m1 -= imbloc;
505  n1 -= tmp1;
506  lcmt00 -= ( iupp - upp + pmb );
507  mblks--;
508  Aoffi += imbloc;
509  Xoffi += imbloc;
510 /*
511 * Operate with the lower triangular part of sub( A ).
512 */
513  if( lower && ( m1 > 0 ) && ( tmp1 > 0 ) )
514  SYM( TYPE, SIDE, ALL, m1, tmp1, K, 0,
515  ALPHA, Mptr( A, Aoffi+1, Ajj, Ald, size ), Ald,
516  XC+(Xoffi+1)*XCinc, LDXC, XR+Xjj*XRinc, LDXR,
517  YC+(Xoffi+1)*YCinc, LDYC, YR+Xjj*YRinc, LDYR );
518  Aii = Aoffi + 1; Ajj = Aoffj + 1;
519  Xii = Xoffi + 1; Xjj = Xoffj + 1;
520  }
521 /*
522 * Loop over the remaining columns of the LCM table.
523 */
524  nbloc = Anb;
525  while( nblks > 0 )
526  {
527  if( nblks == 1 ) nbloc = lnbloc;
528 /*
529 * While there are blocks remaining that own upper entries, keep going south.
530 * Adjust the current LCM value as well as the local row index in A and XC.
531 */
532  while( ( mblks > 0 ) && ( lcmt00 > upp ) )
533  { lcmt00 -= pmb; mblks--; Aoffi += Amb; Xoffi += Amb; }
534 /*
535 * Operate with the upper triangular part of sub( A ).
536 */
537  tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
538  if( upper && ( tmp1 > 0 ) )
539  {
540  SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
541  ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
542  XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
543  YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
544  Aii += tmp1;
545  Xii += tmp1;
546  m1 -= tmp1;
547  }
548 /*
549 * Return if no more row in the LCM table.
550 */
551  if( mblks <= 0 ) return;
552 /*
553 * lcmt00 <= upp. The current block owns either diagonals or lower entries.
554 * Save the current position in the LCM table. After this column has been
555 * completely taken care of, re-start from this row and the next column of
556 * the LCM table.
557 */
558  lcmt = lcmt00; mblkd = mblks; ioffd = Aoffi; ioffx = Xoffi;
559 
560  mbloc = Amb;
561  while( ( mblkd > 0 ) && ( lcmt >= low ) )
562  {
563 /*
564 * A block owning diagonals lcmt00 >= low && lcmt00 <= upp has been found.
565 */
566  if( mblkd == 1 ) mbloc = lmbloc;
567  SYM( TYPE, SIDE, UPLO, mbloc, nbloc, K, lcmt,
568  ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
569  XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
570  YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
571  lcmt00 = lcmt; lcmt -= pmb;
572  mblks = mblkd; mblkd--;
573  Aoffi = ioffd; Xoffi = ioffx;
574  ioffd += mbloc; ioffx += mbloc;
575  }
576 /*
577 * Operate with the lower triangular part of sub( A ).
578 */
579  tmp1 = m1 - ioffd + Aii - 1;
580  if( lower && ( tmp1 > 0 ) )
581  SYM( TYPE, SIDE, ALL, tmp1, nbloc, K, 0,
582  ALPHA, Mptr( A, ioffd+1, Aoffj+1, Ald, size ), Ald,
583  XC+(ioffx+1)*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
584  YC+(ioffx+1)*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
585 
586  tmp1 = MIN( Aoffi, iimax ) - Aii + 1;
587  m1 -= tmp1;
588  n1 -= nbloc;
589  lcmt00 += qnb;
590  nblks--;
591  Aoffj += nbloc;
592  Xoffj += nbloc;
593 /*
594 * Operate with the upper triangular part of sub( A ).
595 */
596  if( upper && ( tmp1 > 0 ) && ( n1 > 0 ) )
597  SYM( TYPE, SIDE, ALL, tmp1, n1, K, 0,
598  ALPHA, Mptr( A, Aii, Aoffj+1, Ald, size ), Ald,
599  XC+Xii*XCinc, LDXC, XR+(Xoffj+1)*XRinc, LDXR,
600  YC+Xii*YCinc, LDYC, YR+(Xoffj+1)*YRinc, LDYR );
601  Aii = Aoffi + 1; Ajj = Aoffj + 1;
602  Xii = Xoffi + 1; Xjj = Xoffj + 1;
603  }
604 /*
605 * End of PB_Cpsym
606 */
607 }
TYPE
#define TYPE
Definition: clamov.c:7
MB_
#define MB_
Definition: PBtools.h:43
NB_
#define NB_
Definition: PBtools.h:44
PB_Cainfog2l
void PB_Cainfog2l()
TZSYM_T
void(* TZSYM_T)()
Definition: pblas.h:428
LLD_
#define LLD_
Definition: PBtools.h:47
PB_Cpsym
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: PB_Cpsym.c:26
PB_Cbinfo
void PB_Cbinfo()
CLOWER
#define CLOWER
Definition: PBblas.h:25
PBTYP_T::size
int size
Definition: pblas.h:329
ALL
#define ALL
Definition: PBblas.h:50
MIN
#define MIN(a_, b_)
Definition: PBtools.h:76
Cblacs_gridinfo
void Cblacs_gridinfo()
PBTYP_T
Definition: pblas.h:325
Mupcase
#define Mupcase(C)
Definition: PBtools.h:83
CUPPER
#define CUPPER
Definition: PBblas.h:26
Mptr
#define Mptr(a_, i_, j_, lda_, siz_)
Definition: PBtools.h:132
CTXT_
#define CTXT_
Definition: PBtools.h:38