ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pslansy.f
Go to the documentation of this file.
1  REAL FUNCTION PSLANSY( 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  REAL a( * ), work( * )
17 * ..
18 *
19 * Purpose
20 * =======
21 *
22 * PSLANSY 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 * real symmetric distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
25 *
26 * PSLANSY 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 PSLANSY as described
101 * above.
102 *
103 * UPLO (global input) CHARACTER
104 * Specifies whether the upper or lower triangular part of the
105 * symmetric 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, PSLANSY is set to zero. N >= 0.
113 *
114 * A (local input) REAL pointer into the local memory
115 * to an array of dimension (LLD_A, LOCc(JA+N-1)) containing the
116 * local pieces of the symmetric 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 sum, value
178 * ..
179 * .. Local Arrays ..
180  REAL ssq( 2 ), colssq( 2 )
181 * ..
182 * .. External Subroutines ..
183  EXTERNAL blacs_gridinfo, pscol2row, pstreecomb,
184  $ saxpy, scombssq, sgamx2d, sgsum2d,
185  $ sgebr2d, sgebs2d, slassq
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, 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 symmetric, 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 ************************************************************************
268 * max norm
269 *
270  ELSE IF( lsame( norm, 'M' ) ) THEN
271 *
272 * Find max(abs(A(i,j))).
273 *
274  VALUE = zero
275 *
276  IF( lsame( uplo, 'U' ) ) THEN
277 *
278 * Handle first block separately
279 *
280  ib = in-ia+1
281 *
282 * Find COLMAXS
283 *
284  IF( mycol.EQ.iacol ) THEN
285  DO 20 k = (jj-1)*lda, (jj+ib-2)*lda, lda
286  IF( ii.GT.iia ) THEN
287  DO 10 ll = iia, ii-1
288  VALUE = max( VALUE, abs( a( ll+k ) ) )
289  10 CONTINUE
290  END IF
291  IF( myrow.EQ.iarow )
292  $ ii = ii + 1
293  20 CONTINUE
294 *
295 * Reset local indices so we can find ROWMAXS
296 *
297  IF( myrow.EQ.iarow )
298  $ ii = ii - ib
299 *
300  END IF
301 *
302 * Find ROWMAXS
303 *
304  IF( myrow.EQ.iarow ) THEN
305  DO 40 k = ii, ii+ib-1
306  IF( jj.LE.jja+nq-1 ) THEN
307  DO 30 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
308  VALUE = max( VALUE, abs( a( k+ll ) ) )
309  30 CONTINUE
310  END IF
311  IF( mycol.EQ.iacol )
312  $ jj = jj + 1
313  40 CONTINUE
314  ii = ii + ib
315  ELSE IF( mycol.EQ.iacol ) THEN
316  jj = jj + ib
317  END IF
318 *
319  icurrow = mod( iarow+1, nprow )
320  icurcol = mod( iacol+1, npcol )
321 *
322 * Loop over the remaining rows/columns of the matrix.
323 *
324  DO 90 i = in+1, ia+n-1, desca( mb_ )
325  ib = min( desca( mb_ ), ia+n-i )
326 *
327 * Find COLMAXS
328 *
329  IF( mycol.EQ.icurcol ) THEN
330  DO 60 k = (jj-1)*lda, (jj+ib-2)*lda, lda
331  IF( ii.GT.iia ) THEN
332  DO 50 ll = iia, ii-1
333  VALUE = max( VALUE, abs( a( ll+k ) ) )
334  50 CONTINUE
335  END IF
336  IF( myrow.EQ.icurrow )
337  $ ii = ii + 1
338  60 CONTINUE
339 *
340 * Reset local indices so we can find ROWMAXS
341 *
342  IF( myrow.EQ.icurrow )
343  $ ii = ii - ib
344  END IF
345 *
346 * Find ROWMAXS
347 *
348  IF( myrow.EQ.icurrow ) THEN
349  DO 80 k = ii, ii+ib-1
350  IF( jj.LE.jja+nq-1 ) THEN
351  DO 70 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
352  VALUE = max( VALUE, abs( a( k+ll ) ) )
353  70 CONTINUE
354  END IF
355  IF( mycol.EQ.icurcol )
356  $ jj = jj + 1
357  80 CONTINUE
358  ii = ii + ib
359  ELSE IF( mycol.EQ.icurcol ) THEN
360  jj = jj + ib
361  END IF
362  icurrow = mod( icurrow+1, nprow )
363  icurcol = mod( icurcol+1, npcol )
364  90 CONTINUE
365 *
366  ELSE
367 *
368 * Handle first block separately
369 *
370  ib = in-ia+1
371 *
372 * Find COLMAXS
373 *
374  IF( mycol.EQ.iacol ) THEN
375  DO 110 k = (jj-1)*lda, (jj+ib-2)*lda, lda
376  IF( ii.LE.iia+np-1 ) THEN
377  DO 100 ll = ii, iia+np-1
378  VALUE = max( VALUE, abs( a( ll+k ) ) )
379  100 CONTINUE
380  END IF
381  IF( myrow.EQ.iarow )
382  $ ii = ii + 1
383  110 CONTINUE
384 *
385 * Reset local indices so we can find ROWMAXS
386 *
387  IF( myrow.EQ.iarow )
388  $ ii = ii - ib
389  END IF
390 *
391 * Find ROWMAXS
392 *
393  IF( myrow.EQ.iarow ) THEN
394  DO 130 k = 0, ib-1
395  IF( jj.GT.jja ) THEN
396  DO 120 ll = (jja-1)*lda, (jj-2)*lda, lda
397  VALUE = max( VALUE, abs( a( ii+ll ) ) )
398  120 CONTINUE
399  END IF
400  ii = ii + 1
401  IF( mycol.EQ.iacol )
402  $ jj = jj + 1
403  130 CONTINUE
404  ELSE IF( mycol.EQ.iacol ) THEN
405  jj = jj + ib
406  END IF
407 *
408  icurrow = mod( iarow+1, nprow )
409  icurcol = mod( iacol+1, npcol )
410 *
411 * Loop over rows/columns of global matrix.
412 *
413  DO 180 i = in+1, ia+n-1, desca( mb_ )
414  ib = min( desca( mb_ ), ia+n-i )
415 *
416 * Find COLMAXS
417 *
418  IF( mycol.EQ.icurcol ) THEN
419  DO 150 k = (jj-1)*lda, (jj+ib-2)*lda, lda
420  IF( ii.LE.iia+np-1 ) THEN
421  DO 140 ll = ii, iia+np-1
422  VALUE = max( VALUE, abs( a( ll+k ) ) )
423  140 CONTINUE
424  END IF
425  IF( myrow.EQ.icurrow )
426  $ ii = ii + 1
427  150 CONTINUE
428 *
429 * Reset local indices so we can find ROWMAXS
430 *
431  IF( myrow.EQ.icurrow )
432  $ ii = ii - ib
433  END IF
434 *
435 * Find ROWMAXS
436 *
437  IF( myrow.EQ.icurrow ) THEN
438  DO 170 k = 0, ib-1
439  IF( jj.GT.jja ) THEN
440  DO 160 ll = (jja-1)*lda, (jj-2)*lda, lda
441  VALUE = max( VALUE, abs( a( ii+ll ) ) )
442  160 CONTINUE
443  END IF
444  ii = ii + 1
445  IF( mycol.EQ.icurcol )
446  $ jj = jj + 1
447  170 CONTINUE
448  ELSE IF( mycol.EQ.icurcol ) THEN
449  jj = jj + ib
450  END IF
451  icurrow = mod( icurrow+1, nprow )
452  icurcol = mod( icurcol+1, npcol )
453 *
454  180 CONTINUE
455 *
456  END IF
457 *
458 * Gather the result on process (IAROW,IACOL).
459 *
460  CALL sgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, i, k, -1,
461  $ iarow, iacol )
462 *
463 ************************************************************************
464 * one or inf norm
465 *
466  ELSE IF( lsame( norm, 'I' ) .OR. lsame( norm, 'O' ) .OR.
467  $ norm.EQ.'1' ) THEN
468 *
469 * Find normI( sub( A ) ) ( = norm1( sub( A ) ), since sub( A ) is
470 * symmetric).
471 *
472  IF( lsame( uplo, 'U' ) ) THEN
473 *
474 * Handle first block separately
475 *
476  ib = in-ia+1
477 *
478 * Find COLSUMS
479 *
480  IF( mycol.EQ.iacol ) THEN
481  ioffa = ( jj - 1 ) * lda
482  DO 200 k = 0, ib-1
483  sum = zero
484  IF( ii.GT.iia ) THEN
485  DO 190 ll = iia, ii-1
486  sum = sum + abs( a( ll+ioffa ) )
487  190 CONTINUE
488  END IF
489  ioffa = ioffa + lda
490  work( jj+k-jja+icsr0 ) = sum
491  IF( myrow.EQ.iarow )
492  $ ii = ii + 1
493  200 CONTINUE
494 *
495 * Reset local indices so we can find ROWSUMS
496 *
497  IF( myrow.EQ.iarow )
498  $ ii = ii - ib
499 *
500  END IF
501 *
502 * Find ROWSUMS
503 *
504  IF( myrow.EQ.iarow ) THEN
505  DO 220 k = ii, ii+ib-1
506  sum = zero
507  IF( jja+nq.GT.jj ) THEN
508  DO 210 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
509  sum = sum + abs( a( k+ll ) )
510  210 CONTINUE
511  END IF
512  work( k-iia+irsc0 ) = sum
513  IF( mycol.EQ.iacol )
514  $ jj = jj + 1
515  220 CONTINUE
516  ii = ii + ib
517  ELSE IF( mycol.EQ.iacol ) THEN
518  jj = jj + ib
519  END IF
520 *
521  icurrow = mod( iarow+1, nprow )
522  icurcol = mod( iacol+1, npcol )
523 *
524 * Loop over remaining rows/columns of global matrix.
525 *
526  DO 270 i = in+1, ia+n-1, desca( mb_ )
527  ib = min( desca( mb_ ), ia+n-i )
528 *
529 * Find COLSUMS
530 *
531  IF( mycol.EQ.icurcol ) THEN
532  ioffa = ( jj - 1 ) * lda
533  DO 240 k = 0, ib-1
534  sum = zero
535  IF( ii.GT.iia ) THEN
536  DO 230 ll = iia, ii-1
537  sum = sum + abs( a( ioffa+ll ) )
538  230 CONTINUE
539  END IF
540  ioffa = ioffa + lda
541  work( jj+k-jja+icsr0 ) = sum
542  IF( myrow.EQ.icurrow )
543  $ ii = ii + 1
544  240 CONTINUE
545 *
546 * Reset local indices so we can find ROWSUMS
547 *
548  IF( myrow.EQ.icurrow )
549  $ ii = ii - ib
550 *
551  END IF
552 *
553 * Find ROWSUMS
554 *
555  IF( myrow.EQ.icurrow ) THEN
556  DO 260 k = ii, ii+ib-1
557  sum = zero
558  IF( jja+nq.GT.jj ) THEN
559  DO 250 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
560  sum = sum + abs( a( k+ll ) )
561  250 CONTINUE
562  END IF
563  work( k-iia+irsc0 ) = sum
564  IF( mycol.EQ.icurcol )
565  $ jj = jj + 1
566  260 CONTINUE
567  ii = ii + ib
568  ELSE IF( mycol.EQ.icurcol ) THEN
569  jj = jj + ib
570  END IF
571 *
572  icurrow = mod( icurrow+1, nprow )
573  icurcol = mod( icurcol+1, npcol )
574 *
575  270 CONTINUE
576 *
577  ELSE
578 *
579 * Handle first block separately
580 *
581  ib = in-ia+1
582 *
583 * Find COLSUMS
584 *
585  IF( mycol.EQ.iacol ) THEN
586  ioffa = (jj-1)*lda
587  DO 290 k = 0, ib-1
588  sum = zero
589  IF( iia+np.GT.ii ) THEN
590  DO 280 ll = ii, iia+np-1
591  sum = sum + abs( a( ioffa+ll ) )
592  280 CONTINUE
593  END IF
594  ioffa = ioffa + lda
595  work( jj+k-jja+icsr0 ) = sum
596  IF( myrow.EQ.iarow )
597  $ ii = ii + 1
598  290 CONTINUE
599 *
600 * Reset local indices so we can find ROWSUMS
601 *
602  IF( myrow.EQ.iarow )
603  $ ii = ii - ib
604 *
605  END IF
606 *
607 * Find ROWSUMS
608 *
609  IF( myrow.EQ.iarow ) THEN
610  DO 310 k = ii, ii+ib-1
611  sum = zero
612  IF( jj.GT.jja ) THEN
613  DO 300 ll = (jja-1)*lda, (jj-2)*lda, lda
614  sum = sum + abs( a( k+ll ) )
615  300 CONTINUE
616  END IF
617  work( k-iia+irsc0 ) = sum
618  IF( mycol.EQ.iacol )
619  $ jj = jj + 1
620  310 CONTINUE
621  ii = ii + ib
622  ELSE IF( mycol.EQ.iacol ) THEN
623  jj = jj + ib
624  END IF
625 *
626  icurrow = mod( iarow+1, nprow )
627  icurcol = mod( iacol+1, npcol )
628 *
629 * Loop over rows/columns of global matrix.
630 *
631  DO 360 i = in+1, ia+n-1, desca( mb_ )
632  ib = min( desca( mb_ ), ia+n-i )
633 *
634 * Find COLSUMS
635 *
636  IF( mycol.EQ.icurcol ) THEN
637  ioffa = ( jj - 1 ) * lda
638  DO 330 k = 0, ib-1
639  sum = zero
640  IF( iia+np.GT.ii ) THEN
641  DO 320 ll = ii, iia+np-1
642  sum = sum + abs( a( ll+ioffa ) )
643  320 CONTINUE
644  END IF
645  ioffa = ioffa + lda
646  work( jj+k-jja+icsr0 ) = sum
647  IF( myrow.EQ.icurrow )
648  $ ii = ii + 1
649  330 CONTINUE
650 *
651 * Reset local indices so we can find ROWSUMS
652 *
653  IF( myrow.EQ.icurrow )
654  $ ii = ii - ib
655 *
656  END IF
657 *
658 * Find ROWSUMS
659 *
660  IF( myrow.EQ.icurrow ) THEN
661  DO 350 k = ii, ii+ib-1
662  sum = zero
663  IF( jj.GT.jja ) THEN
664  DO 340 ll = (jja-1)*lda, (jj-2)*lda, lda
665  sum = sum + abs( a( k+ll ) )
666  340 CONTINUE
667  END IF
668  work(k-iia+irsc0) = sum
669  IF( mycol.EQ.icurcol )
670  $ jj = jj + 1
671  350 CONTINUE
672  ii = ii + ib
673  ELSE IF( mycol.EQ.icurcol ) THEN
674  jj = jj + ib
675  END IF
676 *
677  icurrow = mod( icurrow+1, nprow )
678  icurcol = mod( icurcol+1, npcol )
679 *
680  360 CONTINUE
681  END IF
682 *
683 * After calls to SGSUM2D, process row 0 will have global
684 * COLSUMS and process column 0 will have global ROWSUMS.
685 * Transpose ROWSUMS and add to COLSUMS to get global row/column
686 * sum, the max of which is the infinity or 1 norm.
687 *
688  IF( mycol.EQ.iacol )
689  $ nq = nq + icoff
690  CALL sgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work( icsr ), 1,
691  $ iarow, mycol )
692  IF( myrow.EQ.iarow )
693  $ np = np + iroff
694  CALL sgsum2d( ictxt, 'Rowwise', ' ', np, 1, work( irsc ),
695  $ max( 1, np ), myrow, iacol )
696 *
697  CALL pscol2row( ictxt, n, 1, desca( mb_ ), work( irsc ),
698  $ max( 1, np ), work( irsr ), max( 1, nq ),
699  $ iarow, iacol, iarow, iacol, work( irsc+np ) )
700 *
701  IF( myrow.EQ.iarow ) THEN
702  IF( mycol.EQ.iacol )
703  $ nq = nq - icoff
704  CALL saxpy( nq, one, work( irsr0 ), 1, work( icsr0 ), 1 )
705  IF( nq.LT.1 ) THEN
706  VALUE = zero
707  ELSE
708  VALUE = work( isamax( nq, work( icsr0 ), 1 ) )
709  END IF
710  CALL sgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, i, k,
711  $ -1, iarow, iacol )
712  END IF
713 *
714 ************************************************************************
715 * Frobenius norm
716 * SSQ(1) is scale
717 * SSQ(2) is sum-of-squares
718 *
719  ELSE IF( lsame( norm, 'F' ) .OR. lsame( norm, 'E' ) ) THEN
720 *
721 * Find normF( sub( A ) ).
722 *
723  ssq(1) = zero
724  ssq(2) = one
725 *
726 * Add off-diagonal entries, first
727 *
728  IF( lsame( uplo, 'U' ) ) THEN
729 *
730 * Handle first block separately
731 *
732  ib = in-ia+1
733 *
734  IF( mycol.EQ.iacol ) THEN
735  DO 370 k = (jj-1)*lda, (jj+ib-2)*lda, lda
736  colssq(1) = zero
737  colssq(2) = one
738  CALL slassq( ii-iia, a( iia+k ), 1,
739  $ colssq(1), colssq(2) )
740  IF( myrow.EQ.iarow )
741  $ ii = ii + 1
742  CALL slassq( ii-iia, a( iia+k ), 1,
743  $ colssq(1), colssq(2) )
744  CALL scombssq( ssq, colssq )
745  370 CONTINUE
746 *
747  jj = jj + ib
748  ELSE IF( myrow.EQ.iarow ) THEN
749  ii = ii + ib
750  END IF
751 *
752  icurrow = mod( iarow+1, nprow )
753  icurcol = mod( iacol+1, npcol )
754 *
755 * Loop over rows/columns of global matrix.
756 *
757  DO 390 i = in+1, ia+n-1, desca( mb_ )
758  ib = min( desca( mb_ ), ia+n-i )
759 *
760  IF( mycol.EQ.icurcol ) THEN
761  DO 380 k = (jj-1)*lda, (jj+ib-2)*lda, lda
762  colssq(1) = zero
763  colssq(2) = one
764  CALL slassq( ii-iia, a( iia+k ), 1,
765  $ colssq(1), colssq(2) )
766  IF( myrow.EQ.icurrow )
767  $ ii = ii + 1
768  CALL slassq( ii-iia, a(iia+k ), 1,
769  $ colssq(1), colssq(2) )
770  CALL scombssq( ssq, colssq )
771  380 CONTINUE
772 *
773  jj = jj + ib
774  ELSE IF( myrow.EQ.icurrow ) THEN
775  ii = ii + ib
776  END IF
777 *
778  icurrow = mod( icurrow+1, nprow )
779  icurcol = mod( icurcol+1, npcol )
780 *
781  390 CONTINUE
782 *
783  ELSE
784 *
785 * Handle first block separately
786 *
787  ib = in-ia+1
788 *
789  IF( mycol.EQ.iacol ) THEN
790  DO 400 k = (jj-1)*lda, (jj+ib-2)*lda, lda
791  colssq(1) = zero
792  colssq(2) = one
793  CALL slassq( iia+np-ii, a( ii+k ), 1,
794  $ colssq(1), colssq(2) )
795  IF( myrow.EQ.iarow )
796  $ ii = ii + 1
797  CALL slassq( iia+np-ii, a( ii+k ), 1,
798  $ colssq(1), colssq(2) )
799  CALL scombssq( ssq, colssq )
800  400 CONTINUE
801 *
802  jj = jj + ib
803  ELSE IF( myrow.EQ.iarow ) THEN
804  ii = ii + ib
805  END IF
806 *
807  icurrow = mod( iarow+1, nprow )
808  icurcol = mod( iacol+1, npcol )
809 *
810 * Loop over rows/columns of global matrix.
811 *
812  DO 420 i = in+1, ia+n-1, desca( mb_ )
813  ib = min( desca( mb_ ), ia+n-i )
814 *
815  IF( mycol.EQ.icurcol ) THEN
816  DO 410 k = (jj-1)*lda, (jj+ib-2)*lda, lda
817  colssq(1) = zero
818  colssq(2) = one
819  CALL slassq( iia+np-ii, a( ii+k ), 1,
820  $ colssq(1), colssq(2) )
821  IF( myrow.EQ.icurrow )
822  $ ii = ii + 1
823  CALL slassq( iia+np-ii, a( ii+k ), 1,
824  $ colssq(1), colssq(2) )
825  CALL scombssq( ssq, colssq )
826  410 CONTINUE
827 *
828  jj = jj + ib
829  ELSE IF( myrow.EQ.icurrow ) THEN
830  ii = ii + ib
831  END IF
832 *
833  icurrow = mod( icurrow+1, nprow )
834  icurcol = mod( icurcol+1, npcol )
835 *
836  420 CONTINUE
837 *
838  END IF
839 *
840 * Perform the global scaled sum
841 *
842  CALL pstreecomb( ictxt, 'All', 2, ssq, iarow, iacol,
843  $ scombssq )
844  VALUE = ssq( 1 ) * sqrt( ssq( 2 ) )
845 *
846  END IF
847 *
848 * Broadcast the result to the other processes
849 *
850  IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) THEN
851  CALL sgebs2d( ictxt, 'All', ' ', 1, 1, VALUE, 1 )
852  ELSE
853  CALL sgebr2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, iarow,
854  $ iacol )
855  END IF
856 *
857  pslansy = VALUE
858 *
859  RETURN
860 *
861 * End of PSLANSY
862 *
863  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
lsame
logical function lsame(CA, CB)
Definition: tools.f:1724
pslansy
real function pslansy(NORM, UPLO, N, A, IA, JA, DESCA, WORK)
Definition: pslansy.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