ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pdsyr_.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 pdsyr_( 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 pdsyr_( 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 * PDSYR performs the symmetric rank 1 operation
43 *
44 * sub( A ) := alpha*sub( X )*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 scalar, sub( X ) is an n element subvector and sub( A ) is
54 * an n by n symmetric 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 symmetric 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 * symmetric 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 * symmetric 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) DOUBLE PRECISION 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) DOUBLE PRECISION 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 symmetric 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 symmetric 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 *
192 * IA (global input) INTEGER
193 * On entry, IA specifies A's global row index, which points to
194 * the beginning of the submatrix sub( A ).
195 *
196 * JA (global input) INTEGER
197 * On entry, JA specifies A's global column index, which points
198 * to the beginning of the submatrix sub( A ).
199 *
200 * DESCA (global and local input) INTEGER array
201 * On entry, DESCA is an integer array of dimension DLEN_. This
202 * is the array descriptor for the matrix A.
203 *
204 * -- Written on April 1, 1998 by
205 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
206 *
207 * ---------------------------------------------------------------------
208 */
209 /*
210 * .. Local Scalars ..
211 */
212  char UploA;
213  int Acol, Ai, Aii, Aimb1, Ainb1, Aj, Ajj, Akp, Akq, Ald, Amb,
214  Amp, Amp0, Anb, Anq, Anq0, Arow, XCfr, XCld, XRfr, XRld,
215  Xi, Xj, ctxt, info, ione=1, k, kb, ktmp, mycol, myrow, nb,
216  npcol, nprow, size, upper;
217  PBTYP_T * type;
218 /*
219 * .. Local Arrays ..
220 */
221  int Ad[DLEN_], Ad0[DLEN_], XCd0[DLEN_], XRd0[DLEN_], Xd[DLEN_];
222  char * Aptr = NULL, * XC = NULL, * XR = NULL;
223 /* ..
224 * .. Executable Statements ..
225 *
226 */
227  upper = ( ( UploA = Mupcase( F2C_CHAR( UPLO )[0] ) ) == CUPPER );
228  PB_CargFtoC( *IA, *JA, DESCA, &Ai, &Aj, Ad );
229  PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
230 #ifndef NO_ARGCHK
231 /*
232 * Test the input parameters
233 */
234  Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
235  if( !( info = ( ( nprow == -1 ) ? -( 701 + CTXT_ ) : 0 ) ) )
236  {
237  if( ( !upper ) && ( UploA != CLOWER ) )
238  {
239  PB_Cwarn( ctxt, __LINE__, "PDSYR", "Illegal UPLO = %c\n", UploA );
240  info = -1;
241  }
242  PB_Cchkvec( ctxt, "PDSYR", "X", *N, 2, Xi, Xj, Xd, *INCX, 7, &info );
243  PB_Cchkmat( ctxt, "PDSYR", "A", *N, 2, *N, 2, Ai, Aj, Ad, 12, &info );
244  }
245  if( info ) { PB_Cabort( ctxt, "PDSYR", info ); return; }
246 #endif
247 /*
248 * Quick return if possible
249 */
250  if( ( *N == 0 ) || ( ALPHA[REAL_PART] == ZERO ) )
251  return;
252 /*
253 * Retrieve process grid information
254 */
255 #ifdef NO_ARGCHK
256  Cblacs_gridinfo( ( ctxt = Ad[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
257 #endif
258 /*
259 * Get type structure
260 */
261  type = PB_Cdtypeset();
262 /*
263 * Compute descriptor Ad0 for sub( A )
264 */
265  PB_Cdescribe( *N, *N, Ai, Aj, Ad, nprow, npcol, myrow, mycol, &Aii, &Ajj,
266  &Ald, &Aimb1, &Ainb1, &Amb, &Anb, &Arow, &Acol, Ad0 );
267 /*
268 * Replicate sub( X ) in process rows (XR) and process columns (XC) spanned by
269 * sub( A )
270 */
271  if( *INCX == Xd[M_] )
272  {
273  PB_CInV( type, NOCONJG, ROW, *N, *N, Ad0, 1, ((char *) X), Xi, Xj,
274  Xd, ROW, &XR, XRd0, &XRfr );
275  PB_CInV( type, NOCONJG, COLUMN, *N, *N, Ad0, 1, XR, 0, 0,
276  XRd0, ROW, &XC, XCd0, &XCfr );
277  }
278  else
279  {
280  PB_CInV( type, NOCONJG, COLUMN, *N, *N, Ad0, 1, ((char *) X), Xi, Xj,
281  Xd, COLUMN, &XC, XCd0, &XCfr );
282  PB_CInV( type, NOCONJG, ROW, *N, *N, Ad0, 1, XC, 0, 0,
283  XCd0, COLUMN, &XR, XRd0, &XRfr );
284  }
285 /*
286 * Local rank-1 update if I own some data
287 */
288  Amp = PB_Cnumroc( *N, 0, Aimb1, Amb, myrow, Arow, nprow );
289  Anq = PB_Cnumroc( *N, 0, Ainb1, Anb, mycol, Acol, npcol );
290 
291  if( ( Amp > 0 ) && ( Anq > 0 ) )
292  {
293  size = type->size;
294  Aptr = Mptr( ((char *) A), Aii, Ajj, Ald, size );
295 /*
296 * Computational partitioning size is computed as the product of the logical
297 * value returned by pilaenv_ and 2 * lcm( nprow, npcol ).
298 */
299  nb = 2 * pilaenv_( &ctxt, C2F_CHAR( &type->type ) ) *
300  PB_Clcm( ( Arow >= 0 ? nprow : 1 ), ( Acol >= 0 ? npcol : 1 ) );
301 
302  XCld = XCd0[LLD_]; XRld = XRd0[LLD_];
303 
304  if( upper )
305  {
306  for( k = 0; k < *N; k += nb )
307  {
308  kb = *N - k; kb = MIN( kb, nb );
309  Akp = PB_Cnumroc( k, 0, Aimb1, Amb, myrow, Arow, nprow );
310  Akq = PB_Cnumroc( k, 0, Ainb1, Anb, mycol, Acol, npcol );
311  Anq0 = PB_Cnumroc( kb, k, Ainb1, Anb, mycol, Acol, npcol );
312  if( Akp > 0 && Anq0 > 0 )
313  dger_( &Akp, &Anq0, ((char *) ALPHA), XC, &ione,
314  Mptr( XR, 0, Akq, XRld, size ), &XRld, Mptr( Aptr, 0, Akq,
315  Ald, size ), &Ald );
316  PB_Cpsyr( type, UPPER, kb, 1, ((char *) ALPHA), Mptr( XC, Akp, 0,
317  XCld, size ), XCld, Mptr( XR, 0, Akq, XRld, size ), XRld,
318  Aptr, k, k, Ad0, PB_Ctzsyr );
319  }
320  }
321  else
322  {
323  for( k = 0; k < *N; k += nb )
324  {
325  kb = *N - k; ktmp = k + ( kb = MIN( kb, nb ) );
326  Akp = PB_Cnumroc( k, 0, Aimb1, Amb, myrow, Arow, nprow );
327  Akq = PB_Cnumroc( k, 0, Ainb1, Anb, mycol, Acol, npcol );
328  PB_Cpsyr( type, LOWER, kb, 1, ((char *) ALPHA), Mptr( XC, Akp, 0,
329  XCld, size ), XCld, Mptr( XR, 0, Akq, XRld, size ), XRld,
330  Aptr, k, k, Ad0, PB_Ctzsyr );
331  Akp = PB_Cnumroc( ktmp, 0, Aimb1, Amb, myrow, Arow, nprow );
332  Amp0 = Amp - Akp;
333  Anq0 = PB_Cnumroc( kb, k, Ainb1, Anb, mycol, Acol, npcol );
334  if( Amp0 > 0 && Anq0 > 0 )
335  dger_( &Amp0, &Anq0, ((char *) ALPHA), Mptr( XC, Akp,
336  0, XCld, size ), &ione, Mptr( XR, 0, Akq, XRld, size ),
337  &XRld, Mptr( Aptr, Akp, Akq, Ald, size ), &Ald );
338  }
339  }
340  }
341  if( XRfr ) free( XR );
342  if( XCfr ) free( XC );
343 /*
344 * End of PDSYR
345 */
346 }
M_
#define M_
Definition: PBtools.h:39
ROW
#define ROW
Definition: PBblacs.h:46
dger_
F_VOID_FCT dger_()
PB_Cwarn
void PB_Cwarn()
COLUMN
#define COLUMN
Definition: PBblacs.h:45
pdsyr_
void pdsyr_(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: pdsyr_.c:24
PBblacs.h
PB_Cpsyr
void PB_Cpsyr()
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
LLD_
#define LLD_
Definition: PBtools.h:47
PB_Cdescribe
void PB_Cdescribe()
PB_Cdtypeset
PBTYP_T * PB_Cdtypeset()
Definition: PB_Cdtypeset.c:19
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_()
PB_Cabort
void PB_Cabort()
CLOWER
#define CLOWER
Definition: PBblas.h:25
F2C_CHAR
#define F2C_CHAR(a)
Definition: pblas.h:120
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
PB_Ctzsyr
void PB_Ctzsyr()
Mptr
#define Mptr(a_, i_, j_, lda_, siz_)
Definition: PBtools.h:132
CTXT_
#define CTXT_
Definition: PBtools.h:38
PB_Clcm
int PB_Clcm()