ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pcahemv_.c
Go to the documentation of this file.
1 /* ---------------------------------------------------------------------
2 *
3 * -- PBLAS 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 pcahemv_( F_CHAR_T UPLO, int * N, float * ALPHA,
21  float * A, int * IA, int * JA, int * DESCA,
22  float * X, int * IX, int * JX, int * DESCX, int * INCX,
23  float * BETA,
24  float * Y, int * IY, int * JY, int * DESCY, int * INCY )
25 #else
26 void pcahemv_( UPLO, N, ALPHA, A, IA, JA, DESCA, X, IX, JX, DESCX,
27  INCX, BETA, Y, IY, JY, DESCY, INCY )
28 /*
29 * .. Scalar Arguments ..
30 */
31  F_CHAR_T UPLO;
32  int * IA, * INCX, * INCY, * IX, * IY, * JA, * JX, * JY,
33  * N;
34  float * ALPHA, * BETA;
35 /*
36 * .. Array Arguments ..
37 */
38  int * DESCA, * DESCX, * DESCY;
39  float * A, * X, * Y;
40 #endif
41 {
42 /*
43 * Purpose
44 * =======
45 *
46 * PCAHEMV performs the matrix-vector operation
47 *
48 * sub( Y ) := abs( alpha )*abs( sub( A ) )*abs( sub( X ) ) +
49 * abs( beta*sub( Y ) ),
50 *
51 * where
52 *
53 * sub( A ) denotes A(IA:IA+M-1,JA:JA+N-1),
54 *
55 * sub( X ) denotes X(IX,JX:JX+N-1) if INCX = M_X,
56 * X(IX:IX+N-1,JX) if INCX = 1 and INCX <> M_X, and,
57 *
58 * sub( Y ) denotes Y(IY,JY:JY+N-1) if INCY = M_Y,
59 * Y(IY:IY+N-1,JY) if INCY = 1 and INCY <> M_Y.
60 *
61 * Alpha and beta are real scalars, sub( Y ) is a n element real subvec-
62 * tor, sub( X ) is an n element subvector and sub( A ) is an n by n
63 * Hermitian submatrix.
64 *
65 * Notes
66 * =====
67 *
68 * A description vector is associated with each 2D block-cyclicly dis-
69 * tributed matrix. This vector stores the information required to
70 * establish the mapping between a matrix entry and its corresponding
71 * process and memory location.
72 *
73 * In the following comments, the character _ should be read as
74 * "of the distributed matrix". Let A be a generic term for any 2D
75 * block cyclicly distributed matrix. Its description vector is DESC_A:
76 *
77 * NOTATION STORED IN EXPLANATION
78 * ---------------- --------------- ------------------------------------
79 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
80 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
81 * the NPROW x NPCOL BLACS process grid
82 * A is distributed over. The context
83 * itself is global, but the handle
84 * (the integer value) may vary.
85 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
86 * ted matrix A, M_A >= 0.
87 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
88 * buted matrix A, N_A >= 0.
89 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
90 * block of the matrix A, IMB_A > 0.
91 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
92 * left block of the matrix A,
93 * INB_A > 0.
94 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
95 * bute the last M_A-IMB_A rows of A,
96 * MB_A > 0.
97 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
98 * bute the last N_A-INB_A columns of
99 * A, NB_A > 0.
100 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
101 * row of the matrix A is distributed,
102 * NPROW > RSRC_A >= 0.
103 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
104 * first column of A is distributed.
105 * NPCOL > CSRC_A >= 0.
106 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
107 * array storing the local blocks of
108 * the distributed matrix A,
109 * IF( Lc( 1, N_A ) > 0 )
110 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
111 * ELSE
112 * LLD_A >= 1.
113 *
114 * Let K be the number of rows of a matrix A starting at the global in-
115 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
116 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
117 * receive if these K rows were distributed over NPROW processes. If K
118 * is the number of columns of a matrix A starting at the global index
119 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
120 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
121 * these K columns were distributed over NPCOL processes.
122 *
123 * The values of Lr() and Lc() may be determined via a call to the func-
124 * tion PB_Cnumroc:
125 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
126 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
127 *
128 * Arguments
129 * =========
130 *
131 * UPLO (global input) CHARACTER*1
132 * On entry, UPLO specifies whether the local pieces of
133 * the array A containing the upper or lower triangular part
134 * of the Hermitian submatrix sub( A ) are to be referenced as
135 * follows:
136 *
137 * UPLO = 'U' or 'u' Only the local pieces corresponding to
138 * the upper triangular part of the
139 * Hermitian submatrix sub( A ) are to be
140 * referenced,
141 *
142 * UPLO = 'L' or 'l' Only the local pieces corresponding to
143 * the lower triangular part of the
144 * Hermitian submatrix sub( A ) are to be
145 * referenced.
146 *
147 * N (global input) INTEGER
148 * On entry, N specifies the order of the submatrix sub( A ).
149 * N must be at least zero.
150 *
151 * ALPHA (global input) REAL
152 * On entry, ALPHA specifies the scalar alpha. When ALPHA is
153 * supplied as zero then the local entries of the arrays A
154 * and X corresponding to the entries of the submatrix sub( A )
155 * and the subvector sub( X ) need not be set on input.
156 *
157 * A (local input) COMPLEX array
158 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
159 * at least Lc( 1, JA+N-1 ). Before entry, this array contains
160 * the local entries of the matrix A.
161 * Before entry with UPLO = 'U' or 'u', this array contains
162 * the local entries of the upper triangular part of the
163 * Hermitian submatrix sub( A ), and the local entries of the
164 * strictly lower triangular of sub( A ) are not referenced.
165 * Before entry with UPLO = 'L' or 'l', this array contains
166 * the local entries of the lower triangular part of the
167 * Hermitian submatrix sub( A ), and the local entries of the
168 * strictly upper triangular of sub( A ) are not referenced.
169 * Note that the imaginary parts of the local entries corres-
170 * ponding to the diagonal elements of sub( A ) need not be
171 * set and assumed to be zero.
172 *
173 * IA (global input) INTEGER
174 * On entry, IA specifies A's global row index, which points to
175 * the beginning of the submatrix sub( A ).
176 *
177 * JA (global input) INTEGER
178 * On entry, JA specifies A's global column index, which points
179 * to the beginning of the submatrix sub( A ).
180 *
181 * DESCA (global and local input) INTEGER array
182 * On entry, DESCA is an integer array of dimension DLEN_. This
183 * is the array descriptor for the matrix A.
184 *
185 * X (local input) COMPLEX array
186 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
187 * is at least MAX( 1, Lr( 1, IX ) ) when INCX = M_X and
188 * MAX( 1, Lr( 1, IX+N-1 ) ) otherwise, and, Kx is at least
189 * Lc( 1, JX+N-1 ) when INCX = M_X and Lc( 1, JX ) otherwise.
190 * Before entry, this array contains the local entries of the
191 * matrix X.
192 *
193 * IX (global input) INTEGER
194 * On entry, IX specifies X's global row index, which points to
195 * the beginning of the submatrix sub( X ).
196 *
197 * JX (global input) INTEGER
198 * On entry, JX specifies X's global column index, which points
199 * to the beginning of the submatrix sub( X ).
200 *
201 * DESCX (global and local input) INTEGER array
202 * On entry, DESCX is an integer array of dimension DLEN_. This
203 * is the array descriptor for the matrix X.
204 *
205 * INCX (global input) INTEGER
206 * On entry, INCX specifies the global increment for the
207 * elements of X. Only two values of INCX are supported in
208 * this version, namely 1 and M_X. INCX must not be zero.
209 *
210 * BETA (global input) REAL
211 * On entry, BETA specifies the scalar beta. When BETA is
212 * supplied as zero then the local entries of the array Y
213 * corresponding to the entries of the subvector sub( Y ) need
214 * not be set on input.
215 *
216 * Y (local input/local output) REAL array
217 * On entry, Y is an array of dimension (LLD_Y, Ky), where LLD_Y
218 * is at least MAX( 1, Lr( 1, IY ) ) when INCY = M_Y and
219 * MAX( 1, Lr( 1, IY+N-1 ) ) otherwise, and, Ky is at least
220 * Lc( 1, JY+N-1 ) when INCY = M_Y and Lc( 1, JY ) otherwise.
221 * Before entry, this array contains the local entries of the
222 * matrix Y. On exit, sub( Y ) is overwritten by the updated
223 * subvector.
224 *
225 * IY (global input) INTEGER
226 * On entry, IY specifies Y's global row index, which points to
227 * the beginning of the submatrix sub( Y ).
228 *
229 * JY (global input) INTEGER
230 * On entry, JY specifies Y's global column index, which points
231 * to the beginning of the submatrix sub( Y ).
232 *
233 * DESCY (global and local input) INTEGER array
234 * On entry, DESCY is an integer array of dimension DLEN_. This
235 * is the array descriptor for the matrix Y.
236 *
237 * INCY (global input) INTEGER
238 * On entry, INCY specifies the global increment for the
239 * elements of Y. Only two values of INCY are supported in
240 * this version, namely 1 and M_Y. INCY must not be zero.
241 *
242 * -- Written on April 1, 1998 by
243 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
244 *
245 * ---------------------------------------------------------------------
246 */
247 /*
248 * .. Local Scalars ..
249 */
250  char UploA, * one, top;
251  int Acol, Ai, Aii, Aimb1, Ainb1, Aj, Ajj, Akp, Akq, Ald, Amb, Amp,
252  Amp0, Anb, Anq, Anq0, Arow, XCfr, XCld, XRfr, XRld, Xi, Xj,
253  YCfr, YCpbY, YCsum, YCld, YRfr, YRpbY, YRsum, YRld, Ycol, Yi,
254  Yii, YisRow, Yj, Yjj, Yld, Ynp, Ynq, Yrow, ctxt, info, ione=1,
255  k, kb, ktmp, mycol, myrow, nb, npcol, nprow, size, upper,
256  usiz;
257  float * tbeta;
258  PBTYP_T * type, * utyp;
259 /*
260 * .. Local Arrays ..
261 */
262  int Ad [DLEN_], Ad0[DLEN_], XCd[DLEN_], XRd[DLEN_], Xd[DLEN_],
263  YCd[DLEN_], YRd[DLEN_], Yd [DLEN_];
264  char * Aptr = NULL, * XC = NULL, * XR = NULL, * YC = NULL,
265  * YR = NULL;
266 /* ..
267 * .. Executable Statements ..
268 *
269 */
270  upper = ( ( UploA = Mupcase( F2C_CHAR( UPLO )[0] ) ) == CUPPER );
271  PB_CargFtoC( *IA, *JA, DESCA, &Ai, &Aj, Ad );
272  PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
273  PB_CargFtoC( *IY, *JY, DESCY, &Yi, &Yj, Yd );
274 #ifndef NO_ARGCHK
275 /*
276 * Test the input parameters
277 */
278  Cblacs_gridinfo( ( ctxt = Ad[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
279  if( !( info = ( ( nprow == -1 ) ? -( 701 + CTXT_ ) : 0 ) ) )
280  {
281  if( ( !upper ) && ( UploA != CLOWER ) )
282  {
283  PB_Cwarn( ctxt, __LINE__, __FILE__, "Illegal UPLO = %c\n", UploA );
284  info = -1;
285  }
286  PB_Cchkmat( ctxt, "PCAHEMV", "A", *N, 2, *N, 2, Ai, Aj, Ad, 7, &info );
287  PB_Cchkvec( ctxt, "PCAHEMV", "X", *N, 2, Xi, Xj, Xd, *INCX, 11, &info );
288  PB_Cchkvec( ctxt, "PCAHEMV", "Y", *N, 2, Yi, Yj, Yd, *INCY, 17, &info );
289  }
290  if( info ) { PB_Cabort( ctxt, "PCAHEMV", info ); return; }
291 #endif
292 /*
293 * Quick return if possible
294 */
295  if( ( *N == 0 ) ||
296  ( ( ALPHA[REAL_PART] == ZERO ) && ( BETA[REAL_PART] == ONE ) ) )
297  return;
298 /*
299 * Retrieve process grid information
300 */
301 #ifdef NO_ARGCHK
302  Cblacs_gridinfo( ( ctxt = Ad[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
303 #endif
304 /*
305 * Get type structure
306 */
307  type = PB_Cctypeset(); size = type->size;
308  utyp = PB_Cstypeset(); usiz = type->usiz;
309 /*
310 * and when alpha is zero
311 */
312  if( ALPHA[REAL_PART] == ZERO )
313  {
314 /*
315 * Retrieve sub( Y )'s local information: Yii, Yjj, Yrow, Ycol
316 */
317  PB_Cinfog2l( Yi, Yj, Yd, nprow, npcol, myrow, mycol, &Yii, &Yjj,
318  &Yrow, &Ycol );
319 
320  if( *INCY == Yd[M_] )
321  {
322 /*
323 * sub( Y ) resides in (a) process row(s)
324 */
325  if( ( myrow == Yrow ) || ( Yrow < 0 ) )
326  {
327 /*
328 * Make sure I own some data and scale sub( Y )
329 */
330  Ynq = PB_Cnumroc( *N, Yj, Yd[INB_], Yd[NB_], mycol, Yd[CSRC_],
331  npcol );
332  if( Ynq > 0 )
333  {
334  Yld = Yd[LLD_];
335  sascal_( &Ynq, ((char *) BETA), Mptr( ((char *) Y), Yii,
336  Yjj, Yld, usiz ), &Yld );
337  }
338  }
339  }
340  else
341  {
342 /*
343 * sub( Y ) resides in (a) process column(s)
344 */
345  if( ( mycol == Ycol ) || ( Ycol < 0 ) )
346  {
347 /*
348 * Make sure I own some data and scale sub( Y )
349 */
350  Ynp = PB_Cnumroc( *N, Yi, Yd[IMB_], Yd[MB_], myrow, Yd[RSRC_],
351  nprow );
352  if( Ynp > 0 )
353  {
354  sascal_( &Ynp, ((char *) BETA), Mptr( ((char *) Y), Yii,
355  Yjj, Yd[LLD_], usiz ), INCY );
356  }
357  }
358  }
359  return;
360  }
361 /*
362 * Compute descriptor Ad0 for sub( A )
363 */
364  PB_Cdescribe( *N, *N, Ai, Aj, Ad, nprow, npcol, myrow, mycol, &Aii, &Ajj,
365  &Ald, &Aimb1, &Ainb1, &Amb, &Anb, &Arow, &Acol, Ad0 );
366 /*
367 * Reuse sub( Y ) and/or create vectors YR in process rows and YC in process
368 * columns spanned by sub( A )
369 */
370  if( ( YisRow = ( *INCY == Yd[M_] ) ) != 0 )
371  {
372  PB_CInOutV( utyp, ROW, *N, *N, Ad0, 1, ((char *)BETA), ((char *) Y),
373  Yi, Yj, Yd, ROW, ((char**)(&tbeta)), &YR, YRd, &YRfr,
374  &YRsum, &YRpbY );
375  PB_COutV( utyp, COLUMN, INIT, *N, *N, Ad0, 1, &YC, YCd, &YCfr, &YCsum );
376  }
377  else
378  {
379  PB_CInOutV( utyp, COLUMN, *N, *N, Ad0, 1, ((char *)BETA), ((char *) Y),
380  Yi, Yj, Yd, COLUMN, ((char**)(&tbeta)), &YC, YCd, &YCfr,
381  &YCsum, &YCpbY );
382  PB_COutV( utyp, ROW, INIT, *N, *N, Ad0, 1, &YR, YRd, &YRfr, &YRsum );
383  }
384 /*
385 * Replicate sub( X ) in process rows (XR) and process columns (XC) spanned by
386 * sub( A )
387 */
388  if( *INCX == Xd[M_] )
389  {
390  PB_CInV( type, NOCONJG, ROW, *N, *N, Ad0, 1, ((char *) X), Xi, Xj, Xd,
391  ROW, &XR, XRd, &XRfr );
392  PB_CInV( type, NOCONJG, COLUMN, *N, *N, Ad0, 1, XR, 0, 0, XRd,
393  ROW, &XC, XCd, &XCfr );
394  }
395  else
396  {
397  PB_CInV( type, NOCONJG, COLUMN, *N, *N, Ad0, 1, ((char *) X), Xi, Xj, Xd,
398  COLUMN, &XC, XCd, &XCfr );
399  PB_CInV( type, NOCONJG, ROW, *N, *N, Ad0, 1, XC, 0, 0, XCd,
400  COLUMN, &XR, XRd, &XRfr );
401  }
402  one = type->one;
403 /*
404 * Local matrix-vector multiply iff I own some data
405 */
406  Aimb1 = Ad0[IMB_ ]; Ainb1 = Ad0[INB_ ]; Amb = Ad0[MB_]; Anb = Ad0[NB_];
407  Acol = Ad0[CSRC_]; Arow = Ad0[RSRC_];
408  Amp = PB_Cnumroc( *N, 0, Aimb1, Amb, myrow, Arow, nprow );
409  Anq = PB_Cnumroc( *N, 0, Ainb1, Anb, mycol, Acol, npcol );
410 
411  if( ( Amp > 0 ) && ( Anq > 0 ) )
412  {
413  Aptr = Mptr( ((char *) A), Aii, Ajj, Ald, size );
414 
415  XCld = XCd[LLD_]; XRld = XRd[LLD_]; YCld = YCd[LLD_]; YRld = YRd[LLD_];
416 /*
417 * Scale YR or YC in the case sub( Y ) has been reused
418 */
419  if( YisRow )
420  {
421 /*
422 * YR resides in (a) process row(s)
423 */
424  if( !YRpbY )
425  {
426  if( ( myrow == YRd[RSRC_] ) || ( YRd[RSRC_] < 0 ) )
427  {
428 /*
429 * Make sure I own some data and scale YR
430 */
431  if( Anq > 0 )
432  sascal_( &Anq, ((char *) tbeta), YR, &YRld );
433  }
434  }
435  }
436  else
437  {
438 /*
439 * YC resides in (a) process column(s)
440 */
441  if( !YCpbY )
442  {
443  if( ( mycol == YCd[CSRC_] ) || ( YCd[CSRC_] < 0 ) )
444  {
445 /*
446 * Make sure I own some data and scale YC
447 */
448  if( Amp > 0 )
449  sascal_( &Amp, ((char *) tbeta), YC, &ione );
450  }
451  }
452  }
453 /*
454 * Computational partitioning size is computed as the product of the logical
455 * value returned by pilaenv_ and 2 * lcm( nprow, npcol )
456 */
457  nb = 2 * pilaenv_( &ctxt, C2F_CHAR( &utyp->type ) ) *
458  PB_Clcm( ( Arow >= 0 ? nprow : 1 ), ( Acol >= 0 ? npcol : 1 ) );
459 
460  if( upper )
461  {
462  for( k = 0; k < *N; k += nb )
463  {
464  kb = *N - k; kb = MIN( kb, nb );
465  Akp = PB_Cnumroc( k, 0, Aimb1, Amb, myrow, Arow, nprow );
466  Akq = PB_Cnumroc( k, 0, Ainb1, Anb, mycol, Acol, npcol );
467  Anq0 = PB_Cnumroc( kb, k, Ainb1, Anb, mycol, Acol, npcol );
468  if( Akp > 0 && Anq0 > 0 )
469  {
470  cagemv_( C2F_CHAR( NOTRAN ), &Akp, &Anq0, ((char *)ALPHA),
471  Mptr( Aptr, 0, Akq, Ald, size ), &Ald, Mptr( XR, 0, Akq,
472  XRld, size ), &XRld, one, YC, &ione );
473  cagemv_( C2F_CHAR( COTRAN ), &Akp, &Anq0, ((char *)ALPHA),
474  Mptr( Aptr, 0, Akq, Ald, size ), &Ald, XC, &ione, one,
475  Mptr( YR, 0, Akq, YRld, usiz ), &YRld );
476  }
477  PB_Cpsym( type, utyp, LEFT, UPPER, kb, 1, ((char *) ALPHA), Aptr, k,
478  k, Ad0, Mptr( XC, Akp, 0, XCld, size ), XCld, Mptr( XR, 0,
479  Akq, XRld, size ), XRld, Mptr( YC, Akp, 0, YCld, usiz ),
480  YCld, Mptr( YR, 0, Akq, YRld, usiz ), YRld, PB_Ctzahemv );
481  }
482  }
483  else
484  {
485  for( k = 0; k < *N; k += nb )
486  {
487  kb = *N - k; ktmp = k + ( kb = MIN( kb, nb ) );
488  Akp = PB_Cnumroc( k, 0, Aimb1, Amb, myrow, Arow, nprow );
489  Akq = PB_Cnumroc( k, 0, Ainb1, Anb, mycol, Acol, npcol );
490  PB_Cpsym( type, utyp, LEFT, LOWER, kb, 1, ((char *) ALPHA), Aptr, k,
491  k, Ad0, Mptr( XC, Akp, 0, XCld, size ), XCld, Mptr( XR, 0,
492  Akq, XRld, size ), XRld, Mptr( YC, Akp, 0, YCld, usiz ),
493  YCld, Mptr( YR, 0, Akq, YRld, usiz ), YRld, PB_Ctzahemv );
494  Akp = PB_Cnumroc( ktmp, 0, Aimb1, Amb, myrow, Arow, nprow );
495  Amp0 = Amp - Akp;
496  Anq0 = PB_Cnumroc( kb, k, Ainb1, Anb, mycol, Acol, npcol );
497  if( Amp0 > 0 && Anq0 > 0 )
498  {
499  cagemv_( C2F_CHAR( NOTRAN ), &Amp0, &Anq0, ((char *) ALPHA),
500  Mptr( Aptr, Akp, Akq, Ald, size ), &Ald, Mptr( XR, 0,
501  Akq, XRld, size ), &XRld, one, Mptr( YC, Akp, 0, YCld,
502  usiz ), &ione );
503  cagemv_( C2F_CHAR( COTRAN ), &Amp0, &Anq0, ((char *) ALPHA),
504  Mptr( Aptr, Akp, Akq, Ald, size ), &Ald, Mptr( XC, Akp,
505  0, XCld, size ), &ione, one, Mptr( YR, 0, Akq, YRld,
506  usiz ), &YRld );
507  }
508  }
509  }
510  }
511  if( XCfr ) free( XC );
512  if( XRfr ) free( XR );
513 
514  if( YisRow )
515  {
516 /*
517 * Combine the partial column results into YC
518 */
519  if( YCsum )
520  {
521  YCd[CSRC_] = 0;
522  if( Amp > 0 )
523  {
524  top = *PB_Ctop( &ctxt, COMBINE, ROW, TOP_GET );
525  Csgsum2d( ctxt, ROW, &top, Amp, 1, YC, YCd[LLD_], myrow, 0 );
526  }
527  }
528 /*
529 * Combine the partial row results into YR
530 */
531  if( YRsum && ( Anq > 0 ) )
532  {
533  top = *PB_Ctop( &ctxt, COMBINE, COLUMN, TOP_GET );
534  Csgsum2d( ctxt, COLUMN, &top, 1, Anq, YR, YRd[LLD_], YRd[RSRC_],
535  mycol );
536  }
537 /*
538 * YR := YR + YC
539 */
540  PB_Cpaxpby( utyp, NOCONJG, *N, 1, one, YC, 0, 0, YCd, COLUMN, one,
541  YR, 0, 0, YRd, ROW );
542 /*
543 * sub( Y ) := beta * sub( Y ) + YR (if necessary)
544 */
545  if( YRpbY )
546  {
547 /*
548 * Retrieve sub( Y )'s local information: Yii, Yjj, Yrow, Ycol
549 */
550  PB_Cinfog2l( Yi, Yj, Yd, nprow, npcol, myrow, mycol, &Yii, &Yjj, &Yrow,
551  &Ycol );
552 /*
553 * sub( Y ) resides in (a) process row(s)
554 */
555  if( ( myrow == Yrow ) || Yrow < 0 )
556  {
557 /*
558 * Make sure I own some data and scale sub( Y )
559 */
560  Ynq = PB_Cnumroc( *N, Yj, Yd[INB_], Yd[NB_], mycol, Yd[CSRC_],
561  npcol );
562  if( Ynq > 0 )
563  {
564  Yld = Yd[LLD_];
565  sascal_( &Ynq, ((char *) BETA), Mptr( ((char *) Y), Yii,
566  Yjj, Yld, usiz ), &Yld );
567  }
568  }
569  PB_Cpaxpby( utyp, NOCONJG, 1, *N, one, YR, 0, 0, YRd, ROW, one,
570  ((char *) Y), Yi, Yj, Yd, ROW );
571  }
572  }
573  else
574  {
575 /*
576 * Combine the partial row results into YR
577 */
578  if( YRsum )
579  {
580  YRd[RSRC_] = 0;
581  if( Anq > 0 )
582  {
583  top = *PB_Ctop( &ctxt, COMBINE, COLUMN, TOP_GET );
584  Csgsum2d( ctxt, COLUMN, &top, 1, Anq, YR, YRd[LLD_], 0,
585  mycol );
586  }
587  }
588 /*
589 * Combine the partial column results into YC
590 */
591  if( YCsum && ( Amp > 0 ) )
592  {
593  top = *PB_Ctop( &ctxt, COMBINE, ROW, TOP_GET );
594  Csgsum2d( ctxt, ROW, &top, Amp, 1, YC, YCd[LLD_], myrow,
595  YCd[CSRC_] );
596  }
597 /*
598 * YC := YR + YC
599 */
600  PB_Cpaxpby( utyp, NOCONJG, 1, *N, one, YR, 0, 0, YRd, ROW, one,
601  YC, 0, 0, YCd, COLUMN );
602 /*
603 * sub( Y ) := beta * sub( Y ) + YC (if necessary)
604 */
605  if( YCpbY )
606  {
607 /*
608 * Retrieve sub( Y )'s local information: Yii, Yjj, Yrow, Ycol
609 */
610  PB_Cinfog2l( Yi, Yj, Yd, nprow, npcol, myrow, mycol, &Yii, &Yjj, &Yrow,
611  &Ycol );
612 /*
613 * sub( Y ) resides in (a) process column(s)
614 */
615  if( ( mycol == Ycol ) || Ycol < 0 )
616  {
617 /*
618 * Make sure I own some data and scale sub( Y )
619 */
620  Ynp = PB_Cnumroc( *N, Yi, Yd[IMB_], Yd[MB_], myrow, Yd[RSRC_],
621  nprow );
622  if( Ynp > 0 )
623  {
624  sascal_( &Ynp, ((char *) BETA), Mptr( ((char *) Y), Yii,
625  Yjj, Yd[LLD_], usiz ), INCY );
626  }
627  }
628  PB_Cpaxpby( utyp, NOCONJG, *N, 1, one, YC, 0, 0, YCd, COLUMN, one,
629  ((char *) Y), Yi, Yj, Yd, COLUMN );
630  }
631  }
632  if( YCfr ) free( YC );
633  if( YRfr ) free( YR );
634 /*
635 * End of PCAHEMV
636 */
637 }
M_
#define M_
Definition: PBtools.h:39
ROW
#define ROW
Definition: PBblacs.h:46
MB_
#define MB_
Definition: PBtools.h:43
PB_Cpaxpby
void PB_Cpaxpby()
PB_Cwarn
void PB_Cwarn()
NB_
#define NB_
Definition: PBtools.h:44
COLUMN
#define COLUMN
Definition: PBblacs.h:45
CSRC_
#define CSRC_
Definition: PBtools.h:46
PB_Cpsym
void PB_Cpsym()
PBblacs.h
PBtools.h
PBblas.h
NOCONJG
#define NOCONJG
Definition: PBblas.h:45
REAL_PART
#define REAL_PART
Definition: pblas.h:135
PBTYP_T::type
char type
Definition: pblas.h:327
PBpblas.h
DLEN_
#define DLEN_
Definition: PBtools.h:48
cagemv_
F_VOID_FCT cagemv_()
PBTYP_T::usiz
int usiz
Definition: pblas.h:328
NOTRAN
#define NOTRAN
Definition: PBblas.h:44
LLD_
#define LLD_
Definition: PBtools.h:47
PB_Cdescribe
void PB_Cdescribe()
F_CHAR_T
char * F_CHAR_T
Definition: pblas.h:118
ZERO
#define ZERO
Definition: PBtools.h:66
PB_Cchkvec
void PB_Cchkvec()
UPPER
#define UPPER
Definition: PBblas.h:52
IMB_
#define IMB_
Definition: PBtools.h:41
pilaenv_
int pilaenv_()
INIT
#define INIT
Definition: PBblas.h:61
PB_Cabort
void PB_Cabort()
CLOWER
#define CLOWER
Definition: PBblas.h:25
LEFT
#define LEFT
Definition: PBblas.h:55
F2C_CHAR
#define F2C_CHAR(a)
Definition: pblas.h:120
TOP_GET
#define TOP_GET
Definition: PBblacs.h:50
PB_Ctop
char * PB_Ctop()
ONE
#define ONE
Definition: PBtools.h:64
RSRC_
#define RSRC_
Definition: PBtools.h:45
PBTYP_T::one
char * one
Definition: pblas.h:331
PB_CargFtoC
void PB_CargFtoC()
COMBINE
#define COMBINE
Definition: PBblacs.h:49
PBTYP_T::size
int size
Definition: pblas.h:329
PB_Cinfog2l
void PB_Cinfog2l()
COTRAN
#define COTRAN
Definition: PBblas.h:48
PB_Cchkmat
void PB_Cchkmat()
PB_Cnumroc
int PB_Cnumroc()
PB_Cstypeset
PBTYP_T * PB_Cstypeset()
Definition: PB_Cstypeset.c:19
PB_CInV
void PB_CInV()
sascal_
F_VOID_FCT sascal_()
PB_CInOutV
void PB_CInOutV()
MIN
#define MIN(a_, b_)
Definition: PBtools.h:76
INB_
#define INB_
Definition: PBtools.h:42
LOWER
#define LOWER
Definition: PBblas.h:51
PB_COutV
void PB_COutV()
C2F_CHAR
#define C2F_CHAR(a)
Definition: pblas.h:121
Cblacs_gridinfo
void Cblacs_gridinfo()
PBTYP_T
Definition: pblas.h:325
Mupcase
#define Mupcase(C)
Definition: PBtools.h:83
pblas.h
CUPPER
#define CUPPER
Definition: PBblas.h:26
pcahemv_
void pcahemv_(F_CHAR_T UPLO, int *N, float *ALPHA, float *A, int *IA, int *JA, int *DESCA, float *X, int *IX, int *JX, int *DESCX, int *INCX, float *BETA, float *Y, int *IY, int *JY, int *DESCY, int *INCY)
Definition: pcahemv_.c:26
Mptr
#define Mptr(a_, i_, j_, lda_, siz_)
Definition: PBtools.h:132
CTXT_
#define CTXT_
Definition: PBtools.h:38
PB_Ctzahemv
void PB_Ctzahemv()
PB_Cctypeset
PBTYP_T * PB_Cctypeset()
Definition: PB_Cctypeset.c:19
Csgsum2d
void Csgsum2d()
PB_Clcm
int PB_Clcm()