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