ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
pzher_.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 pzher_( F_CHAR_T UPLO, int * N, double * ALPHA,
21  double * X, int * IX, int * JX, int * DESCX, int * INCX,
22  double * A, int * IA, int * JA, int * DESCA )
23 #else
24 void pzher_( 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  double * ALPHA;
31 /*
32 * .. Array Arguments ..
33 */
34  int * DESCA, * DESCX;
35  double * A, * X;
36 #endif
37 {
38 /*
39 * Purpose
40 * =======
41 *
42 * PZHER 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) DOUBLE PRECISION
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*16 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*16 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  cmplx16 Calpha;
222  PBTYP_T * type;
223 /*
224 * .. Local Arrays ..
225 */
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__, "PZHER", "Illegal UPLO = %c\n", UploA );
245  info = -1;
246  }
247  PB_Cchkvec( ctxt, "PZHER", "X", *N, 2, Xi, Xj, Xd, *INCX, 7, &info );
248  PB_Cchkmat( ctxt, "PZHER", "A", *N, 2, *N, 2, Ai, Aj, Ad, 12, &info );
249  }
250  if( info ) { PB_Cabort( ctxt, "PZHER", 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_Cztypeset();
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  zgerc_( &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  zgerc_( &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 PZHER
352 */
353 }
M_
#define M_
Definition: PBtools.h:39
ROW
#define ROW
Definition: PBblacs.h:46
PB_Cwarn
void PB_Cwarn()
COLUMN
#define COLUMN
Definition: PBblacs.h:45
PBblacs.h
PB_Cpsyr
void PB_Cpsyr()
PBtools.h
PBblas.h
PB_Ctzher
void PB_Ctzher()
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
PB_Cztypeset
PBTYP_T * PB_Cztypeset()
Definition: PB_Cztypeset.c:19
DLEN_
#define DLEN_
Definition: PBtools.h:48
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
pilaenv_
int pilaenv_()
zgerc_
F_VOID_FCT zgerc_()
PB_Cabort
void PB_Cabort()
CLOWER
#define CLOWER
Definition: PBblas.h:25
F2C_CHAR
#define F2C_CHAR(a)
Definition: pblas.h:120
cmplx16
double cmplx16[2]
Definition: pblas.h:133
PB_CargFtoC
void PB_CargFtoC()
PBTYP_T::size
int size
Definition: pblas.h:329
PB_Cchkmat
void PB_Cchkmat()
PB_Cnumroc
int PB_Cnumroc()
PB_CInV
void PB_CInV()
MIN
#define MIN(a_, b_)
Definition: PBtools.h:76
LOWER
#define LOWER
Definition: PBblas.h:51
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
pzher_
void pzher_(F_CHAR_T UPLO, int *N, double *ALPHA, double *X, int *IX, int *JX, int *DESCX, int *INCX, double *A, int *IA, int *JA, int *DESCA)
Definition: pzher_.c:24
Mptr
#define Mptr(a_, i_, j_, lda_, siz_)
Definition: PBtools.h:132
CTXT_
#define CTXT_
Definition: PBtools.h:38
IMAG_PART
#define IMAG_PART
Definition: pblas.h:136
PB_Clcm
int PB_Clcm()