SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pzrot.c
Go to the documentation of this file.
1/* ---------------------------------------------------------------------
2*
3* -- Mark R. Fahey
4* June 28, 2000
5*
6* ---------------------------------------------------------------------
7*/
8/*
9* Include files
10*/
11#include "pblas.h"
12
13void pzrot_( Int *n, complex16 X[], Int *ix, Int *jx, Int desc_X[], Int *incx, complex16 Y[], Int *iy, Int *jy, Int desc_Y[], Int *incy, double *c, complex16 *s )
14/*
15* Mark Fahey
16* June 22, 2000
17*/
18{
19/*
20* Purpose
21* =======
22*
23* PZROT applies a plane rotation, where the cos (C) is real and the
24* sin (S) is complex, and the vectors CX and CY are complex, i.e.,
25*
26* [ sub( X ) ] := [ C S ] [ sub( X ) ]
27* [ sub( Y ) ] := [ -conjg(S) C ] [ sub( Y ) ]
28*
29* where sub( X ) denotes X(IX,JX:JX+N-1) if INCX = M_X,
30* X(IX:IX+N-1,JX) if INCX = 1 and INCX <> M_X,
31*
32* sub( Y ) denotes Y(IY,JY:JY+N-1) if INCY = M_Y,
33* Y(IY:IY+N-1,JY) if INCY = 1 and INCY <> M_Y,
34*
35* and where C*C + S*CONJG(S) = 1.0.
36*
37* Notes
38* =====
39*
40* Each global data object is described by an associated description
41* vector. This vector stores the information required to establish
42* the mapping between an object element and its corresponding process
43* and memory location.
44*
45* Let A be a generic term for any 2D block cyclicly distributed array.
46* Such a global array has an associated description vector DESCA.
47* In the following comments, the character _ should be read as
48* "of the global array".
49*
50* NOTATION STORED IN EXPLANATION
51* --------------- -------------- --------------------------------------
52* DT_A (global) descA[ DT_ ] The descriptor type. In this case,
53* DT_A = 1.
54* CTXT_A (global) descA[ CTXT_ ] The BLACS context handle, indicating
55* the BLACS process grid A is distribu-
56* ted over. The context itself is glo-
57* bal, but the handle (the integer
58* value) may vary.
59* M_A (global) descA[ M_ ] The number of rows in the global
60* array A.
61* N_A (global) descA[ N_ ] The number of columns in the global
62* array A.
63* MB_A (global) descA[ MB_ ] The blocking factor used to distribu-
64* te the rows of the array.
65* NB_A (global) descA[ NB_ ] The blocking factor used to distribu-
66* te the columns of the array.
67* RSRC_A (global) descA[ RSRC_ ] The process row over which the first
68* row of the array A is distributed.
69* CSRC_A (global) descA[ CSRC_ ] The process column over which the
70* first column of the array A is
71* distributed.
72* LLD_A (local) descA[ LLD_ ] The leading dimension of the local
73* array. LLD_A >= MAX(1,LOCr(M_A)).
74*
75* Let K be the number of rows or columns of a distributed matrix,
76* and assume that its process grid has dimension p x q.
77* LOCr( K ) denotes the number of elements of K that a process
78* would receive if K were distributed over the p processes of its
79* process column.
80* Similarly, LOCc( K ) denotes the number of elements of K that a
81* process would receive if K were distributed over the q processes of
82* its process row.
83* The values of LOCr() and LOCc() may be determined via a call to the
84* ScaLAPACK tool function, NUMROC:
85* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
86* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
87* An upper bound for these quantities may be computed by:
88* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
89* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
90*
91* Because vectors may be seen as particular matrices, a distributed
92* vector is considered to be a distributed matrix.
93*
94* If INCX = M_X and INCY = M_Y, NB_X must be equal to NB_Y, and the
95* process column having the first entries of sub( Y ) must also contain
96* the first entries of sub( X ). Moreover, the quantity
97* MOD( JX-1, NB_X ) must be equal to MOD( JY-1, NB_Y ).
98*
99* If INCX = M_X, INCY = 1 and INCY <> M_Y, NB_X must be equal to MB_Y.
100* Moreover, the quantity MOD( JX-1, NB_X ) must be equal to
101* MOD( IY-1, MB_Y ).
102*
103* If INCX = 1, INCX <> M_X and INCY = M_Y, MB_X must be equal to NB_Y.
104* Moreover, the quantity MOD( IX-1, MB_X ) must be equal to
105* MOD( JY-1, NB_Y ).
106*
107* If INCX = 1, INCX <> M_X, INCY = 1 and INCY <> M_Y, MB_X must be
108* equal to MB_Y, and the process row having the first entries of
109* sub( Y ) must also contain the first entries of sub( X ). Moreover,
110* the quantity MOD( IX-1, MB_X ) must be equal to MOD( IY-1, MB_Y ).
111*
112* Arguments
113* =========
114*
115* N (input) INTEGER
116* The number of elements in the vectors CX and CY.
117*
118* X (local input) COMPLEX array containing the local
119* pieces of a distributed matrix of dimension of at least
120* ( (JX-1)*M_X + IX + ( N - 1 )*abs( INCX ) )
121* This array contains the entries of the distributed vector
122* sub( X ).
123* On output, CX is overwritten with C*X + S*Y.
124*
125* IX (global input) pointer to INTEGER
126* The global row index of the submatrix of the distributed
127* matrix X to operate on.
128*
129* JX (global input) pointer to INTEGER
130* The global column index of the submatrix of the distributed
131* matrix X to operate on.
132*
133* DESCX (global and local input) INTEGER array of dimension 8.
134* The array descriptor of the distributed matrix X.
135*
136* INCX (global input) pointer to INTEGER
137* The global increment for the elements of X. Only two values
138* of INCX are supported in this version, namely 1 and M_X.
139*
140* Y (local input) COMPLEX array containing the local
141* pieces of a distributed matrix of dimension of at least
142* ( (JY-1)*M_Y + IY + ( N - 1 )*abs( INCY ) )
143* This array contains the entries of the distributed vector
144* sub( Y ).
145* On output, CY is overwritten with -CONJG(S)*X + C*Y.
146*
147* IY (global input) pointer to INTEGER
148* The global row index of the submatrix of the distributed
149* matrix Y to operate on.
150*
151* JY (global input) pointer to INTEGER
152* The global column index of the submatrix of the distributed
153* matrix Y to operate on.
154*
155* DESCY (global and local input) INTEGER array of dimension 8.
156* The array descriptor of the distributed matrix Y.
157*
158* INCY (global input) pointer to INTEGER
159* The global increment for the elements of Y. Only two values
160* of INCY are supported in this version, namely 1 and M_Y.
161*
162* C (input) pointer to DOUBLE
163* S (input) pointer COMPLEX
164* C and S define a rotation
165* [ C S ]
166* [ -conjg(S) C ]
167* where C*C + S*CONJG(S) = 1.0.
168*
169* =====================================================================
170*
171* .. Local Scalars ..
172*/
173 Int ictxt, iix, iiy, info, ixcol, ixrow, iycol, iyrow, jjx,
174 jjy, lcm, lcmp, mycol, myrow, nn, np, np0,
175 nprow, npcol, nq, nz, ione=1, tmp1, wksz;
176 complex16 xwork[1], ywork[1], zero;
177/* ..
178* .. PBLAS Buffer ..
179*/
180 complex16 * buff;
181/* ..
182* .. External Functions ..
183*/
184 void blacs_gridinfo_();
185 void zgerv2d_();
186 void zgesd2d_();
187 void pbchkvect();
188 void PB_Cabort();
189 char * getpbbuf();
190 F_INTG_FCT pbztrnv_();
191 F_INTG_FCT zrot_();
192 F_INTG_FCT ilcm_();
193/* ..
194* .. Executable Statements ..
195*
196* Get grid parameters
197*/
198 ictxt = desc_X[CTXT_];
199 blacs_gridinfo_( &ictxt, &nprow, &npcol, &myrow, &mycol );
200/*
201* Test the input parameters
202*/
203 info = 0;
204 if( nprow == -1 )
205 info = -(500+CTXT_+1);
206 else
207 {
208 pbchkvect( *n, 1, *ix, *jx, desc_X, *incx, 5, &iix, &jjx,
209 &ixrow, &ixcol, nprow, npcol, myrow, mycol, &info );
210 pbchkvect( *n, 1, *iy, *jy, desc_Y, *incy, 10, &iiy, &jjy,
211 &iyrow, &iycol, nprow, npcol, myrow, mycol, &info );
212
213 if( info == 0 )
214 {
215 if( *n != 1 )
216 {
217 if( *incx == desc_X[M_] )
218 { /* X is distributed along a process row */
219 if( *incy == desc_Y[M_] )
220 { /* Y is distributed over a process row */
221 if( ( ixcol != iycol ) ||
222 ( ( (*jx-1) % desc_X[NB_] ) !=
223 ( (*jy-1) % desc_Y[NB_] ) ) )
224 info = -9;
225 else if( desc_Y[NB_] != desc_X[NB_] )
226 info = -(1000+NB_+1);
227 }
228 else if( ( *incy == 1 ) && ( *incy != desc_Y[M_] ) )
229 { /* Y is distributed over a process column */
230 if( ( (*jx-1) % desc_X[NB_] ) != ( (*iy-1) % desc_Y[MB_] ) )
231 info = -8;
232 else if( desc_Y[MB_] != desc_X[NB_] )
233 info = -(1000+MB_+1);
234 }
235 else
236 {
237 info = -11;
238 }
239 }
240 else if( ( *incx == 1 ) && ( *incx != desc_X[M_] ) )
241 { /* X is distributed along a process column */
242 if( *incy == desc_Y[M_] )
243 { /* Y is distributed over a process row */
244 if( ( (*ix-1) % desc_X[MB_] ) != ( (*jy-1) % desc_Y[NB_] ) )
245 info = -9;
246 else if( desc_Y[NB_] != desc_X[MB_] )
247 info = -(1000+NB_+1);
248 }
249 else if( ( *incy == 1 ) && ( *incy != desc_Y[M_] ) )
250 { /* Y is distributed over a process column */
251 if( ( ixrow != iyrow ) ||
252 ( ( (*ix-1) % desc_X[MB_] ) !=
253 ( (*iy-1) % desc_Y[MB_] ) ) )
254 info = -8;
255 else if( desc_Y[MB_] != desc_X[MB_] )
256 info = -(1000+MB_+1);
257 }
258 else
259 {
260 info = -11;
261 }
262 }
263 else
264 {
265 info = -6;
266 }
267 }
268 if( ictxt != desc_Y[CTXT_] )
269 info = -(1000+CTXT_+1);
270 }
271 }
272 if( info ) { PB_Cabort( ictxt, "PZROT", info ); return; }
273/*
274 if( info )
275 {
276 pberror_( &ictxt, "PZROT", &info );
277 return;
278 }
279*/
280
281/*
282* Quick return if possible.
283*/
284 zero.re = ZERO;
285 zero.im = ZERO;
286 if( *n == 0 ) return;
287/*
288* rotation
289*/
290 if( *n == 1 )
291 {
292 if( ( myrow == ixrow ) && ( mycol == ixcol ) )
293 {
294 buff = &X[iix-1+(jjx-1)*desc_X[LLD_]];
295 if( ( myrow != iyrow ) || ( mycol != iycol ) )
296 {
297 zgesd2d_( &ictxt, n, n, buff, n, &iyrow, &iycol );
298 zgerv2d_( &ictxt, n, n, ywork, n, &iyrow, &iycol );
299 }
300 else
301 *ywork = Y[iiy-1+(jjy-1)*desc_Y[LLD_]];
302 zrot_( n, buff, n, ywork, n, c, s );
303 X[iix-1+(jjx-1)*desc_X[LLD_]] = *buff;
304 if( ( myrow == iyrow ) && ( mycol == iycol ) )
305 Y[iiy-1+(jjy-1)*desc_Y[LLD_]] = *ywork;
306 }
307 else if( ( myrow == iyrow ) && ( mycol == iycol ) )
308 {
309 zgesd2d_( &ictxt, n, n, &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], n,
310 &ixrow, &ixcol );
311 zgerv2d_( &ictxt, n, n, xwork, n, &ixrow, &ixcol );
312 zrot_( n, xwork, n, &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], n, c, s );
313 }
314 return;
315 }
316
317 if( ( *incx == desc_X[M_] ) && ( *incy == desc_Y[M_] ) )
318 { /* X and Y are both distributed over a process row */
319 nz = (*jx-1) % desc_Y[NB_];
320 nn = *n + nz;
321 nq = numroc_( &nn, &desc_X[NB_], &mycol, &ixcol, &npcol );
322 if( mycol == ixcol )
323 nq -= nz;
324 if( ixrow == iyrow )
325 {
326 if( myrow == ixrow )
327 {
328 zrot_( &nq, &X[iix-1+(jjx-1)*desc_X[LLD_]], &desc_X[LLD_],
329 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], &desc_Y[LLD_], c, s );
330 }
331 }
332 else
333 {
334 if( myrow == ixrow )
335 {
336 zgesd2d_( &ictxt, &ione, &nq,
337 &X[iix-1+(jjx-1)*desc_X[LLD_]], &desc_X[LLD_],
338 &iyrow, &mycol );
339 buff = (complex16 *)getpbbuf( "PZROT", nq*sizeof(complex16) );
340 zgerv2d_( &ictxt, &nq, &ione, buff, &nq, &iyrow, &mycol );
341 zrot_( &nq, &X[iix-1+(jjx-1)*desc_X[LLD_]], &desc_X[LLD_],
342 buff, &ione, c, s );
343 }
344 else if( myrow == iyrow )
345 {
346 zgesd2d_( &ictxt, &ione, &nq,
347 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], &desc_Y[LLD_],
348 &ixrow, &mycol );
349 buff = (complex16 *)getpbbuf( "PZROT", nq*sizeof(complex16) );
350 zgerv2d_( &ictxt, &nq, &ione, buff, &nq, &ixrow, &mycol );
351 zrot_( &nq, buff, &ione,
352 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], &desc_Y[LLD_], c, s );
353 }
354 }
355 }
356 else if( ( *incx == 1 ) && ( *incx != desc_X[M_] ) &&
357 ( *incy == 1 ) && ( *incy != desc_Y[M_] ) )
358 { /* X and Y are both distributed over a process column */
359 nz = (*ix-1) % desc_X[MB_];
360 nn = *n + nz;
361 np = numroc_( &nn, &desc_X[MB_], &myrow, &ixrow, &nprow );
362 if( myrow == ixrow )
363 np -= nz;
364 if( ixcol == iycol )
365 {
366 if( mycol == ixcol )
367 {
368 zrot_( &np, &X[iix-1+(jjx-1)*desc_X[LLD_]], incx,
369 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], incy, c, s );
370 }
371 }
372 else
373 {
374 if( mycol == ixcol )
375 {
376 zgesd2d_( &ictxt, &np, &ione,
377 &X[iix-1+(jjx-1)*desc_X[LLD_]], &desc_X[LLD_],
378 &myrow, &iycol );
379 buff = (complex16 *)getpbbuf( "PZROT", np*sizeof(complex16) );
380 zgerv2d_( &ictxt, &np, &ione, buff, &np, &myrow, &iycol );
381 zrot_( &np, &X[iix-1+(jjx-1)*desc_X[LLD_]], incx,
382 buff, &ione, c, s );
383 }
384 else if( mycol == iycol )
385 {
386 zgesd2d_( &ictxt, &np, &ione,
387 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], &desc_Y[LLD_],
388 &myrow, &ixcol );
389 buff = (complex16 *)getpbbuf( "PZROT", np*sizeof(complex16) );
390 zgerv2d_( &ictxt, &np, &ione, buff, &np, &myrow, &ixcol );
391 zrot_( &np, buff, &ione,
392 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], incy, c, s );
393 }
394 }
395 }
396 else /* X and Y are not distributed along the same direction */
397 {
398 lcm = ilcm_( &nprow, &npcol );
399 if( ( *incx == 1 ) && ( *incx != desc_X[M_] ) )
400 { /* X is distributed over a process column */
401 lcmp = lcm / nprow;
402 nz = (*jy-1) % desc_Y[NB_];
403 nn = *n + nz;
404 tmp1 = nn / desc_Y[MB_];
405 np = numroc_( &nn, &desc_X[MB_], &myrow, &ixrow, &nprow );
406 np0 = MYROC0( tmp1, nn, desc_X[MB_], nprow );
407 tmp1 = np0 / desc_X[MB_];
408 wksz = MYROC0( tmp1, np0, desc_X[MB_], lcmp );
409 wksz = np + wksz;
410
411 buff = (complex16 *)getpbbuf( "PZROT", wksz*sizeof(complex16) );
412
413 if( mycol == iycol )
414 jjy -= nz;
415 if( myrow == ixrow )
416 np -= nz;
417 pbztrnv_( &ictxt, C2F_CHAR( "R" ), C2F_CHAR( "T" ), n,
418 &desc_Y[NB_], &nz, &Y[iiy-1+(jjy-1)*desc_Y[LLD_]],
419 &desc_Y[LLD_], &zero, buff, &ione, &iyrow, &iycol,
420 &ixrow, &ixcol, buff+np );
421 if( mycol == ixcol )
422 {
423 zrot_( &np, &X[iix-1+(jjx-1)*desc_X[LLD_]],
424 incx, buff, &ione, c, s );
425 }
426 pbztrnv_( &ictxt, C2F_CHAR( "R" ), C2F_CHAR( "T" ), n,
427 &desc_Y[NB_], &nz, buff, &ione, &zero,
428 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], &desc_Y[LLD_],
429 &ixrow, &ixcol, &iyrow, &iycol, buff+np );
430 }
431 else /* Y is distributed over a process column */
432 {
433 lcmp = lcm / nprow;
434 nz = (*jx-1) % desc_X[NB_];
435 nn = *n + nz;
436 tmp1 = nn / desc_X[MB_];
437 np = numroc_( &nn, desc_Y+MB_, &myrow, &iyrow, &nprow );
438 np0 = MYROC0( tmp1, nn, desc_Y[MB_], nprow );
439 tmp1 = np0 / desc_Y[MB_];
440 wksz = MYROC0( tmp1, np0, desc_Y[MB_], lcmp );
441 wksz = np + wksz;
442
443 buff = (complex16 *)getpbbuf( "PZROT", wksz*sizeof(complex16) );
444
445 if( myrow == iyrow )
446 np -= nz;
447 pbztrnv_( &ictxt, C2F_CHAR( "R" ), C2F_CHAR( "T" ), n,
448 &desc_X[NB_], &nz, &X[iix-1+(jjx-1)*desc_X[LLD_]],
449 &desc_X[LLD_], &zero, buff, &ione, &ixrow, &ixcol,
450 &iyrow, &iycol, buff+np );
451 if( mycol == iycol )
452 {
453 zrot_( &np, buff, &ione,
454 &Y[iiy-1+(jjy-1)*desc_Y[LLD_]], incy, c, s );
455 }
456 pbztrnv_( &ictxt, C2F_CHAR( "R" ), C2F_CHAR( "T" ), n,
457 &desc_X[NB_], &nz, buff, &ione, &zero,
458 &X[iix-1+(jjx-1)*desc_X[LLD_]], &desc_X[LLD_],
459 &iyrow, &iycol, &ixrow, &ixcol, buff+np );
460 }
461 }
462}
Int
#define Int
Definition Bconfig.h:22
F_INTG_FCT
#define F_INTG_FCT
Definition pblas.h:128
C2F_CHAR
#define C2F_CHAR(a)
Definition pblas.h:125
CTXT_
#define CTXT_
Definition PBtools.h:38
MB_
#define MB_
Definition PBtools.h:43
PB_Cabort
void PB_Cabort()
LLD_
#define LLD_
Definition PBtools.h:47
M_
#define M_
Definition PBtools.h:39
ZERO
#define ZERO
Definition PBtools.h:66
NB_
#define NB_
Definition PBtools.h:44
MYROC0
#define MYROC0(nblocks, n, nb, nprocs)
Definition pblas.h:195
blacs_gridinfo_
F_VOID_FUNC blacs_gridinfo_(Int *ConTxt, Int *nprow, Int *npcol, Int *myrow, Int *mycol)
Definition blacs_info_.c:6
getpbbuf
char * getpbbuf(char *mess, Int length)
Definition getpbbuf.c:3
pbchkvect
void pbchkvect(Int n, Int npos0, Int ix, Int jx, Int desc_X[], Int incx, Int dpos0, Int *iix, Int *jjx, Int *ixrow, Int *ixcol, Int nprow, Int npcol, Int myrow, Int mycol, Int *info)
Definition pbchkvect.c:15
pzrot_
void pzrot_(Int *n, complex16 X[], Int *ix, Int *jx, Int desc_X[], Int *incx, complex16 Y[], Int *iy, Int *jy, Int desc_Y[], Int *incy, double *c, complex16 *s)
Definition pzrot.c:13
complex16::re
double re
Definition pblas.h:97
complex16::im
double im
Definition pblas.h:97
zgerv2d_
F_VOID_FUNC zgerv2d_(Int *ConTxt, Int *m, Int *n, double *A, Int *lda, Int *rsrc, Int *csrc)
Definition zgerv2d_.c:6
zgesd2d_
F_VOID_FUNC zgesd2d_(Int *ConTxt, Int *m, Int *n, double *A, Int *lda, Int *rdest, Int *cdest)
Definition zgesd2d_.c:7