ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pclanhe.f
Go to the documentation of this file.
1  REAL FUNCTION PCLANHE( 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  REAL work( * )
16  COMPLEX a( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PCLANHE 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 * PCLANHE 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 PCLANHE 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, PCLANHE is set to zero. N >= 0.
113 *
114 * A (local input) COMPLEX 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) REAL 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  REAL one, zero
170  parameter( one = 1.0e+0, zero = 0.0e+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  REAL absa, scale, sum, value
178 * ..
179 * .. Local Arrays ..
180  REAL rwork( 2 )
181 * ..
182 * .. External Subroutines ..
183  EXTERNAL blacs_gridinfo, classq, pscol2row,
184  $ pstreecomb, saxpy, scombssq,
185  $ sgamx2d, sgsum2d, sgebr2d, sgebs2d
186 * ..
187 * .. External Functions ..
188  LOGICAL lsame
189  INTEGER iceil, isamax, numroc
190  EXTERNAL iceil, isamax, lsame, numroc
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC abs, max, min, mod, real, sqrt
194 * ..
195 * .. Executable Statements ..
196 *
197 * Get grid parameters and local indexes.
198 *
199  ictxt = desca( ctxt_ )
200  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
201  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol,
202  $ iia, jja, iarow, iacol )
203 *
204  iroff = mod( ia-1, desca( mb_ ) )
205  icoff = mod( ja-1, desca( nb_ ) )
206  np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
207  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
208  icsr = 1
209  irsr = icsr + nq
210  irsc = irsr + nq
211  IF( myrow.EQ.iarow ) THEN
212  irsc0 = irsc + iroff
213  np = np - iroff
214  ELSE
215  irsc0 = irsc
216  END IF
217  IF( mycol.EQ.iacol ) THEN
218  icsr0 = icsr + icoff
219  irsr0 = irsr + icoff
220  nq = nq - icoff
221  ELSE
222  icsr0 = icsr
223  irsr0 = irsr
224  END IF
225  in = min( iceil( ia, desca( mb_ ) ) * desca( mb_ ), ia+n-1 )
226  lda = desca( lld_ )
227 *
228 * If the matrix is Hermitian, we address only a triangular portion
229 * of the matrix. A sum of row (column) i of the complete matrix
230 * can be obtained by adding along row i and column i of the the
231 * triangular matrix, stopping/starting at the diagonal, which is
232 * the point of reflection. The pictures below demonstrate this.
233 * In the following code, the row sums created by --- rows below are
234 * refered to as ROWSUMS, and the column sums shown by | are refered
235 * to as COLSUMS. Infinity-norm = 1-norm = ROWSUMS+COLSUMS.
236 *
237 * UPLO = 'U' UPLO = 'L'
238 * ____i______ ___________
239 * |\ | | |\ |
240 * | \ | | | \ |
241 * | \ | | | \ |
242 * | \|------| i i|---\ |
243 * | \ | | |\ |
244 * | \ | | | \ |
245 * | \ | | | \ |
246 * | \ | | | \ |
247 * | \ | | | \ |
248 * | \ | | | \ |
249 * |__________\| |___|______\|
250 * i
251 *
252 * II, JJ : local indices into array A
253 * ICURROW : process row containing diagonal block
254 * ICURCOL : process column containing diagonal block
255 * IRSC0 : pointer to part of work used to store the ROWSUMS while
256 * they are stored along a process column
257 * IRSR0 : pointer to part of work used to store the ROWSUMS after
258 * they have been transposed to be along a process row
259 *
260  ii = iia
261  jj = jja
262 *
263  IF( n.EQ.0 ) THEN
264 *
265  VALUE = zero
266 *
267  ELSE IF( lsame( norm, 'M' ) ) THEN
268 *
269 * Find max(abs(A(i,j))).
270 *
271  VALUE = zero
272 *
273  IF( lsame( uplo, 'U' ) ) THEN
274 *
275 * Handle first block separately
276 *
277  ib = in-ia+1
278 *
279 * Find COLMAXS
280 *
281  IF( mycol.EQ.iacol ) THEN
282  DO 20 k = (jj-1)*lda, (jj+ib-2)*lda, lda
283  IF( ii.GT.iia ) THEN
284  DO 10 ll = iia, ii-1
285  VALUE = max( VALUE, abs( a( ll+k ) ) )
286  10 CONTINUE
287  END IF
288  IF( myrow.EQ.iarow )
289  $ ii = ii + 1
290  20 CONTINUE
291 *
292 * Reset local indices so we can find ROWMAXS
293 *
294  IF( myrow.EQ.iarow )
295  $ ii = ii - ib
296 *
297  END IF
298 *
299 * Find ROWMAXS
300 *
301  IF( myrow.EQ.iarow ) THEN
302  DO 40 k = ii, ii+ib-1
303  IF( mycol.EQ.iacol ) THEN
304  IF( jj.LE.jja+nq-1 ) THEN
305  VALUE = max( VALUE,
306  $ abs( real( a( k+(jj-1)*lda ) ) ) )
307  DO 30 ll = jj*lda, (jja+nq-2)*lda, lda
308  VALUE = max( VALUE, abs( a( k+ll ) ) )
309  30 CONTINUE
310  END IF
311  ELSE
312  IF( jj.LE.jja+nq-1 ) THEN
313  DO 35 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
314  VALUE = max( VALUE, abs( a( k+ll ) ) )
315  35 CONTINUE
316  END IF
317  END IF
318  IF( mycol.EQ.iacol )
319  $ jj = jj + 1
320  40 CONTINUE
321  ii = ii + ib
322  ELSE IF( mycol.EQ.iacol ) THEN
323  jj = jj + ib
324  END IF
325 *
326  icurrow = mod( iarow+1, nprow )
327  icurcol = mod( iacol+1, npcol )
328 *
329 * Loop over the remaining rows/columns of the matrix.
330 *
331  DO 90 i = in+1, ia+n-1, desca( mb_ )
332  ib = min( desca( mb_ ), ia+n-i )
333 *
334 * Find COLMAXS
335 *
336  IF( mycol.EQ.icurcol ) THEN
337  DO 60 k = (jj-1)*lda, (jj+ib-2)*lda, lda
338  IF( ii.GT.iia ) THEN
339  DO 50 ll = iia, ii-1
340  VALUE = max( VALUE, abs( a( ll+k ) ) )
341  50 CONTINUE
342  END IF
343  IF( myrow.EQ.icurrow )
344  $ ii = ii + 1
345  60 CONTINUE
346 *
347 * Reset local indices so we can find ROWMAXS
348 *
349  IF( myrow.EQ.icurrow )
350  $ ii = ii - ib
351  END IF
352 *
353 * Find ROWMAXS
354 *
355  IF( myrow.EQ.icurrow ) THEN
356  DO 80 k = ii, ii+ib-1
357  IF( mycol.EQ.icurcol ) THEN
358  IF( jj.LE.jja+nq-1 ) THEN
359  VALUE = max( VALUE,
360  $ abs( real( a( k+(jj-1)*lda ) ) ) )
361  DO 70 ll = jj*lda, (jja+nq-2)*lda, lda
362  VALUE = max( VALUE, abs( a( k+ll ) ) )
363  70 CONTINUE
364  END IF
365  ELSE
366  IF( jj.LE.jja+nq-1 ) THEN
367  DO 75 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
368  VALUE = max( VALUE, abs( a( k+ll ) ) )
369  75 CONTINUE
370  END IF
371  END IF
372  IF( mycol.EQ.icurcol )
373  $ jj = jj + 1
374  80 CONTINUE
375  ii = ii + ib
376  ELSE IF( mycol.EQ.icurcol ) THEN
377  jj = jj + ib
378  END IF
379  icurrow = mod( icurrow+1, nprow )
380  icurcol = mod( icurcol+1, npcol )
381  90 CONTINUE
382 *
383  ELSE
384 *
385 * Handle first block separately
386 *
387  ib = in-ia+1
388 *
389 * Find COLMAXS
390 *
391  IF( mycol.EQ.iacol ) THEN
392  DO 110 k = (jj-1)*lda, (jj+ib-2)*lda, lda
393  IF( myrow.EQ.iarow ) THEN
394  IF( ii.LE.iia+np-1 ) THEN
395  VALUE = max( VALUE, abs( real( a( ii+k ) ) ) )
396  DO 100 ll = ii+1, iia+np-1
397  VALUE = max( VALUE, abs( a( ll+k ) ) )
398  100 CONTINUE
399  END IF
400  ELSE
401  IF( ii.LE.iia+np-1 ) THEN
402  DO 105 ll = ii, iia+np-1
403  VALUE = max( VALUE, abs( a( ll+k ) ) )
404  105 CONTINUE
405  END IF
406  END IF
407  IF( myrow.EQ.iarow )
408  $ ii = ii + 1
409  110 CONTINUE
410 *
411 * Reset local indices so we can find ROWMAXS
412 *
413  IF( myrow.EQ.iarow )
414  $ ii = ii - ib
415  END IF
416 *
417 * Find ROWMAXS
418 *
419  IF( myrow.EQ.iarow ) THEN
420  DO 130 k = 0, ib-1
421  IF( jj.GT.jja ) THEN
422  DO 120 ll = (jja-1)*lda, (jj-2)*lda, lda
423  VALUE = max( VALUE, abs( a( ii+ll ) ) )
424  120 CONTINUE
425  END IF
426  ii = ii + 1
427  IF( mycol.EQ.iacol )
428  $ jj = jj + 1
429  130 CONTINUE
430  ELSE IF( mycol.EQ.iacol ) THEN
431  jj = jj + ib
432  END IF
433 *
434  icurrow = mod( iarow+1, nprow )
435  icurcol = mod( iacol+1, npcol )
436 *
437 * Loop over rows/columns of global matrix.
438 *
439  DO 180 i = in+1, ia+n-1, desca( mb_ )
440  ib = min( desca( mb_ ), ia+n-i )
441 *
442 * Find COLMAXS
443 *
444  IF( mycol.EQ.icurcol ) THEN
445  DO 150 k = (jj-1)*lda, (jj+ib-2)*lda, lda
446  IF( myrow.EQ.icurrow ) THEN
447  IF( ii.LE.iia+np-1 ) THEN
448  VALUE = max( VALUE,
449  $ abs( real( a( ii+k ) ) ) )
450  DO 140 ll = ii+1, iia+np-1
451  VALUE = max( VALUE, abs( a( ll+k ) ) )
452  140 CONTINUE
453  END IF
454  ELSE
455  IF( ii.LE.iia+np-1 ) THEN
456  DO 145 ll = ii, iia+np-1
457  VALUE = max( VALUE, abs( a( ll+k ) ) )
458  145 CONTINUE
459  END IF
460  END IF
461  IF( myrow.EQ.icurrow )
462  $ ii = ii + 1
463  150 CONTINUE
464 *
465 * Reset local indices so we can find ROWMAXS
466 *
467  IF( myrow.EQ.icurrow )
468  $ ii = ii - ib
469  END IF
470 *
471 * Find ROWMAXS
472 *
473  IF( myrow.EQ.icurrow ) THEN
474  DO 170 k = 0, ib-1
475  IF( jj.GT.jja ) THEN
476  DO 160 ll = (jja-1)*lda, (jj-2)*lda, lda
477  VALUE = max( VALUE, abs( a( ii+ll ) ) )
478  160 CONTINUE
479  END IF
480  ii = ii + 1
481  IF( mycol.EQ.icurcol )
482  $ jj = jj + 1
483  170 CONTINUE
484  ELSE IF( mycol.EQ.icurcol ) THEN
485  jj = jj + ib
486  END IF
487  icurrow = mod( icurrow+1, nprow )
488  icurcol = mod( icurcol+1, npcol )
489 *
490  180 CONTINUE
491 *
492  END IF
493 *
494 * Gather the result on process (IAROW,IACOL).
495 *
496  CALL sgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, i, k, -1,
497  $ iarow, iacol )
498 *
499  ELSE IF( lsame( norm, 'I' ) .OR. lsame( norm, 'O' ) .OR.
500  $ norm.EQ.'1' ) THEN
501 *
502 * Find normI( sub( A ) ) ( = norm1( sub( A ) ), since sub( A ) is
503 * hermitian).
504 *
505  IF( lsame( uplo, 'U' ) ) THEN
506 *
507 * Handle first block separately
508 *
509  ib = in-ia+1
510 *
511 * Find COLSUMS
512 *
513  IF( mycol.EQ.iacol ) THEN
514  ioffa = ( jj - 1 ) * lda
515  DO 200 k = 0, ib-1
516  sum = zero
517  IF( ii.GT.iia ) THEN
518  DO 190 ll = iia, ii-1
519  sum = sum + abs( a( ll+ioffa ) )
520  190 CONTINUE
521  END IF
522  ioffa = ioffa + lda
523  work( jj+k-jja+icsr0 ) = sum
524  IF( myrow.EQ.iarow )
525  $ ii = ii + 1
526  200 CONTINUE
527 *
528 * Reset local indices so we can find ROWSUMS
529 *
530  IF( myrow.EQ.iarow )
531  $ ii = ii - ib
532 *
533  END IF
534 *
535 * Find ROWSUMS
536 *
537  IF( myrow.EQ.iarow ) THEN
538  DO 220 k = ii, ii+ib-1
539  sum = zero
540  IF( mycol.EQ.iacol ) THEN
541  IF( jja+nq.GT.jj ) THEN
542  sum = abs( real( a( k+(jj-1)*lda ) ) )
543  DO 210 ll = jj*lda, (jja+nq-2)*lda, lda
544  sum = sum + abs( a( k+ll ) )
545  210 CONTINUE
546  END IF
547  ELSE
548  IF( jja+nq.GT.jj ) THEN
549  DO 215 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
550  sum = sum + abs( a( k+ll ) )
551  215 CONTINUE
552  END IF
553  END IF
554  work( k-iia+irsc0 ) = sum
555  IF( mycol.EQ.iacol )
556  $ jj = jj + 1
557  220 CONTINUE
558  ii = ii + ib
559  ELSE IF( mycol.EQ.iacol ) THEN
560  jj = jj + ib
561  END IF
562 *
563  icurrow = mod( iarow+1, nprow )
564  icurcol = mod( iacol+1, npcol )
565 *
566 * Loop over remaining rows/columns of global matrix.
567 *
568  DO 270 i = in+1, ia+n-1, desca( mb_ )
569  ib = min( desca( mb_ ), ia+n-i )
570 *
571 * Find COLSUMS
572 *
573  IF( mycol.EQ.icurcol ) THEN
574  ioffa = ( jj - 1 ) * lda
575  DO 240 k = 0, ib-1
576  sum = zero
577  IF( ii.GT.iia ) THEN
578  DO 230 ll = iia, ii-1
579  sum = sum + abs( a( ioffa+ll ) )
580  230 CONTINUE
581  END IF
582  ioffa = ioffa + lda
583  work( jj+k-jja+icsr0 ) = sum
584  IF( myrow.EQ.icurrow )
585  $ ii = ii + 1
586  240 CONTINUE
587 *
588 * Reset local indices so we can find ROWSUMS
589 *
590  IF( myrow.EQ.icurrow )
591  $ ii = ii - ib
592 *
593  END IF
594 *
595 * Find ROWSUMS
596 *
597  IF( myrow.EQ.icurrow ) THEN
598  DO 260 k = ii, ii+ib-1
599  sum = zero
600  IF( mycol.EQ.icurcol ) THEN
601  IF( jja+nq.GT.jj ) THEN
602  sum = abs( real( a( k+(jj-1)*lda ) ) )
603  DO 250 ll = jj*lda, (jja+nq-2)*lda, lda
604  sum = sum + abs( a( k+ll ) )
605  250 CONTINUE
606  END IF
607  ELSE
608  IF( jja+nq.GT.jj ) THEN
609  DO 255 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
610  sum = sum + abs( a( k+ll ) )
611  255 CONTINUE
612  END IF
613  END IF
614  work( k-iia+irsc0 ) = sum
615  IF( mycol.EQ.icurcol )
616  $ jj = jj + 1
617  260 CONTINUE
618  ii = ii + ib
619  ELSE IF( mycol.EQ.icurcol ) THEN
620  jj = jj + ib
621  END IF
622 *
623  icurrow = mod( icurrow+1, nprow )
624  icurcol = mod( icurcol+1, npcol )
625 *
626  270 CONTINUE
627 *
628  ELSE
629 *
630 * Handle first block separately
631 *
632  ib = in-ia+1
633 *
634 * Find COLSUMS
635 *
636  IF( mycol.EQ.iacol ) THEN
637  ioffa = (jj-1)*lda
638  DO 290 k = 0, ib-1
639  sum = zero
640  IF( myrow.EQ.iarow ) THEN
641  IF( iia+np.GT.ii ) THEN
642  sum = abs( real( a( ioffa+ii ) ) )
643  DO 280 ll = ii+1, iia+np-1
644  sum = sum + abs( a( ioffa+ll ) )
645  280 CONTINUE
646  END IF
647  ELSE
648  DO 285 ll = ii, iia+np-1
649  sum = sum + abs( a( ioffa+ll ) )
650  285 CONTINUE
651  END IF
652  ioffa = ioffa + lda
653  work( jj+k-jja+icsr0 ) = sum
654  IF( myrow.EQ.iarow )
655  $ ii = ii + 1
656  290 CONTINUE
657 *
658 * Reset local indices so we can find ROWSUMS
659 *
660  IF( myrow.EQ.iarow )
661  $ ii = ii - ib
662 *
663  END IF
664 *
665 * Find ROWSUMS
666 *
667  IF( myrow.EQ.iarow ) THEN
668  DO 310 k = ii, ii+ib-1
669  sum = zero
670  IF( jj.GT.jja ) THEN
671  DO 300 ll = (jja-1)*lda, (jj-2)*lda, lda
672  sum = sum + abs( a( k+ll ) )
673  300 CONTINUE
674  END IF
675  work( k-iia+irsc0 ) = sum
676  IF( mycol.EQ.iacol )
677  $ jj = jj + 1
678  310 CONTINUE
679  ii = ii + ib
680  ELSE IF( mycol.EQ.iacol ) THEN
681  jj = jj + ib
682  END IF
683 *
684  icurrow = mod( iarow+1, nprow )
685  icurcol = mod( iacol+1, npcol )
686 *
687 * Loop over rows/columns of global matrix.
688 *
689  DO 360 i = in+1, ia+n-1, desca( mb_ )
690  ib = min( desca( mb_ ), ia+n-i )
691 *
692 * Find COLSUMS
693 *
694  IF( mycol.EQ.icurcol ) THEN
695  ioffa = ( jj - 1 ) * lda
696  DO 330 k = 0, ib-1
697  sum = zero
698  IF( myrow.EQ.icurrow ) THEN
699  IF( iia+np.GT.ii ) THEN
700  sum = abs( real( a( ii+ioffa ) ) )
701  DO 320 ll = ii+1, iia+np-1
702  sum = sum + abs( a( ll+ioffa ) )
703  320 CONTINUE
704  ELSE IF( ii.EQ.iia+np-1 ) THEN
705  sum = abs( real( a( ii+ioffa ) ) )
706  END IF
707  ELSE
708  DO 325 ll = ii, iia+np-1
709  sum = sum + abs( a( ll+ioffa ) )
710  325 CONTINUE
711  END IF
712  ioffa = ioffa + lda
713  work( jj+k-jja+icsr0 ) = sum
714  IF( myrow.EQ.icurrow )
715  $ ii = ii + 1
716  330 CONTINUE
717 *
718 * Reset local indices so we can find ROWSUMS
719 *
720  IF( myrow.EQ.icurrow )
721  $ ii = ii - ib
722 *
723  END IF
724 *
725 * Find ROWSUMS
726 *
727  IF( myrow.EQ.icurrow ) THEN
728  DO 350 k = ii, ii+ib-1
729  sum = zero
730  IF( jj.GT.jja ) THEN
731  DO 340 ll = (jja-1)*lda, (jj-2)*lda, lda
732  sum = sum + abs( a( k+ll ) )
733  340 CONTINUE
734  END IF
735  work(k-iia+irsc0) = sum
736  IF( mycol.EQ.icurcol )
737  $ jj = jj + 1
738  350 CONTINUE
739  ii = ii + ib
740  ELSE IF( mycol.EQ.icurcol ) THEN
741  jj = jj + ib
742  END IF
743 *
744  icurrow = mod( icurrow+1, nprow )
745  icurcol = mod( icurcol+1, npcol )
746 *
747  360 CONTINUE
748  END IF
749 *
750 * After calls to SGSUM2D, process row 0 will have global
751 * COLSUMS and process column 0 will have global ROWSUMS.
752 * Transpose ROWSUMS and add to COLSUMS to get global row/column
753 * sum, the max of which is the infinity or 1 norm.
754 *
755  IF( mycol.EQ.iacol )
756  $ nq = nq + icoff
757  CALL sgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work( icsr ), 1,
758  $ iarow, mycol )
759  IF( myrow.EQ.iarow )
760  $ np = np + iroff
761  CALL sgsum2d( ictxt, 'Rowwise', ' ', np, 1, work( irsc ),
762  $ max( 1, np ), myrow, iacol )
763 *
764  CALL pscol2row( ictxt, n, 1, desca( mb_ ), work( irsc ),
765  $ max( 1, np ), work( irsr ), max( 1, nq ),
766  $ iarow, iacol, iarow, iacol, work( irsc+np ) )
767 *
768  IF( myrow.EQ.iarow ) THEN
769  IF( mycol.EQ.iacol )
770  $ nq = nq - icoff
771  CALL saxpy( nq, one, work( irsr0 ), 1, work( icsr0 ), 1 )
772  IF( nq.LT.1 ) THEN
773  VALUE = zero
774  ELSE
775  VALUE = work( isamax( nq, work( icsr0 ), 1 ) )
776  END IF
777  CALL sgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, i, k,
778  $ -1, iarow, iacol )
779  END IF
780 *
781  ELSE IF( lsame( norm, 'F' ) .OR. lsame( norm, 'E' ) ) THEN
782 *
783 * Find normF( sub( A ) ).
784 *
785  scale = zero
786  sum = one
787 *
788 * Add off-diagonal entries, first
789 *
790  IF( lsame( uplo, 'U' ) ) THEN
791 *
792 * Handle first block separately
793 *
794  ib = in-ia+1
795 *
796  IF( mycol.EQ.iacol ) THEN
797  DO 370 k = (jj-1)*lda, (jj+ib-2)*lda, lda
798  CALL classq( ii-iia, a( iia+k ), 1, scale, sum )
799  CALL classq( ii-iia, a( iia+k ), 1, scale, sum )
800  IF( myrow.EQ.iarow ) THEN
801  IF( real( a( ii+k ) ).NE.zero ) THEN
802  absa = abs( real( a( ii+k ) ) )
803  IF( scale.LT.absa ) THEN
804  sum = one + sum * ( scale / absa )**2
805  scale = absa
806  ELSE
807  sum = sum + ( absa / scale )**2
808  END IF
809  END IF
810  ii = ii + 1
811  END IF
812  370 CONTINUE
813 *
814  jj = jj + ib
815  ELSE IF( myrow.EQ.iarow ) THEN
816  ii = ii + ib
817  END IF
818 *
819  icurrow = mod( iarow+1, nprow )
820  icurcol = mod( iacol+1, npcol )
821 *
822 * Loop over rows/columns of global matrix.
823 *
824  DO 390 i = in+1, ia+n-1, desca( mb_ )
825  ib = min( desca( mb_ ), ia+n-i )
826 *
827  IF( mycol.EQ.icurcol ) THEN
828  DO 380 k = (jj-1)*lda, (jj+ib-2)*lda, lda
829  CALL classq( ii-iia, a( iia+k ), 1, scale, sum )
830  CALL classq( ii-iia, a( iia+k ), 1, scale, sum )
831  IF( myrow.EQ.icurrow ) THEN
832  IF( real( a( ii+k ) ).NE.zero ) THEN
833  absa = abs( real( a( ii+k ) ) )
834  IF( scale.LT.absa ) THEN
835  sum = one + sum*( scale / absa )**2
836  scale = absa
837  ELSE
838  sum = sum + ( absa / scale )**2
839  END IF
840  END IF
841  ii = ii + 1
842  END IF
843  380 CONTINUE
844 *
845  jj = jj + ib
846  ELSE IF( myrow.EQ.icurrow ) THEN
847  ii = ii + ib
848  END IF
849 *
850  icurrow = mod( icurrow+1, nprow )
851  icurcol = mod( icurcol+1, npcol )
852 *
853  390 CONTINUE
854 *
855  ELSE
856 *
857 * Handle first block separately
858 *
859  ib = in-ia+1
860 *
861  IF( mycol.EQ.iacol ) THEN
862  DO 400 k = (jj-1)*lda, (jj+ib-2)*lda, lda
863  IF( myrow.EQ.iarow ) THEN
864  IF( real( a( ii+k ) ).NE.zero ) THEN
865  absa = abs( real( a( ii+k ) ) )
866  IF( scale.LT.absa ) THEN
867  sum = one + sum * ( scale / absa )**2
868  scale = absa
869  ELSE
870  sum = sum + ( absa / scale )**2
871  END IF
872  END IF
873  ii = ii + 1
874  END IF
875  CALL classq( iia+np-ii, a( ii+k ), 1, scale, sum )
876  CALL classq( iia+np-ii, a( ii+k ), 1, scale, sum )
877  400 CONTINUE
878 *
879  jj = jj + ib
880  ELSE IF( myrow.EQ.iarow ) THEN
881  ii = ii + ib
882  END IF
883 *
884  icurrow = mod( iarow+1, nprow )
885  icurcol = mod( iacol+1, npcol )
886 *
887 * Loop over rows/columns of global matrix.
888 *
889  DO 420 i = in+1, ia+n-1, desca( mb_ )
890  ib = min( desca( mb_ ), ia+n-i )
891 *
892  IF( mycol.EQ.icurcol ) THEN
893  DO 410 k = (jj-1)*lda, (jj+ib-2)*lda, lda
894  IF( myrow.EQ.icurrow ) THEN
895  IF( real( a( ii+k ) ).NE.zero ) THEN
896  absa = abs( real( a( ii+k ) ) )
897  IF( scale.LT.absa ) THEN
898  sum = one + sum * ( scale / absa )**2
899  scale = absa
900  ELSE
901  sum = sum + ( absa / scale )**2
902  END IF
903  END IF
904  ii = ii + 1
905  END IF
906  CALL classq( iia+np-ii, a( ii+k ), 1, scale, sum )
907  CALL classq( iia+np-ii, a( ii+k ), 1, scale, sum )
908  410 CONTINUE
909 *
910  jj = jj + ib
911  ELSE IF( myrow.EQ.icurrow ) THEN
912  ii = ii + ib
913  END IF
914 *
915  icurrow = mod( icurrow+1, nprow )
916  icurcol = mod( icurcol+1, npcol )
917 *
918  420 CONTINUE
919 *
920  END IF
921 *
922 * Perform the global scaled sum
923 *
924  rwork( 1 ) = scale
925  rwork( 2 ) = sum
926 *
927  CALL pstreecomb( ictxt, 'All', 2, rwork, iarow, iacol,
928  $ scombssq )
929  VALUE = rwork( 1 ) * sqrt( rwork( 2 ) )
930 *
931  END IF
932 *
933 * Broadcast the result to the other processes
934 *
935  IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) THEN
936  CALL sgebs2d( ictxt, 'All', ' ', 1, 1, VALUE, 1 )
937  ELSE
938  CALL sgebr2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, iarow,
939  $ iacol )
940  END IF
941 *
942  pclanhe = VALUE
943 *
944  RETURN
945 *
946 * End of PCLANHE
947 *
948  END
max
#define max(A, B)
Definition: pcgemr.c:180
infog2l
subroutine infog2l(GRINDX, GCINDX, DESC, NPROW, NPCOL, MYROW, MYCOL, LRINDX, LCINDX, RSRC, CSRC)
Definition: infog2l.f:3
pscol2row
subroutine pscol2row(ICTXT, M, N, NB, VS, LDVS, VD, LDVD, RSRC, CSRC, RDEST, CDEST, WORK)
Definition: pscol2row.f:3
pstreecomb
subroutine pstreecomb(ICTXT, SCOPE, N, MINE, RDEST0, CDEST0, SUBPTR)
Definition: pstreecomb.f:3
numroc
integer function numroc(N, NB, IPROC, ISRCPROC, NPROCS)
Definition: numroc.f:2
scombssq
subroutine scombssq(V1, V2)
Definition: pstreecomb.f:258
min
#define min(A, B)
Definition: pcgemr.c:181
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2