ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pscnrm2_.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 pscnrm2_( int * N, float * NORM2,
21  float * X, int * IX, int * JX, int * DESCX, int * INCX )
22 #else
23 void pscnrm2_( N, NORM2, X, IX, JX, DESCX, INCX )
24 /*
25 * .. Scalar Arguments ..
26 */
27  int * INCX, * IX, * JX, * N;
28  float * NORM2;
29 /*
30 * .. Array Arguments ..
31 */
32  int * DESCX;
33  float * X;
34 #endif
35 {
36 /*
37 * Purpose
38 * =======
39 *
40 * PSCNRM2 computes the 2-norm of a subvector sub( X ),
41 *
42 * where
43 *
44 * sub( X ) denotes X(IX,JX:JX+N-1) if INCX = M_X,
45 * X(IX:IX+N-1,JX) if INCX = 1 and INCX <> M_X.
46 *
47 * Notes
48 * =====
49 *
50 * A description vector is associated with each 2D block-cyclicly dis-
51 * tributed matrix. This vector stores the information required to
52 * establish the mapping between a matrix entry and its corresponding
53 * process and memory location.
54 *
55 * In the following comments, the character _ should be read as
56 * "of the distributed matrix". Let A be a generic term for any 2D
57 * block cyclicly distributed matrix. Its description vector is DESC_A:
58 *
59 * NOTATION STORED IN EXPLANATION
60 * ---------------- --------------- ------------------------------------
61 * DTYPE_A (global) DESCA[ DTYPE_ ] The descriptor type.
62 * CTXT_A (global) DESCA[ CTXT_ ] The BLACS context handle, indicating
63 * the NPROW x NPCOL BLACS process grid
64 * A is distributed over. The context
65 * itself is global, but the handle
66 * (the integer value) may vary.
67 * M_A (global) DESCA[ M_ ] The number of rows in the distribu-
68 * ted matrix A, M_A >= 0.
69 * N_A (global) DESCA[ N_ ] The number of columns in the distri-
70 * buted matrix A, N_A >= 0.
71 * IMB_A (global) DESCA[ IMB_ ] The number of rows of the upper left
72 * block of the matrix A, IMB_A > 0.
73 * INB_A (global) DESCA[ INB_ ] The number of columns of the upper
74 * left block of the matrix A,
75 * INB_A > 0.
76 * MB_A (global) DESCA[ MB_ ] The blocking factor used to distri-
77 * bute the last M_A-IMB_A rows of A,
78 * MB_A > 0.
79 * NB_A (global) DESCA[ NB_ ] The blocking factor used to distri-
80 * bute the last N_A-INB_A columns of
81 * A, NB_A > 0.
82 * RSRC_A (global) DESCA[ RSRC_ ] The process row over which the first
83 * row of the matrix A is distributed,
84 * NPROW > RSRC_A >= 0.
85 * CSRC_A (global) DESCA[ CSRC_ ] The process column over which the
86 * first column of A is distributed.
87 * NPCOL > CSRC_A >= 0.
88 * LLD_A (local) DESCA[ LLD_ ] The leading dimension of the local
89 * array storing the local blocks of
90 * the distributed matrix A,
91 * IF( Lc( 1, N_A ) > 0 )
92 * LLD_A >= MAX( 1, Lr( 1, M_A ) )
93 * ELSE
94 * LLD_A >= 1.
95 *
96 * Let K be the number of rows of a matrix A starting at the global in-
97 * dex IA,i.e, A( IA:IA+K-1, : ). Lr( IA, K ) denotes the number of rows
98 * that the process of row coordinate MYROW ( 0 <= MYROW < NPROW ) would
99 * receive if these K rows were distributed over NPROW processes. If K
100 * is the number of columns of a matrix A starting at the global index
101 * JA, i.e, A( :, JA:JA+K-1, : ), Lc( JA, K ) denotes the number of co-
102 * lumns that the process MYCOL ( 0 <= MYCOL < NPCOL ) would receive if
103 * these K columns were distributed over NPCOL processes.
104 *
105 * The values of Lr() and Lc() may be determined via a call to the func-
106 * tion PB_Cnumroc:
107 * Lr( IA, K ) = PB_Cnumroc( K, IA, IMB_A, MB_A, MYROW, RSRC_A, NPROW )
108 * Lc( JA, K ) = PB_Cnumroc( K, JA, INB_A, NB_A, MYCOL, CSRC_A, NPCOL )
109 *
110 * Arguments
111 * =========
112 *
113 * N (global input) INTEGER
114 * On entry, N specifies the length of the subvector sub( X ).
115 * N must be at least zero.
116 *
117 * NORM2 (local output) REAL
118 * On exit, NORM2 specifies the 2-norm of the subvector sub( X )
119 * only in its scope (See below for further details).
120 *
121 * X (local input) COMPLEX array
122 * On entry, X is an array of dimension (LLD_X, Kx), where LLD_X
123 * is at least MAX( 1, Lr( 1, IX ) ) when INCX = M_X and
124 * MAX( 1, Lr( 1, IX+N-1 ) ) otherwise, and, Kx is at least
125 * Lc( 1, JX+N-1 ) when INCX = M_X and Lc( 1, JX ) otherwise.
126 * Before entry, this array contains the local entries of the
127 * matrix X.
128 *
129 * IX (global input) INTEGER
130 * On entry, IX specifies X's global row index, which points to
131 * the beginning of the submatrix sub( X ).
132 *
133 * JX (global input) INTEGER
134 * On entry, JX specifies X's global column index, which points
135 * to the beginning of the submatrix sub( X ).
136 *
137 * DESCX (global and local input) INTEGER array
138 * On entry, DESCX is an integer array of dimension DLEN_. This
139 * is the array descriptor for the matrix X.
140 *
141 * INCX (global input) INTEGER
142 * On entry, INCX specifies the global increment for the
143 * elements of X. Only two values of INCX are supported in
144 * this version, namely 1 and M_X. INCX must not be zero.
145 *
146 * Further Details
147 * ===============
148 *
149 * When the result of a vector-oriented PBLAS call is a scalar, this
150 * scalar is set only within the process scope which owns the vector(s)
151 * being operated on. Let sub( X ) be a generic term for the input vec-
152 * tor(s). Then, the processes owning the correct the answer is determi-
153 * ned as follows: if an operation involves more than one vector, the
154 * processes receiving the result will be the union of the following set
155 * of processes for each vector:
156 *
157 * If N = 1, M_X = 1 and INCX = 1, then one cannot determine if a pro-
158 * cess row or process column owns the vector operand, therefore only
159 * the process owning sub( X ) receives the correct result;
160 *
161 * If INCX = M_X, then sub( X ) is a vector distributed over a process
162 * row. Each process in this row receives the result;
163 *
164 * If INCX = 1, then sub( X ) is a vector distributed over a process
165 * column. Each process in this column receives the result;
166 *
167 * -- Written on April 1, 1998 by
168 * Antoine Petitet, University of Tennessee, Knoxville 37996, USA.
169 *
170 * ---------------------------------------------------------------------
171 */
172 /*
173 * .. Local Scalars ..
174 */
175  char * Xptr = NULL, top;
176  int Xcol, Xi, Xii, Xj, Xjj, Xld, Xnp, Xnq, Xrow, ctxt, dst, dist,
177  info, k, mycol, mydist, myrow, npcol, nprow, src, size;
178  float Xtmp, scale, ssq, temp1, temp2;
179  PBTYP_T * type;
180 /*
181 * .. Local Arrays ..
182 */
183  int Xd[DLEN_];
184  float work[4];
185 /* ..
186 * .. Executable Statements ..
187 *
188 */
189  PB_CargFtoC( *IX, *JX, DESCX, &Xi, &Xj, Xd );
190 #ifndef NO_ARGCHK
191 /*
192 * Test the input parameters
193 */
194  Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
195  if( !( info = ( ( nprow == -1 ) ? -( 601 + CTXT_ ) : 0 ) ) )
196  PB_Cchkvec( ctxt, "PSCNRM2", "X", *N, 1, Xi, Xj, Xd, *INCX, 6, &info );
197  if( info ) { PB_Cabort( ctxt, "PSCNRM2", info ); return; }
198 #endif
199 /*
200 * Initialize NORM2
201 */
202  *NORM2 = ZERO;
203 /*
204 * Quick return if possible
205 */
206  if( *N == 0 ) return;
207 /*
208 * Retrieve process grid information
209 */
210 #ifdef NO_ARGCHK
211  Cblacs_gridinfo( ( ctxt = Xd[CTXT_] ), &nprow, &npcol, &myrow, &mycol );
212 #endif
213 /*
214 * Retrieve sub( X )'s local information: Xii, Xjj, Xrow, Xcol
215 */
216  PB_Cinfog2l( Xi, Xj, Xd, nprow, npcol, myrow, mycol, &Xii, &Xjj,
217  &Xrow, &Xcol );
218 /*
219 * Handle degenerate case separately, sub( X )'s scope is just one process
220 */
221  if( ( *N == 1 ) && ( *INCX == 1 ) && ( Xd[M_] == 1 ) )
222  {
223 /*
224 * Make sure I own some data and compute NORM2
225 */
226  if( ( ( myrow == Xrow ) || ( Xrow < 0 ) ) &&
227  ( ( mycol == Xcol ) || ( Xcol < 0 ) ) )
228  {
229  scale = ZERO;
230  ssq = ONE;
231  type = PB_Cctypeset();
232  Xptr = Mptr( ((char *) X), Xii, Xjj, Xd[LLD_], type->size );
233  Xtmp = ((float *) Xptr)[REAL_PART];
234  if( Xtmp != ZERO )
235  {
236  temp1 = ABS( Xtmp );
237  if( scale < temp1 )
238  {
239  temp2 = scale / temp1;
240  ssq = ONE + ssq * ( temp2 * temp2 );
241  scale = temp1;
242  }
243  else
244  {
245  temp2 = temp1 / scale;
246  ssq = ssq + ( temp2 * temp2 );
247  }
248  }
249  Xtmp = ((float *) Xptr)[IMAG_PART];
250  if( Xtmp != ZERO )
251  {
252  temp1 = ABS( Xtmp );
253  if( scale < temp1 )
254  {
255  temp2 = scale / temp1;
256  ssq = ONE + ssq * ( temp2 * temp2 );
257  scale = temp1;
258  }
259  else
260  {
261  temp2 = temp1 / scale;
262  ssq = ssq + ( temp2 * temp2 );
263  }
264  }
265 /*
266 * Compute NORM2 = SCALE * SQRT( SSQ )
267 */
268  sasqrtb_( &scale, &ssq, NORM2 );
269  }
270  return;
271  }
272  else if( *INCX == Xd[M_] )
273  {
274 /*
275 * sub( X ) resides in (a) process row(s)
276 */
277  if( ( myrow == Xrow ) || ( Xrow < 0 ) )
278  {
279 /*
280 * Initialize SCALE and SSQ
281 */
282  scale = ZERO;
283  ssq = ONE;
284 /*
285 * Make sure I own some data and compute local sum of squares
286 */
287  Xnq = PB_Cnumroc( *N, Xj, Xd[INB_], Xd[NB_], mycol, Xd[CSRC_], npcol );
288  if( Xnq > 0 )
289  {
290  Xld = Xd[LLD_];
291  type = PB_Cctypeset(); size = type->size;
292  Xptr = Mptr( ((char *) X), Xii, Xjj, Xld, size );
293 
294  for( k = 0; k < Xnq; k++ )
295  {
296  Xtmp = ((float *) Xptr)[REAL_PART];
297  if( Xtmp != ZERO )
298  {
299  temp1 = ABS( Xtmp );
300  if( scale < temp1 )
301  {
302  temp2 = scale / temp1;
303  ssq = ONE + ssq * ( temp2 * temp2 );
304  scale = temp1;
305  }
306  else
307  {
308  temp2 = temp1 / scale;
309  ssq = ssq + ( temp2 * temp2 );
310  }
311  }
312  Xtmp = ((float *) Xptr)[IMAG_PART];
313  if( Xtmp != ZERO )
314  {
315  temp1 = ABS( Xtmp );
316  if( scale < temp1 )
317  {
318  temp2 = scale / temp1;
319  ssq = ONE + ssq * ( temp2 * temp2 );
320  scale = temp1;
321  }
322  else
323  {
324  temp2 = temp1 / scale;
325  ssq = ssq + ( temp2 * temp2 );
326  }
327  }
328  Xptr += Xld * size;
329  }
330  }
331 /*
332 * If Xnq <= 0, SCALE is zero and SSQ is one (see initialization above)
333 */
334  if( ( npcol >= 2 ) && ( Xcol >= 0 ) )
335  {
336 /*
337 * Combine the local sum of squares using a 1-tree topology within process row
338 * 0 if npcol > 1 and Xcol >= 0, i.e sub( X ) is distributed.
339 */
340  work[0] = scale;
341  work[1] = ssq;
342 
343  mydist = mycol;
344  k = 1;
345 l_10:
346  if( mydist & 1 )
347  {
348  dist = k * ( mydist - 1 );
349  dst = MPosMod( dist, npcol );
350  Csgesd2d( ctxt, 2, 1, ((char*) work), 2, myrow, dst );
351  goto l_20;
352  }
353  else
354  {
355  dist = mycol + k;
356  src = MPosMod( dist, npcol );
357 
358  if( mycol < src )
359  {
360  Csgerv2d( ctxt, 2, 1, ((char*)&work[2]), 2, myrow, src );
361  if( work[0] >= work[2] )
362  {
363  if( work[0] != ZERO )
364  {
365  temp1 = work[2] / work[0];
366  work[1] = work[1] + ( temp1 * temp1 ) * work[3];
367  }
368  }
369  else
370  {
371  temp1 = work[0] / work[2];
372  work[1] = work[3] + ( temp1 * temp1 ) * work[1];
373  work[0] = work[2];
374  }
375  }
376  mydist >>= 1;
377  }
378  k <<= 1;
379 
380  if( k < npcol ) goto l_10;
381 l_20:
382 /*
383 * Process column 0 broadcasts the combined values of SCALE and SSQ within their
384 * process row.
385 */
386  top = *PB_Ctop( &ctxt, BCAST, ROW, TOP_GET );
387  if( mycol == 0 )
388  {
389  Csgebs2d( ctxt, ROW, &top, 2, 1, ((char*)work), 2 );
390  }
391  else
392  {
393  Csgebr2d( ctxt, ROW, &top, 2, 1, ((char*)work), 2,
394  myrow, 0 );
395  }
396 /*
397 * Compute NORM2 redundantly NORM2 = WORK( 1 ) * SQRT( WORK( 2 ) )
398 */
399  sasqrtb_( &work[0], &work[1], NORM2 );
400  }
401  else
402  {
403 /*
404 * Compute NORM2 redundantly ( sub( X ) is not distributed )
405 */
406  sasqrtb_( &scale, &ssq, NORM2 );
407  }
408  }
409  return;
410  }
411  else
412  {
413 /*
414 * sub( X ) resides in (a) process column(s)
415 */
416  if( ( mycol == Xcol ) || ( Xcol < 0 ) )
417  {
418 /*
419 * Initialize SCALE and SSQ
420 */
421  scale = ZERO;
422  ssq = ONE;
423 /*
424 * Make sure I own some data and compute local sum of squares
425 */
426  Xnp = PB_Cnumroc( *N, Xi, Xd[IMB_], Xd[MB_], myrow, Xd[RSRC_], nprow );
427  if( Xnp > 0 )
428  {
429  type = PB_Cctypeset(); size = type->size;
430  Xptr = Mptr( ((char *) X), Xii, Xjj, Xd[LLD_], size );
431 
432  for( k = 0; k < Xnp; k++ )
433  {
434  Xtmp = ((float *) Xptr)[REAL_PART];
435  if( Xtmp != ZERO )
436  {
437  temp1 = ABS( Xtmp );
438  if( scale < temp1 )
439  {
440  temp2 = scale / temp1;
441  ssq = ONE + ssq * ( temp2 * temp2 );
442  scale = temp1;
443  }
444  else
445  {
446  temp2 = temp1 / scale;
447  ssq = ssq + ( temp2 * temp2 );
448  }
449  }
450  Xtmp = ((float *) Xptr)[IMAG_PART];
451  if( Xtmp != ZERO )
452  {
453  temp1 = ABS( Xtmp );
454  if( scale < temp1 )
455  {
456  temp2 = scale / temp1;
457  ssq = ONE + ssq * ( temp2 * temp2 );
458  scale = temp1;
459  }
460  else
461  {
462  temp2 = temp1 / scale;
463  ssq = ssq + ( temp2 * temp2 );
464  }
465  }
466  Xptr += size;
467  }
468  }
469 /*
470 * If Xnp <= 0, SCALE is zero and SSQ is one (see initialization above)
471 */
472  if( ( nprow >= 2 ) && ( Xrow >= 0 ) )
473  {
474 /*
475 * Combine the local sum of squares using a 1-tree topology within process
476 * column 0 if nprow > 1 and Xrow >= 0, i.e sub( X ) is distributed.
477 */
478  work[0] = scale;
479  work[1] = ssq;
480 
481  mydist = myrow;
482  k = 1;
483 l_30:
484  if( mydist & 1 )
485  {
486  dist = k * ( mydist - 1 );
487  dst = MPosMod( dist, nprow );
488  Csgesd2d( ctxt, 2, 1, ((char*)work), 2, dst, mycol );
489  goto l_40;
490  }
491  else
492  {
493  dist = myrow + k;
494  src = MPosMod( dist, nprow );
495 
496  if( myrow < src )
497  {
498  Csgerv2d( ctxt, 2, 1, ((char*)&work[2]), 2, src, mycol );
499  if( work[0] >= work[2] )
500  {
501  if( work[0] != ZERO )
502  {
503  temp1 = work[2] / work[0];
504  work[1] = work[1] + ( temp1 * temp1 ) * work[3];
505  }
506  }
507  else
508  {
509  temp1 = work[0] / work[2];
510  work[1] = work[3] + ( temp1 * temp1 ) * work[1];
511  work[0] = work[2];
512  }
513  }
514  mydist >>= 1;
515  }
516  k <<= 1;
517 
518  if( k < nprow ) goto l_30;
519 l_40:
520 /*
521 * Process column 0 broadcasts the combined values of SCALE and SSQ within their
522 * process column
523 */
524  top = *PB_Ctop( &ctxt, BCAST, COLUMN, TOP_GET );
525  if( myrow == 0 )
526  {
527  Csgebs2d( ctxt, COLUMN, &top, 2, 1, ((char*)work), 2 );
528  }
529  else
530  {
531  Csgebr2d( ctxt, COLUMN, &top, 2, 1, ((char*)work), 2,
532  0, mycol );
533  }
534 /*
535 * Compute NORM2 redundantly NORM2 = WORK[0] * SQRT( WORK[1] )
536 */
537  sasqrtb_( &work[0], &work[1], NORM2 );
538  }
539  else
540  {
541 /*
542 * Compute NORM2 redundantly ( sub( X ) is not distributed )
543 */
544  sasqrtb_( &scale, &ssq, NORM2 );
545  }
546  }
547  return;
548  }
549 /*
550 * End of PSCNRM2
551 */
552 }
M_
#define M_
Definition: PBtools.h:39
ROW
#define ROW
Definition: PBblacs.h:46
MB_
#define MB_
Definition: PBtools.h:43
NB_
#define NB_
Definition: PBtools.h:44
COLUMN
#define COLUMN
Definition: PBblacs.h:45
CSRC_
#define CSRC_
Definition: PBtools.h:46
PBblacs.h
PBtools.h
Csgesd2d
void Csgesd2d()
Csgebs2d
void Csgebs2d()
PBblas.h
REAL_PART
#define REAL_PART
Definition: pblas.h:135
PBpblas.h
DLEN_
#define DLEN_
Definition: PBtools.h:48
MPosMod
#define MPosMod(I, d)
Definition: PBtools.h:95
LLD_
#define LLD_
Definition: PBtools.h:47
Csgebr2d
void Csgebr2d()
pscnrm2_
void pscnrm2_(int *N, float *NORM2, float *X, int *IX, int *JX, int *DESCX, int *INCX)
Definition: pscnrm2_.c:23
sasqrtb_
F_VOID_FCT sasqrtb_()
ZERO
#define ZERO
Definition: PBtools.h:66
PB_Cchkvec
void PB_Cchkvec()
IMB_
#define IMB_
Definition: PBtools.h:41
PB_Cabort
void PB_Cabort()
TOP_GET
#define TOP_GET
Definition: PBblacs.h:50
Csgerv2d
void Csgerv2d()
PB_Ctop
char * PB_Ctop()
ONE
#define ONE
Definition: PBtools.h:64
RSRC_
#define RSRC_
Definition: PBtools.h:45
PB_CargFtoC
void PB_CargFtoC()
BCAST
#define BCAST
Definition: PBblacs.h:48
PBTYP_T::size
int size
Definition: pblas.h:329
PB_Cinfog2l
void PB_Cinfog2l()
PB_Cnumroc
int PB_Cnumroc()
ABS
#define ABS(a_)
Definition: PBtools.h:75
INB_
#define INB_
Definition: PBtools.h:42
Cblacs_gridinfo
void Cblacs_gridinfo()
PBTYP_T
Definition: pblas.h:325
pblas.h
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_Cctypeset
PBTYP_T * PB_Cctypeset()
Definition: PB_Cctypeset.c:19