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__
pcher_(F_CHAR_T UPLO,int * N,float * ALPHA,float * X,int * IX,int * JX,int * DESCX,int * INCX,float * A,int * IA,int * JA,int * DESCA)20 void pcher_( F_CHAR_T UPLO, int * N, float * ALPHA,
21 float * X, int * IX, int * JX, int * DESCX, int * INCX,
22 float * A, int * IA, int * JA, int * DESCA )
23 #else
24 void pcher_( UPLO, N, ALPHA, X, IX, JX, DESCX, INCX, A, IA, JA, DESCA )
25 /*
26 * .. Scalar Arguments ..
27 */
28 F_CHAR_T UPLO;
29 int * IA, * INCX, * IX, * JA, * JX, * N;
30 float * ALPHA;
31 /*
32 * .. Array Arguments ..
33 */
34 int * DESCA, * DESCX;
35 float * A, * X;
36 #endif
37 {
38 /*
39 * Purpose
40 * =======
41 *
42 * PCHER performs the Hermitian rank 1 operation
43 *
44 * sub( A ) := alpha*sub( X )*conjg( sub( X )' ) + sub( A ),
45 *
46 * where
47 *
48 * sub( A ) denotes A(IA:IA+N-1,JA:JA+N-1), and,
49 *
50 * sub( X ) denotes X(IX,JX:JX+N-1) if INCX = M_X,
51 * X(IX:IX+N-1,JX) if INCX = 1 and INCX <> M_X.
52 *
53 * Alpha is a real scalar, sub( X ) is an n element subvector and
54 * sub( A ) is an n by n Hermitian submatrix.
55 *
56 * Notes
57 * =====
58 *
59 * A description vector is associated with each 2D block-cyclicly dis-
60 * tributed matrix. This vector stores the information required to
61 * establish the mapping between a matrix entry and its corresponding
62 * process and memory location.
63 *
64 * In the following comments, the character _ should be read as
65 * "of the distributed matrix". Let A be a generic term for any 2D
66 * block cyclicly distributed matrix. Its description vector is DESC_A:
67 *
68 * NOTATION STORED IN EXPLANATION
69 * ---------------- --------------- ------------------------------------
70 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
71 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
72 * the NPROW x NPCOL BLACS process grid
73 * A is distributed over. The context
74 * itself is global, but the handle
75 * (the integer value) may vary.
76 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
77 * ted matrix A, M_A >= 0.
78 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
79 * buted matrix A, N_A >= 0.
80 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
81 * block of the matrix A, IMB_A > 0.
82 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
83 * left block of the matrix A,
84 * INB_A > 0.
85 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
86 * bute the last M_A-IMB_A rows of A,
87 * MB_A > 0.
88 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
89 * bute the last N_A-INB_A columns of
90 * A, NB_A > 0.
91 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
92 * row of the matrix A is distributed,
93 * NPROW > RSRC_A >= 0.
94 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
95 * first column of A is distributed.
96 * NPCOL > CSRC_A >= 0.
97 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
98 * array storing the local blocks of
99 * the distributed matrix A,
100 * IF( Lc( 1, N_A ) > 0 )
101 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
102 * ELSE
103 * LLD_A >= 1.
104 *
105 * Let K be the number of rows of a matrix A starting at the global in-
106 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
107 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
108 * receive if these K rows were distributed over NPROW processes. If K
109 * is the number of columns of a matrix A starting at the global index
110 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
111 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
112 * these K columns were distributed over NPCOL processes.
113 *
114 * The values of Lr() and Lc() may be determined via a call to the func-
115 * tion PB_Cnumroc:
116 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
117 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
118 *
119 * Arguments
120 * =========
121 *
122 * UPLO (global input) CHARACTER*1
123 * On entry, UPLO specifies whether the local pieces of
124 * the array A containing the upper or lower triangular part
125 * of the Hermitian submatrix sub( A ) are to be referenced as
126 * follows:
127 *
128 * UPLO = 'U' or 'u' Only the local pieces corresponding to
129 * the upper triangular part of the
130 * Hermitian submatrix sub( A ) are to be
131 * referenced,
132 *
133 * UPLO = 'L' or 'l' Only the local pieces corresponding to
134 * the lower triangular part of the
135 * Hermitian submatrix sub( A ) are to be
136 * referenced.
137 *
138 * N (global input) INTEGER
139 * On entry, N specifies the order of the submatrix sub( A ).
140 * N must be at least zero.
141 *
142 * ALPHA (global input) REAL
143 * On entry, ALPHA specifies the scalar alpha. When ALPHA is
144 * supplied as zero then the local entries of the array X
145 * corresponding to the entries of the subvector sub( X ) need
146 * not be set on input.
147 *
148 * X (local input) COMPLEX array
149 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
150 * is at least MAX( 1, Lr( 1, IX ) ) when INCX = M_X and
151 * MAX( 1, Lr( 1, IX+N-1 ) ) otherwise, and, Kx is at least
152 * Lc( 1, JX+N-1 ) when INCX = M_X and Lc( 1, JX ) otherwise.
153 * Before entry, this array contains the local entries of the
154 * matrix X.
155 *
156 * IX (global input) INTEGER
157 * On entry, IX specifies X's global row index, which points to
158 * the beginning of the submatrix sub( X ).
159 *
160 * JX (global input) INTEGER
161 * On entry, JX specifies X's global column index, which points
162 * to the beginning of the submatrix sub( X ).
163 *
164 * DESCX (global and local input) INTEGER array
165 * On entry, DESCX is an integer array of dimension DLEN_. This
166 * is the array descriptor for the matrix X.
167 *
168 * INCX (global input) INTEGER
169 * On entry, INCX specifies the global increment for the
170 * elements of X. Only two values of INCX are supported in
171 * this version, namely 1 and M_X. INCX must not be zero.
172 *
173 * A (local input/local output) COMPLEX array
174 * On entry, A is an array of dimension (LLD_A, Ka), where Ka is
175 * at least Lc( 1, JA+N-1 ). Before entry, this array contains
176 * the local entries of the matrix A.
177 * Before entry with UPLO = 'U' or 'u', this array contains
178 * the local entries corresponding to the upper triangular part
179 * of the Hermitian submatrix sub( A ), and the local entries
180 * corresponding to the strictly lower triangular of sub( A )
181 * are not referenced. On exit, the upper triangular part of
182 * sub( A ) is overwritten by the upper triangular part of the
183 * updated submatrix.
184 * Before entry with UPLO = 'L' or 'l', this array contains
185 * the local entries corresponding to the lower triangular part
186 * of the Hermitian submatrix sub( A ), and the local entries
187 * corresponding to the strictly upper triangular of sub( A )
188 * are not referenced. On exit, the lower triangular part of
189 * sub( A ) is overwritten by the lower triangular part of the
190 * updated submatrix.
191 * Note that the imaginary parts of the local entries corres-
192 * ponding to the diagonal elements of sub( A ) need not be
193 * set, they are assumed to be zero, and on exit they are set
194 * to zero.
195 *
196 * IA (global input) INTEGER
197 * On entry, IA specifies A's global row index, which points to
198 * the beginning of the submatrix sub( A ).
199 *
200 * JA (global input) INTEGER
201 * On entry, JA specifies A's global column index, which points
202 * to the beginning of the submatrix sub( A ).
203 *
204 * DESCA (global and local input) INTEGER array
205 * On entry, DESCA is an integer array of dimension DLEN_. This
206 * is the array descriptor for the matrix A.
207 *
208 * -- Written on April 1, 1998 by
209 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
210 *
211 * ---------------------------------------------------------------------
212 */
213 /*
214 * .. Local Scalars ..
215 */
216 char UploA;
217 int Acol, Ai, Aii, Aimb1, Ainb1, Aj, Ajj, Akp, Akq, Ald, Amb,
218 Amp, Amp0, Anb, Anq, Anq0, Arow, XCfr, XCld, XRfr, XRld,
219 Xi, Xj, ctxt, info, ione=1, k, kb, ktmp, mycol, myrow, nb,
220 npcol, nprow, size, upper;
221 cmplx Calpha;
222 PBTYP_T * type;
223 /*
224 * .. Local Arrays ..
225 */
226 int Ad[DLEN_], Ad0[DLEN_], XCd0[DLEN_], XRd0[DLEN_], Xd[DLEN_];
227 char * Aptr = NULL, * XC = NULL, * XR = NULL;
228 /* ..
229 * .. Executable Statements ..
230 *
231 */
232 upper = ( ( UploA = Mupcase( F2C_CHAR( UPLO )[0] ) ) == CUPPER );
233 PB_CargFtoC( *IA, *JA, DESCA, &Ai, &Aj, Ad );
234 PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
235 #ifndef NO_ARGCHK
236 /*
237 * Test the input parameters
238 */
239 Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
240 if( !( info = ( ( nprow == -1 ) ? -( 701 + CTXT_ ) : 0 ) ) )
241 {
242 if( ( !upper ) && ( UploA != CLOWER ) )
243 {
244 PB_Cwarn( ctxt, __LINE__, "PCHER", "Illegal UPLO = %c\n", UploA );
245 info = -1;
246 }
247 PB_Cchkvec( ctxt, "PCHER", "X", *N, 2, Xi, Xj, Xd, *INCX, 7, &info );
248 PB_Cchkmat( ctxt, "PCHER", "A", *N, 2, *N, 2, Ai, Aj, Ad, 12, &info );
249 }
250 if( info ) { PB_Cabort( ctxt, "PCHER", info ); return; }
251 #endif
252 /*
253 * Quick return if possible
254 */
255 if( ( *N == 0 ) || ( ALPHA[REAL_PART] == ZERO ) )
256 return;
257 /*
258 * Retrieve process grid information
259 */
260 #ifdef NO_ARGCHK
261 Cblacs_gridinfo( ( ctxt = Ad[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
262 #endif
263 /*
264 * Get type structure
265 */
266 type = PB_Cctypeset();
267 /*
268 * Compute descriptor Ad0 for sub( A )
269 */
270 PB_Cdescribe( *N, *N, Ai, Aj, Ad, nprow, npcol, myrow, mycol, &Aii, &Ajj,
271 &Ald, &Aimb1, &Ainb1, &Amb, &Anb, &Arow, &Acol, Ad0 );
272 /*
273 * Replicate sub( X ) in process rows (XR) and process columns (XC) spanned by
274 * sub( A )
275 */
276 if( *INCX == Xd[M_] )
277 {
278 PB_CInV( type, NOCONJG, ROW, *N, *N, Ad0, 1, ((char *) X), Xi, Xj,
279 Xd, ROW, &XR, XRd0, &XRfr );
280 PB_CInV( type, NOCONJG, COLUMN, *N, *N, Ad0, 1, XR, 0, 0,
281 XRd0, ROW, &XC, XCd0, &XCfr );
282 }
283 else
284 {
285 PB_CInV( type, NOCONJG, COLUMN, *N, *N, Ad0, 1, ((char *) X), Xi, Xj,
286 Xd, COLUMN, &XC, XCd0, &XCfr );
287 PB_CInV( type, NOCONJG, ROW, *N, *N, Ad0, 1, XC, 0, 0,
288 XCd0, COLUMN, &XR, XRd0, &XRfr );
289 }
290 /*
291 * Local rank-1 update if I own some data
292 */
293 Amp = PB_Cnumroc( *N, 0, Aimb1, Amb, myrow, Arow, nprow );
294 Anq = PB_Cnumroc( *N, 0, Ainb1, Anb, mycol, Acol, npcol );
295
296 if( ( Amp > 0 ) && ( Anq > 0 ) )
297 {
298 size = type->size;
299 Aptr = Mptr( ((char *) A), Aii, Ajj, Ald, size );
300 /*
301 * Computational partitioning size is computed as the product of the logical
302 * value returned by pilaenv_ and 2 * lcm( nprow, npcol ).
303 */
304 nb = 2 * pilaenv_( &ctxt, C2F_CHAR( &type->type ) ) *
305 PB_Clcm( ( Arow >= 0 ? nprow : 1 ), ( Acol >= 0 ? npcol : 1 ) );
306
307 XCld = XCd0[LLD_]; XRld = XRd0[LLD_];
308 Calpha[REAL_PART] = ALPHA[REAL_PART];
309 Calpha[IMAG_PART] = ZERO;
310
311 if( upper )
312 {
313 for( k = 0; k < *N; k += nb )
314 {
315 kb = *N - k; kb = MIN( kb, nb );
316 Akp = PB_Cnumroc( k, 0, Aimb1, Amb, myrow, Arow, nprow );
317 Akq = PB_Cnumroc( k, 0, Ainb1, Anb, mycol, Acol, npcol );
318 Anq0 = PB_Cnumroc( kb, k, Ainb1, Anb, mycol, Acol, npcol );
319 if( Akp > 0 && Anq0 > 0 )
320 cgerc_( &Akp, &Anq0, ((char *) Calpha), XC, &ione,
321 Mptr( XR, 0, Akq, XRld, size ), &XRld, Mptr( Aptr, 0,
322 Akq, Ald, size ), &Ald );
323 PB_Cpsyr( type, UPPER, kb, 1, ((char *) Calpha), Mptr( XC, Akp, 0,
324 XCld, size ), XCld, Mptr( XR, 0, Akq, XRld, size ), XRld,
325 Aptr, k, k, Ad0, PB_Ctzher );
326 }
327 }
328 else
329 {
330 for( k = 0; k < *N; k += nb )
331 {
332 kb = *N - k; ktmp = k + ( kb = MIN( kb, nb ) );
333 Akp = PB_Cnumroc( k, 0, Aimb1, Amb, myrow, Arow, nprow );
334 Akq = PB_Cnumroc( k, 0, Ainb1, Anb, mycol, Acol, npcol );
335 PB_Cpsyr( type, LOWER, kb, 1, ((char *) Calpha), Mptr( XC, Akp, 0,
336 XCld, size ), XCld, Mptr( XR, 0, Akq, XRld, size ), XRld,
337 Aptr, k, k, Ad0, PB_Ctzher );
338 Akp = PB_Cnumroc( ktmp, 0, Aimb1, Amb, myrow, Arow, nprow );
339 Amp0 = Amp - Akp;
340 Anq0 = PB_Cnumroc( kb, k, Ainb1, Anb, mycol, Acol, npcol );
341 if( Amp0 > 0 && Anq0 > 0 )
342 cgerc_( &Amp0, &Anq0, ((char *) Calpha), Mptr( XC, Akp,
343 0, XCld, size ), &ione, Mptr( XR, 0, Akq, XRld, size ),
344 &XRld, Mptr( Aptr, Akp, Akq, Ald, size ), &Ald );
345 }
346 }
347 }
348 if( XRfr ) free( XR );
349 if( XCfr ) free( XC );
350 /*
351 * End of PCHER
352 */
353 }
354