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