ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzlansy.f
Go to the documentation of this file.
1  DOUBLE PRECISION FUNCTION pzlansy( NORM, UPLO, N, A, IA, JA,
2  $ DESCA, WORK )
3  IMPLICIT NONE
4 *
5 * -- ScaLAPACK auxiliary routine (version 1.7) --
6 * University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7 * and University of California, Berkeley.
8 * May 1, 1997
9 *
10 * .. Scalar Arguments ..
11  CHARACTER norm, uplo
12  INTEGER ia, ja, n
13 * ..
14 * .. Array Arguments ..
15  INTEGER desca( * )
16  DOUBLE PRECISION work( * )
17  COMPLEX*16 a( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PZLANSY returns the value of the one norm, or the Frobenius norm,
24 * or the infinity norm, or the element of largest absolute value of a
25 * real symmetric distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
26 *
27 * PZLANSY returns the value
28 *
29 * ( max(abs(A(i,j))), NORM = 'M' or 'm' with IA <= i <= IA+N-1,
30 * ( and JA <= j <= JA+N-1,
31 * (
32 * ( norm1( sub( A ) ), NORM = '1', 'O' or 'o'
33 * (
34 * ( normI( sub( A ) ), NORM = 'I' or 'i'
35 * (
36 * ( normF( sub( A ) ), NORM = 'F', 'f', 'E' or 'e'
37 *
38 * where norm1 denotes the one norm of a matrix (maximum column sum),
39 * normI denotes the infinity norm of a matrix (maximum row sum) and
40 * normF denotes the Frobenius norm of a matrix (square root of sum of
41 * squares). Note that max(abs(A(i,j))) is not a matrix norm.
42 *
43 * Notes
44 * =====
45 *
46 * Each global data object is described by an associated description
47 * vector. This vector stores the information required to establish
48 * the mapping between an object element and its corresponding process
49 * and memory location.
50 *
51 * Let A be a generic term for any 2D block cyclicly distributed array.
52 * Such a global array has an associated description vector DESCA.
53 * In the following comments, the character _ should be read as
54 * "of the global array".
55 *
56 * NOTATION STORED IN EXPLANATION
57 * --------------- -------------- --------------------------------------
58 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
59 * DTYPE_A = 1.
60 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
61 * the BLACS process grid A is distribu-
62 * ted over. The context itself is glo-
63 * bal, but the handle (the integer
64 * value) may vary.
65 * M_A (global) DESCA( M_ ) The number of rows in the global
66 * array A.
67 * N_A (global) DESCA( N_ ) The number of columns in the global
68 * array A.
69 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
70 * the rows of the array.
71 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
72 * the columns of the array.
73 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
74 * row of the array A is distributed.
75 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
76 * first column of the array A is
77 * distributed.
78 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
79 * array. LLD_A >= MAX(1,LOCr(M_A)).
80 *
81 * Let K be the number of rows or columns of a distributed matrix,
82 * and assume that its process grid has dimension p x q.
83 * LOCr( K ) denotes the number of elements of K that a process
84 * would receive if K were distributed over the p processes of its
85 * process column.
86 * Similarly, LOCc( K ) denotes the number of elements of K that a
87 * process would receive if K were distributed over the q processes of
88 * its process row.
89 * The values of LOCr() and LOCc() may be determined via a call to the
90 * ScaLAPACK tool function, NUMROC:
91 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
92 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
93 * An upper bound for these quantities may be computed by:
94 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
95 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
96 *
97 * Arguments
98 * =========
99 *
100 * NORM (global input) CHARACTER
101 * Specifies the value to be returned in PZLANSY as described
102 * above.
103 *
104 * UPLO (global input) CHARACTER
105 * Specifies whether the upper or lower triangular part of the
106 * symmetric matrix sub( A ) is to be referenced.
107 * = 'U': Upper triangular part of sub( A ) is referenced,
108 * = 'L': Lower triangular part of sub( A ) is referenced.
109 *
110 * N (global input) INTEGER
111 * The number of rows and columns to be operated on i.e the
112 * number of rows and columns of the distributed submatrix
113 * sub( A ). When N = 0, PZLANSY is set to zero. N >= 0.
114 *
115 * A (local input) COMPLEX*16 pointer into the local memory
116 * to an array of dimension (LLD_A, LOCc(JA+N-1)) containing the
117 * local pieces of the symmetric distributed matrix sub( A ).
118 * If UPLO = 'U', the leading N-by-N upper triangular part of
119 * sub( A ) contains the upper triangular matrix which norm is
120 * to be computed, and the strictly lower triangular part of
121 * this matrix is not referenced. If UPLO = 'L', the leading
122 * N-by-N lower triangular part of sub( A ) contains the lower
123 * triangular matrix which norm is to be computed, and the
124 * strictly upper triangular part of sub( A ) is not referenced.
125 *
126 * IA (global input) INTEGER
127 * The row index in the global array A indicating the first
128 * row of sub( A ).
129 *
130 * JA (global input) INTEGER
131 * The column index in the global array A indicating the
132 * first column of sub( A ).
133 *
134 * DESCA (global and local input) INTEGER array of dimension DLEN_.
135 * The array descriptor for the distributed matrix A.
136 *
137 * WORK (local workspace) DOUBLE PRECISION array dimension (LWORK)
138 * LWORK >= 0 if NORM = 'M' or 'm' (not referenced),
139 * 2*Nq0+Np0+LDW if NORM = '1', 'O', 'o', 'I' or 'i',
140 * where LDW is given by:
141 * IF( NPROW.NE.NPCOL ) THEN
142 * LDW = MB_A*CEIL(CEIL(Np0/MB_A)/(LCM/NPROW))
143 * ELSE
144 * LDW = 0
145 * END IF
146 * 0 if NORM = 'F', 'f', 'E' or 'e' (not referenced),
147 *
148 * where LCM is the least common multiple of NPROW and NPCOL
149 * LCM = ILCM( NPROW, NPCOL ) and CEIL denotes the ceiling
150 * operation (ICEIL).
151 *
152 * IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
153 * IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
154 * IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
155 * Np0 = NUMROC( N+IROFFA, MB_A, MYROW, IAROW, NPROW ),
156 * Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
157 *
158 * ICEIL, ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
159 * MYROW, MYCOL, NPROW and NPCOL can be determined by calling
160 * the subroutine BLACS_GRIDINFO.
161 *
162 * =====================================================================
163 *
164 * .. Parameters ..
165  INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
166  $ lld_, mb_, m_, nb_, n_, rsrc_
167  parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
168  $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
169  $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
170  DOUBLE PRECISION one, zero
171  parameter( one = 1.0d+0, zero = 0.0d+0 )
172 * ..
173 * .. Local Scalars ..
174  INTEGER i, iarow, iacol, ib, icoff, ictxt, icurcol,
175  $ icurrow, ii, iia, in, iroff, icsr, icsr0,
176  $ ioffa, irsc, irsc0, irsr, irsr0, jj, jja, k,
177  $ lda, ll, mycol, myrow, np, npcol, nprow, nq
178  DOUBLE PRECISION sum, value
179 * ..
180 * .. Local Arrays ..
181  DOUBLE PRECISION ssq( 2 ), colssq( 2 )
182 * ..
183 * .. External Subroutines ..
184  EXTERNAL blacs_gridinfo, daxpy, dcombssq,
185  $ dgamx2d, dgsum2d, dgebr2d,
186  $ dgebs2d, pdcol2row, pdtreecomb,
187  $ zlassq
188 * ..
189 * .. External Functions ..
190  LOGICAL lsame
191  INTEGER iceil, idamax, numroc
192  EXTERNAL iceil, idamax, lsame, numroc
193 * ..
194 * .. Intrinsic Functions ..
195  INTRINSIC abs, max, min, mod, sqrt
196 * ..
197 * .. Executable Statements ..
198 *
199 * Get grid parameters and local indexes.
200 *
201  ictxt = desca( ctxt_ )
202  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
203  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol,
204  $ iia, jja, iarow, iacol )
205 *
206  iroff = mod( ia-1, desca( mb_ ) )
207  icoff = mod( ja-1, desca( nb_ ) )
208  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
209  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
210  icsr = 1
211  irsr = icsr + nq
212  irsc = irsr + nq
213  IF( myrow.EQ.iarow ) THEN
214  irsc0 = irsc + iroff
215  np = np - iroff
216  ELSE
217  irsc0 = irsc
218  END IF
219  IF( mycol.EQ.iacol ) THEN
220  icsr0 = icsr + icoff
221  irsr0 = irsr + icoff
222  nq = nq - icoff
223  ELSE
224  icsr0 = icsr
225  irsr0 = irsr
226  END IF
227  in = min( iceil( ia, desca( mb_ ) ) * desca( mb_ ), ia+n-1 )
228  lda = desca( lld_ )
229 *
230 * If the matrix is symmetric, we address only a triangular portion
231 * of the matrix. A sum of row (column) i of the complete matrix
232 * can be obtained by adding along row i and column i of the the
233 * triangular matrix, stopping/starting at the diagonal, which is
234 * the point of reflection. The pictures below demonstrate this.
235 * In the following code, the row sums created by --- rows below are
236 * refered to as ROWSUMS, and the column sums shown by | are refered
237 * to as COLSUMS. Infinity-norm = 1-norm = ROWSUMS+COLSUMS.
238 *
239 * UPLO = 'U' UPLO = 'L'
240 * ____i______ ___________
241 * |\ | | |\ |
242 * | \ | | | \ |
243 * | \ | | | \ |
244 * | \|------| i i|---\ |
245 * | \ | | |\ |
246 * | \ | | | \ |
247 * | \ | | | \ |
248 * | \ | | | \ |
249 * | \ | | | \ |
250 * | \ | | | \ |
251 * |__________\| |___|______\|
252 * i
253 *
254 * II, JJ : local indices into array A
255 * ICURROW : process row containing diagonal block
256 * ICURCOL : process column containing diagonal block
257 * IRSC0 : pointer to part of work used to store the ROWSUMS while
258 * they are stored along a process column
259 * IRSR0 : pointer to part of work used to store the ROWSUMS after
260 * they have been transposed to be along a process row
261 *
262  ii = iia
263  jj = jja
264 *
265  IF( n.EQ.0 ) THEN
266 *
267  VALUE = zero
268 *
269 ************************************************************************
270 * max norm
271 *
272  ELSE IF( lsame( norm, 'M' ) ) THEN
273 *
274 * Find max(abs(A(i,j))).
275 *
276  VALUE = zero
277 *
278  IF( lsame( uplo, 'U' ) ) THEN
279 *
280 * Handle first block separately
281 *
282  ib = in-ia+1
283 *
284 * Find COLMAXS
285 *
286  IF( mycol.EQ.iacol ) THEN
287  DO 20 k = (jj-1)*lda, (jj+ib-2)*lda, lda
288  IF( ii.GT.iia ) THEN
289  DO 10 ll = iia, ii-1
290  VALUE = max( VALUE, abs( a( ll+k ) ) )
291  10 CONTINUE
292  END IF
293  IF( myrow.EQ.iarow )
294  $ ii = ii + 1
295  20 CONTINUE
296 *
297 * Reset local indices so we can find ROWMAXS
298 *
299  IF( myrow.EQ.iarow )
300  $ ii = ii - ib
301 *
302  END IF
303 *
304 * Find ROWMAXS
305 *
306  IF( myrow.EQ.iarow ) THEN
307  DO 40 k = ii, ii+ib-1
308  IF( jj.LE.jja+nq-1 ) THEN
309  DO 30 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
310  VALUE = max( VALUE, abs( a( k+ll ) ) )
311  30 CONTINUE
312  END IF
313  IF( mycol.EQ.iacol )
314  $ jj = jj + 1
315  40 CONTINUE
316  ii = ii + ib
317  ELSE IF( mycol.EQ.iacol ) THEN
318  jj = jj + ib
319  END IF
320 *
321  icurrow = mod( iarow+1, nprow )
322  icurcol = mod( iacol+1, npcol )
323 *
324 * Loop over the remaining rows/columns of the matrix.
325 *
326  DO 90 i = in+1, ia+n-1, desca( mb_ )
327  ib = min( desca( mb_ ), ia+n-i )
328 *
329 * Find COLMAXS
330 *
331  IF( mycol.EQ.icurcol ) THEN
332  DO 60 k = (jj-1)*lda, (jj+ib-2)*lda, lda
333  IF( ii.GT.iia ) THEN
334  DO 50 ll = iia, ii-1
335  VALUE = max( VALUE, abs( a( ll+k ) ) )
336  50 CONTINUE
337  END IF
338  IF( myrow.EQ.icurrow )
339  $ ii = ii + 1
340  60 CONTINUE
341 *
342 * Reset local indices so we can find ROWMAXS
343 *
344  IF( myrow.EQ.icurrow )
345  $ ii = ii - ib
346  END IF
347 *
348 * Find ROWMAXS
349 *
350  IF( myrow.EQ.icurrow ) THEN
351  DO 80 k = ii, ii+ib-1
352  IF( jj.LE.jja+nq-1 ) THEN
353  DO 70 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
354  VALUE = max( VALUE, abs( a( k+ll ) ) )
355  70 CONTINUE
356  END IF
357  IF( mycol.EQ.icurcol )
358  $ jj = jj + 1
359  80 CONTINUE
360  ii = ii + ib
361  ELSE IF( mycol.EQ.icurcol ) THEN
362  jj = jj + ib
363  END IF
364  icurrow = mod( icurrow+1, nprow )
365  icurcol = mod( icurcol+1, npcol )
366  90 CONTINUE
367 *
368  ELSE
369 *
370 * Handle first block separately
371 *
372  ib = in-ia+1
373 *
374 * Find COLMAXS
375 *
376  IF( mycol.EQ.iacol ) THEN
377  DO 110 k = (jj-1)*lda, (jj+ib-2)*lda, lda
378  IF( ii.LE.iia+np-1 ) THEN
379  DO 100 ll = ii, iia+np-1
380  VALUE = max( VALUE, abs( a( ll+k ) ) )
381  100 CONTINUE
382  END IF
383  IF( myrow.EQ.iarow )
384  $ ii = ii + 1
385  110 CONTINUE
386 *
387 * Reset local indices so we can find ROWMAXS
388 *
389  IF( myrow.EQ.iarow )
390  $ ii = ii - ib
391  END IF
392 *
393 * Find ROWMAXS
394 *
395  IF( myrow.EQ.iarow ) THEN
396  DO 130 k = 0, ib-1
397  IF( jj.GT.jja ) THEN
398  DO 120 ll = (jja-1)*lda, (jj-2)*lda, lda
399  VALUE = max( VALUE, abs( a( ii+ll ) ) )
400  120 CONTINUE
401  END IF
402  ii = ii + 1
403  IF( mycol.EQ.iacol )
404  $ jj = jj + 1
405  130 CONTINUE
406  ELSE IF( mycol.EQ.iacol ) THEN
407  jj = jj + ib
408  END IF
409 *
410  icurrow = mod( iarow+1, nprow )
411  icurcol = mod( iacol+1, npcol )
412 *
413 * Loop over rows/columns of global matrix.
414 *
415  DO 180 i = in+1, ia+n-1, desca( mb_ )
416  ib = min( desca( mb_ ), ia+n-i )
417 *
418 * Find COLMAXS
419 *
420  IF( mycol.EQ.icurcol ) THEN
421  DO 150 k = (jj-1)*lda, (jj+ib-2)*lda, lda
422  IF( ii.LE.iia+np-1 ) THEN
423  DO 140 ll = ii, iia+np-1
424  VALUE = max( VALUE, abs( a( ll+k ) ) )
425  140 CONTINUE
426  END IF
427  IF( myrow.EQ.icurrow )
428  $ ii = ii + 1
429  150 CONTINUE
430 *
431 * Reset local indices so we can find ROWMAXS
432 *
433  IF( myrow.EQ.icurrow )
434  $ ii = ii - ib
435  END IF
436 *
437 * Find ROWMAXS
438 *
439  IF( myrow.EQ.icurrow ) THEN
440  DO 170 k = 0, ib-1
441  IF( jj.GT.jja ) THEN
442  DO 160 ll = (jja-1)*lda, (jj-2)*lda, lda
443  VALUE = max( VALUE, abs( a( ii+ll ) ) )
444  160 CONTINUE
445  END IF
446  ii = ii + 1
447  IF( mycol.EQ.icurcol )
448  $ jj = jj + 1
449  170 CONTINUE
450  ELSE IF( mycol.EQ.icurcol ) THEN
451  jj = jj + ib
452  END IF
453  icurrow = mod( icurrow+1, nprow )
454  icurcol = mod( icurcol+1, npcol )
455 *
456  180 CONTINUE
457 *
458  END IF
459 *
460 * Gather the result on process (IAROW,IACOL).
461 *
462  CALL dgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, i, k, -1,
463  $ iarow, iacol )
464 *
465 ************************************************************************
466 * one or inf norm
467 *
468  ELSE IF( lsame( norm, 'I' ) .OR. lsame( norm, 'O' ) .OR.
469  $ norm.EQ.'1' ) THEN
470 *
471 * Find normI( sub( A ) ) ( = norm1( sub( A ) ), since sub( A ) is
472 * symmetric).
473 *
474  IF( lsame( uplo, 'U' ) ) THEN
475 *
476 * Handle first block separately
477 *
478  ib = in-ia+1
479 *
480 * Find COLSUMS
481 *
482  IF( mycol.EQ.iacol ) THEN
483  ioffa = ( jj - 1 ) * lda
484  DO 200 k = 0, ib-1
485  sum = zero
486  IF( ii.GT.iia ) THEN
487  DO 190 ll = iia, ii-1
488  sum = sum + abs( a( ll+ioffa ) )
489  190 CONTINUE
490  END IF
491  ioffa = ioffa + lda
492  work( jj+k-jja+icsr0 ) = sum
493  IF( myrow.EQ.iarow )
494  $ ii = ii + 1
495  200 CONTINUE
496 *
497 * Reset local indices so we can find ROWSUMS
498 *
499  IF( myrow.EQ.iarow )
500  $ ii = ii - ib
501 *
502  END IF
503 *
504 * Find ROWSUMS
505 *
506  IF( myrow.EQ.iarow ) THEN
507  DO 220 k = ii, ii+ib-1
508  sum = zero
509  IF( jja+nq.GT.jj ) THEN
510  DO 210 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
511  sum = sum + abs( a( k+ll ) )
512  210 CONTINUE
513  END IF
514  work( k-iia+irsc0 ) = sum
515  IF( mycol.EQ.iacol )
516  $ jj = jj + 1
517  220 CONTINUE
518  ii = ii + ib
519  ELSE IF( mycol.EQ.iacol ) THEN
520  jj = jj + ib
521  END IF
522 *
523  icurrow = mod( iarow+1, nprow )
524  icurcol = mod( iacol+1, npcol )
525 *
526 * Loop over remaining rows/columns of global matrix.
527 *
528  DO 270 i = in+1, ia+n-1, desca( mb_ )
529  ib = min( desca( mb_ ), ia+n-i )
530 *
531 * Find COLSUMS
532 *
533  IF( mycol.EQ.icurcol ) THEN
534  ioffa = ( jj - 1 ) * lda
535  DO 240 k = 0, ib-1
536  sum = zero
537  IF( ii.GT.iia ) THEN
538  DO 230 ll = iia, ii-1
539  sum = sum + abs( a( ioffa+ll ) )
540  230 CONTINUE
541  END IF
542  ioffa = ioffa + lda
543  work( jj+k-jja+icsr0 ) = sum
544  IF( myrow.EQ.icurrow )
545  $ ii = ii + 1
546  240 CONTINUE
547 *
548 * Reset local indices so we can find ROWSUMS
549 *
550  IF( myrow.EQ.icurrow )
551  $ ii = ii - ib
552 *
553  END IF
554 *
555 * Find ROWSUMS
556 *
557  IF( myrow.EQ.icurrow ) THEN
558  DO 260 k = ii, ii+ib-1
559  sum = zero
560  IF( jja+nq.GT.jj ) THEN
561  DO 250 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
562  sum = sum + abs( a( k+ll ) )
563  250 CONTINUE
564  END IF
565  work( k-iia+irsc0 ) = sum
566  IF( mycol.EQ.icurcol )
567  $ jj = jj + 1
568  260 CONTINUE
569  ii = ii + ib
570  ELSE IF( mycol.EQ.icurcol ) THEN
571  jj = jj + ib
572  END IF
573 *
574  icurrow = mod( icurrow+1, nprow )
575  icurcol = mod( icurcol+1, npcol )
576 *
577  270 CONTINUE
578 *
579  ELSE
580 *
581 * Handle first block separately
582 *
583  ib = in-ia+1
584 *
585 * Find COLSUMS
586 *
587  IF( mycol.EQ.iacol ) THEN
588  ioffa = (jj-1)*lda
589  DO 290 k = 0, ib-1
590  sum = zero
591  IF( iia+np.GT.ii ) THEN
592  DO 280 ll = ii, iia+np-1
593  sum = sum + abs( a( ioffa+ll ) )
594  280 CONTINUE
595  END IF
596  ioffa = ioffa + lda
597  work( jj+k-jja+icsr0 ) = sum
598  IF( myrow.EQ.iarow )
599  $ ii = ii + 1
600  290 CONTINUE
601 *
602 * Reset local indices so we can find ROWSUMS
603 *
604  IF( myrow.EQ.iarow )
605  $ ii = ii - ib
606 *
607  END IF
608 *
609 * Find ROWSUMS
610 *
611  IF( myrow.EQ.iarow ) THEN
612  DO 310 k = ii, ii+ib-1
613  sum = zero
614  IF( jj.GT.jja ) THEN
615  DO 300 ll = (jja-1)*lda, (jj-2)*lda, lda
616  sum = sum + abs( a( k+ll ) )
617  300 CONTINUE
618  END IF
619  work( k-iia+irsc0 ) = sum
620  IF( mycol.EQ.iacol )
621  $ jj = jj + 1
622  310 CONTINUE
623  ii = ii + ib
624  ELSE IF( mycol.EQ.iacol ) THEN
625  jj = jj + ib
626  END IF
627 *
628  icurrow = mod( iarow+1, nprow )
629  icurcol = mod( iacol+1, npcol )
630 *
631 * Loop over rows/columns of global matrix.
632 *
633  DO 360 i = in+1, ia+n-1, desca( mb_ )
634  ib = min( desca( mb_ ), ia+n-i )
635 *
636 * Find COLSUMS
637 *
638  IF( mycol.EQ.icurcol ) THEN
639  ioffa = ( jj - 1 ) * lda
640  DO 330 k = 0, ib-1
641  sum = zero
642  IF( iia+np.GT.ii ) THEN
643  DO 320 ll = ii, iia+np-1
644  sum = sum + abs( a( ll+ioffa ) )
645  320 CONTINUE
646  END IF
647  ioffa = ioffa + lda
648  work( jj+k-jja+icsr0 ) = sum
649  IF( myrow.EQ.icurrow )
650  $ ii = ii + 1
651  330 CONTINUE
652 *
653 * Reset local indices so we can find ROWSUMS
654 *
655  IF( myrow.EQ.icurrow )
656  $ ii = ii - ib
657 *
658  END IF
659 *
660 * Find ROWSUMS
661 *
662  IF( myrow.EQ.icurrow ) THEN
663  DO 350 k = ii, ii+ib-1
664  sum = zero
665  IF( jj.GT.jja ) THEN
666  DO 340 ll = (jja-1)*lda, (jj-2)*lda, lda
667  sum = sum + abs( a( k+ll ) )
668  340 CONTINUE
669  END IF
670  work(k-iia+irsc0) = sum
671  IF( mycol.EQ.icurcol )
672  $ jj = jj + 1
673  350 CONTINUE
674  ii = ii + ib
675  ELSE IF( mycol.EQ.icurcol ) THEN
676  jj = jj + ib
677  END IF
678 *
679  icurrow = mod( icurrow+1, nprow )
680  icurcol = mod( icurcol+1, npcol )
681 *
682  360 CONTINUE
683  END IF
684 *
685 * After calls to DGSUM2D, process row 0 will have global
686 * COLSUMS and process column 0 will have global ROWSUMS.
687 * Transpose ROWSUMS and add to COLSUMS to get global row/column
688 * sum, the max of which is the infinity or 1 norm.
689 *
690  IF( mycol.EQ.iacol )
691  $ nq = nq + icoff
692  CALL dgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work( icsr ), 1,
693  $ iarow, mycol )
694  IF( myrow.EQ.iarow )
695  $ np = np + iroff
696  CALL dgsum2d( ictxt, 'Rowwise', ' ', np, 1, work( irsc ),
697  $ max( 1, np ), myrow, iacol )
698 *
699  CALL pdcol2row( ictxt, n, 1, desca( mb_ ), work( irsc ),
700  $ max( 1, np ), work( irsr ), max( 1, nq ),
701  $ iarow, iacol, iarow, iacol, work( irsc+np ) )
702 *
703  IF( myrow.EQ.iarow ) THEN
704  IF( mycol.EQ.iacol )
705  $ nq = nq - icoff
706  CALL daxpy( nq, one, work( irsr0 ), 1, work( icsr0 ), 1 )
707  IF( nq.LT.1 ) THEN
708  VALUE = zero
709  ELSE
710  VALUE = work( idamax( nq, work( icsr0 ), 1 ) )
711  END IF
712  CALL dgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, i, k,
713  $ -1, iarow, iacol )
714  END IF
715 *
716 ************************************************************************
717 * Frobenius norm
718 * SSQ(1) is scale
719 * SSQ(2) is sum-of-squares
720 *
721  ELSE IF( lsame( norm, 'F' ) .OR. lsame( norm, 'E' ) ) THEN
722 *
723 * Find normF( sub( A ) ).
724 *
725  ssq(1) = zero
726  ssq(2) = one
727 *
728 * Add off-diagonal entries, first
729 *
730  IF( lsame( uplo, 'U' ) ) THEN
731 *
732 * Handle first block separately
733 *
734  ib = in-ia+1
735 *
736  IF( mycol.EQ.iacol ) THEN
737  DO 370 k = (jj-1)*lda, (jj+ib-2)*lda, lda
738  colssq(1) = zero
739  colssq(2) = one
740  CALL zlassq( ii-iia, a( iia+k ), 1,
741  $ colssq(1), colssq(2) )
742  IF( myrow.EQ.iarow )
743  $ ii = ii + 1
744  CALL zlassq( ii-iia, a( iia+k ), 1,
745  $ colssq(1), colssq(2) )
746  CALL dcombssq( ssq, colssq )
747  370 CONTINUE
748 *
749  jj = jj + ib
750  ELSE IF( myrow.EQ.iarow ) THEN
751  ii = ii + ib
752  END IF
753 *
754  icurrow = mod( iarow+1, nprow )
755  icurcol = mod( iacol+1, npcol )
756 *
757 * Loop over rows/columns of global matrix.
758 *
759  DO 390 i = in+1, ia+n-1, desca( mb_ )
760  ib = min( desca( mb_ ), ia+n-i )
761 *
762  IF( mycol.EQ.icurcol ) THEN
763  DO 380 k = (jj-1)*lda, (jj+ib-2)*lda, lda
764  colssq(1) = zero
765  colssq(2) = one
766  CALL zlassq( ii-iia, a( iia+k ), 1,
767  $ colssq(1), colssq(2) )
768  IF( myrow.EQ.icurrow )
769  $ ii = ii + 1
770  CALL zlassq( ii-iia, a(iia+k ), 1,
771  $ colssq(1), colssq(2) )
772  CALL dcombssq( ssq, colssq )
773  380 CONTINUE
774 *
775  jj = jj + ib
776  ELSE IF( myrow.EQ.icurrow ) THEN
777  ii = ii + ib
778  END IF
779 *
780  icurrow = mod( icurrow+1, nprow )
781  icurcol = mod( icurcol+1, npcol )
782 *
783  390 CONTINUE
784 *
785  ELSE
786 *
787 * Handle first block separately
788 *
789  ib = in-ia+1
790 *
791  IF( mycol.EQ.iacol ) THEN
792  DO 400 k = (jj-1)*lda, (jj+ib-2)*lda, lda
793  colssq(1) = zero
794  colssq(2) = one
795  CALL zlassq( iia+np-ii, a( ii+k ), 1,
796  $ colssq(1), colssq(2) )
797  IF( myrow.EQ.iarow )
798  $ ii = ii + 1
799  CALL zlassq( iia+np-ii, a( ii+k ), 1,
800  $ colssq(1), colssq(2) )
801  CALL dcombssq( ssq, colssq )
802  400 CONTINUE
803 *
804  jj = jj + ib
805  ELSE IF( myrow.EQ.iarow ) THEN
806  ii = ii + ib
807  END IF
808 *
809  icurrow = mod( iarow+1, nprow )
810  icurcol = mod( iacol+1, npcol )
811 *
812 * Loop over rows/columns of global matrix.
813 *
814  DO 420 i = in+1, ia+n-1, desca( mb_ )
815  ib = min( desca( mb_ ), ia+n-i )
816 *
817  IF( mycol.EQ.icurcol ) THEN
818  DO 410 k = (jj-1)*lda, (jj+ib-2)*lda, lda
819  colssq(1) = zero
820  colssq(2) = one
821  CALL zlassq( iia+np-ii, a( ii+k ), 1,
822  $ colssq(1), colssq(2) )
823  IF( myrow.EQ.icurrow )
824  $ ii = ii + 1
825  CALL zlassq( iia+np-ii, a( ii+k ), 1,
826  $ colssq(1), colssq(2) )
827  CALL dcombssq( ssq, colssq )
828  410 CONTINUE
829 *
830  jj = jj + ib
831  ELSE IF( myrow.EQ.icurrow ) THEN
832  ii = ii + ib
833  END IF
834 *
835  icurrow = mod( icurrow+1, nprow )
836  icurcol = mod( icurcol+1, npcol )
837 *
838  420 CONTINUE
839 *
840  END IF
841 *
842 * Perform the global scaled sum
843 *
844  CALL pdtreecomb( ictxt, 'All', 2, ssq, iarow, iacol,
845  $ dcombssq )
846  VALUE = ssq( 1 ) * sqrt( ssq( 2 ) )
847 *
848  END IF
849 *
850 * Broadcast the result to the other processes
851 *
852  IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) THEN
853  CALL dgebs2d( ictxt, 'All', ' ', 1, 1, VALUE, 1 )
854  ELSE
855  CALL dgebr2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, iarow,
856  $ iacol )
857  END IF
858 *
859  pzlansy = VALUE
860 *
861  RETURN
862 *
863 * End of PZLANSY
864 *
865  END
max
#define max(A, B)
Definition: pcgemr.c:180
pdtreecomb
subroutine pdtreecomb(ICTXT, SCOPE, N, MINE, RDEST0, CDEST0, SUBPTR)
Definition: pdtreecomb.f:3
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
dcombssq
subroutine dcombssq(V1, V2)
Definition: pdtreecomb.f:259
pzlansy
double precision function pzlansy(NORM, UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pzlansy.f:3
pdcol2row
subroutine pdcol2row(ICTXT, M, N, NB, VS, LDVS, VD, LDVD, RSRC, CSRC, RDEST, CDEST, WORK)
Definition: pdcol2row.f:3
min
#define min(A, B)
Definition: pcgemr.c:181
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2