ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
pzlantr.f
Go to the documentation of this file.
1  DOUBLE PRECISION FUNCTION pzlantr( NORM, UPLO, DIAG, M, N, A,
2  $ IA, JA, 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 diag, norm, uplo
12  INTEGER ia, ja, m, n
13 * ..
14 * .. Array Arguments ..
15  INTEGER desca( * )
16  DOUBLE PRECISION work( * )
17  COMPLEX*16 a( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * PZLANTR returns the value of the one norm, or the Frobenius norm,
24 * or the infinity norm, or the element of largest absolute value of a
25 * trapezoidal or triangular distributed matrix sub( A ) denoting
26 * A(IA:IA+M-1, JA:JA+N-1).
27 *
28 * PZLANTR returns the value
29 *
30 * ( max(abs(A(i,j))), NORM = 'M' or 'm' with ia <= i <= ia+m-1,
31 * ( and ja <= j <= ja+n-1,
32 * (
33 * ( norm1( sub( A ) ), NORM = '1', 'O' or 'o'
34 * (
35 * ( normI( sub( A ) ), NORM = 'I' or 'i'
36 * (
37 * ( normF( sub( A ) ), NORM = 'F', 'f', 'E' or 'e'
38 *
39 * where norm1 denotes the one norm of a matrix (maximum column sum),
40 * normI denotes the infinity norm of a matrix (maximum row sum) and
41 * normF denotes the Frobenius norm of a matrix (square root of sum of
42 * squares). Note that max(abs(A(i,j))) is not a matrix norm.
43 *
44 * Notes
45 * =====
46 *
47 * Each global data object is described by an associated description
48 * vector. This vector stores the information required to establish
49 * the mapping between an object element and its corresponding process
50 * and memory location.
51 *
52 * Let A be a generic term for any 2D block cyclicly distributed array.
53 * Such a global array has an associated description vector DESCA.
54 * In the following comments, the character _ should be read as
55 * "of the global array".
56 *
57 * NOTATION STORED IN EXPLANATION
58 * --------------- -------------- --------------------------------------
59 * DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
60 * DTYPE_A = 1.
61 * CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
62 * the BLACS process grid A is distribu-
63 * ted over. The context itself is glo-
64 * bal, but the handle (the integer
65 * value) may vary.
66 * M_A (global) DESCA( M_ ) The number of rows in the global
67 * array A.
68 * N_A (global) DESCA( N_ ) The number of columns in the global
69 * array A.
70 * MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
71 * the rows of the array.
72 * NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
73 * the columns of the array.
74 * RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
75 * row of the array A is distributed.
76 * CSRC_A (global) DESCA( CSRC_ ) The process column over which the
77 * first column of the array A is
78 * distributed.
79 * LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
80 * array. LLD_A >= MAX(1,LOCr(M_A)).
81 *
82 * Let K be the number of rows or columns of a distributed matrix,
83 * and assume that its process grid has dimension p x q.
84 * LOCr( K ) denotes the number of elements of K that a process
85 * would receive if K were distributed over the p processes of its
86 * process column.
87 * Similarly, LOCc( K ) denotes the number of elements of K that a
88 * process would receive if K were distributed over the q processes of
89 * its process row.
90 * The values of LOCr() and LOCc() may be determined via a call to the
91 * ScaLAPACK tool function, NUMROC:
92 * LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
93 * LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
94 * An upper bound for these quantities may be computed by:
95 * LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
96 * LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
97 *
98 * Arguments
99 * =========
100 *
101 * NORM (global input) CHARACTER
102 * Specifies the value to be returned in PZLANTR as described
103 * above.
104 *
105 * UPLO (global input) CHARACTER
106 * Specifies whether the matrix sub( A ) is upper or lower
107 * trapezoidal.
108 * = 'U': Upper trapezoidal
109 * = 'L': Lower trapezoidal
110 * Note that sub( A ) is triangular instead of trapezoidal
111 * if M = N.
112 *
113 * DIAG (global input) CHARACTER
114 * Specifies whether or not the distributed matrix sub( A ) has
115 * unit diagonal.
116 * = 'N': Non-unit diagonal
117 * = 'U': Unit diagonal
118 *
119 * M (global input) INTEGER
120 * The number of rows to be operated on i.e the number of rows
121 * of the distributed submatrix sub( A ). When M = 0, PZLANTR is
122 * set to zero. M >= 0.
123 *
124 * N (global input) INTEGER
125 * The number of columns to be operated on i.e the number of
126 * columns of the distributed submatrix sub( A ). When N = 0,
127 * PZLANTR is set to zero. N >= 0.
128 *
129 * A (local input) COMPLEX*16 pointer into the local memory
130 * to an array of dimension (LLD_A, LOCc(JA+N-1) ) containing
131 * the local pieces of sub( A ).
132 *
133 * IA (global input) INTEGER
134 * The row index in the global array A indicating the first
135 * row of sub( A ).
136 *
137 * JA (global input) INTEGER
138 * The column index in the global array A indicating the
139 * first column of sub( A ).
140 *
141 * DESCA (global and local input) INTEGER array of dimension DLEN_.
142 * The array descriptor for the distributed matrix A.
143 *
144 * WORK (local workspace) DOUBLE PRECISION array dimension (LWORK)
145 * LWORK >= 0 if NORM = 'M' or 'm' (not referenced),
146 * Nq0 if NORM = '1', 'O' or 'o',
147 * Mp0 if NORM = 'I' or 'i',
148 * 0 if NORM = 'F', 'f', 'E' or 'e' (not referenced),
149 * where
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 * Mp0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
155 * Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
156 *
157 * INDXG2P and NUMROC are ScaLAPACK tool functions; MYROW,
158 * MYCOL, NPROW and NPCOL can be determined by calling the
159 * 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  LOGICAL udiag
174  INTEGER iacol, iarow, ictxt, ii, iia, icoff, ioffa,
175  $ iroff, j, jb, jj, jja, jn, kk, lda, ll, mp,
176  $ 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, dcombssq, dgebr2d,
184  $ dgebs2d, dgamx2d, dgsum2d, infog2l,
185  $ pdtreecomb, zlassq
186 * ..
187 * .. External Functions ..
188  LOGICAL lsame
189  INTEGER iceil, idamax, numroc
190  EXTERNAL lsame, iceil, idamax, numroc
191 * ..
192 * .. Intrinsic Functions ..
193  INTRINSIC abs, dble, max, min, mod, sqrt
194 * ..
195 * .. Executable Statements ..
196 *
197 * Get grid parameters
198 *
199  ictxt = desca( ctxt_ )
200  CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
201 *
202  udiag = lsame( diag, 'U' )
203  CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
204  $ iarow, iacol )
205  iroff = mod( ia-1, desca( mb_ ) )
206  icoff = mod( ja-1, desca( nb_ ) )
207  mp = numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
208  nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
209  IF( myrow.EQ.iarow )
210  $ mp = mp - iroff
211  IF( mycol.EQ.iacol )
212  $ nq = nq - icoff
213  lda = desca( lld_ )
214  ioffa = ( jja - 1 ) * lda
215 *
216  IF( min( m, n ).EQ.0 ) THEN
217 *
218  VALUE = zero
219 *
220 ************************************************************************
221 * max norm
222 *
223  ELSE IF( lsame( norm, 'M' ) ) THEN
224 *
225 * Find max(abs(A(i,j))).
226 *
227  IF( udiag ) THEN
228  VALUE = one
229  ELSE
230  VALUE = zero
231  END IF
232 *
233  IF( lsame( uplo, 'U' ) ) THEN
234 *
235 * Upper triangular matrix
236 *
237  ii = iia
238  jj = jja
239  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
240  jb = jn-ja+1
241 *
242  IF( mycol.EQ.iacol ) THEN
243  IF( myrow.EQ.iarow ) THEN
244  IF( udiag ) THEN
245  DO 20 ll = jj, jj + jb -1
246  DO 10 kk = iia, min(ii+ll-jj-1,iia+mp-1)
247  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
248  10 CONTINUE
249  ioffa = ioffa + lda
250  20 CONTINUE
251  ELSE
252  DO 40 ll = jj, jj + jb -1
253  DO 30 kk = iia, min( ii+ll-jj, iia+mp-1 )
254  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
255  30 CONTINUE
256  ioffa = ioffa + lda
257  40 CONTINUE
258  END IF
259  ELSE
260  DO 60 ll = jj, jj + jb -1
261  DO 50 kk = iia, min( ii-1, iia+mp-1 )
262  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
263  50 CONTINUE
264  ioffa = ioffa + lda
265  60 CONTINUE
266  END IF
267  jj = jj + jb
268  END IF
269 *
270  IF( myrow.EQ.iarow )
271  $ ii = ii + jb
272  iarow = mod( iarow+1, nprow )
273  iacol = mod( iacol+1, npcol )
274 *
275 * Loop over remaining block of columns
276 *
277  DO 130 j = jn+1, ja+n-1, desca( nb_ )
278  jb = min( ja+n-j, desca( nb_ ) )
279 *
280  IF( mycol.EQ.iacol ) THEN
281  IF( myrow.EQ.iarow ) THEN
282  IF( udiag ) THEN
283  DO 80 ll = jj, jj + jb -1
284  DO 70 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
285  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
286  70 CONTINUE
287  ioffa = ioffa + lda
288  80 CONTINUE
289  ELSE
290  DO 100 ll = jj, jj + jb -1
291  DO 90 kk = iia, min( ii+ll-jj, iia+mp-1 )
292  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
293  90 CONTINUE
294  ioffa = ioffa + lda
295  100 CONTINUE
296  END IF
297  ELSE
298  DO 120 ll = jj, jj + jb -1
299  DO 110 kk = iia, min( ii-1, iia+mp-1 )
300  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
301  110 CONTINUE
302  ioffa = ioffa + lda
303  120 CONTINUE
304  END IF
305  jj = jj + jb
306  END IF
307 *
308  IF( myrow.EQ.iarow )
309  $ ii = ii + jb
310  iarow = mod( iarow+1, nprow )
311  iacol = mod( iacol+1, npcol )
312 *
313  130 CONTINUE
314 *
315  ELSE
316 *
317 * Lower triangular matrix
318 *
319  ii = iia
320  jj = jja
321  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
322  jb = jn-ja+1
323 *
324  IF( mycol.EQ.iacol ) THEN
325  IF( myrow.EQ.iarow ) THEN
326  IF( udiag ) THEN
327  DO 150 ll = jj, jj + jb -1
328  DO 140 kk = ii+ll-jj+1, iia+mp-1
329  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
330  140 CONTINUE
331  ioffa = ioffa + lda
332  150 CONTINUE
333  ELSE
334  DO 170 ll = jj, jj + jb -1
335  DO 160 kk = ii+ll-jj, iia+mp-1
336  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
337  160 CONTINUE
338  ioffa = ioffa + lda
339  170 CONTINUE
340  END IF
341  ELSE
342  DO 190 ll = jj, jj + jb -1
343  DO 180 kk = ii, iia+mp-1
344  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
345  180 CONTINUE
346  ioffa = ioffa + lda
347  190 CONTINUE
348  END IF
349  jj = jj + jb
350  END IF
351 *
352  IF( myrow.EQ.iarow )
353  $ ii = ii + jb
354  iarow = mod( iarow+1, nprow )
355  iacol = mod( iacol+1, npcol )
356 *
357 * Loop over remaining block of columns
358 *
359  DO 260 j = jn+1, ja+n-1, desca( nb_ )
360  jb = min( ja+n-j, desca( nb_ ) )
361 *
362  IF( mycol.EQ.iacol ) THEN
363  IF( myrow.EQ.iarow ) THEN
364  IF( udiag ) THEN
365  DO 210 ll = jj, jj + jb -1
366  DO 200 kk = ii+ll-jj+1, iia+mp-1
367  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
368  200 CONTINUE
369  ioffa = ioffa + lda
370  210 CONTINUE
371  ELSE
372  DO 230 ll = jj, jj + jb -1
373  DO 220 kk = ii+ll-jj, iia+mp-1
374  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
375  220 CONTINUE
376  ioffa = ioffa + lda
377  230 CONTINUE
378  END IF
379  ELSE
380  DO 250 ll = jj, jj + jb -1
381  DO 240 kk = ii, iia+mp-1
382  VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
383  240 CONTINUE
384  ioffa = ioffa + lda
385  250 CONTINUE
386  END IF
387  jj = jj + jb
388  END IF
389 *
390  IF( myrow.EQ.iarow )
391  $ ii = ii + jb
392  iarow = mod( iarow+1, nprow )
393  iacol = mod( iacol+1, npcol )
394 *
395  260 CONTINUE
396 *
397  END IF
398 *
399 * Gather the intermediate results to process (0,0).
400 *
401  CALL dgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, kk, ll, -1,
402  $ 0, 0 )
403 *
404 ************************************************************************
405 * one norm
406 *
407  ELSE IF( lsame( norm, 'O' ) .OR. norm.EQ.'1' ) THEN
408 *
409  VALUE = zero
410 *
411  IF( lsame( uplo, 'U' ) ) THEN
412 *
413 * Upper triangular matrix
414 *
415  ii = iia
416  jj = jja
417  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
418  jb = jn-ja+1
419 *
420  IF( mycol.EQ.iacol ) THEN
421  IF( myrow.EQ.iarow ) THEN
422  IF( udiag ) THEN
423  DO 280 ll = jj, jj + jb -1
424  sum = zero
425  DO 270 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
426  sum = sum + abs( a( ioffa+kk ) )
427  270 CONTINUE
428 * Unit diagonal entry
429  kk = ii+ll-jj
430  IF (kk <= iia+mp-1) THEN
431  sum = sum + one
432  ENDIF
433  ioffa = ioffa + lda
434  work( ll-jja+1 ) = sum
435  280 CONTINUE
436  ELSE
437  DO 300 ll = jj, jj + jb -1
438  sum = zero
439  DO 290 kk = iia, min( ii+ll-jj, iia+mp-1 )
440  sum = sum + abs( a( ioffa+kk ) )
441  290 CONTINUE
442  ioffa = ioffa + lda
443  work( ll-jja+1 ) = sum
444  300 CONTINUE
445  END IF
446  ELSE
447  DO 320 ll = jj, jj + jb -1
448  sum = zero
449  DO 310 kk = iia, min( ii-1, iia+mp-1 )
450  sum = sum + abs( a( ioffa+kk ) )
451  310 CONTINUE
452  ioffa = ioffa + lda
453  work( ll-jja+1 ) = sum
454  320 CONTINUE
455  END IF
456  jj = jj + jb
457  END IF
458 *
459  IF( myrow.EQ.iarow )
460  $ ii = ii + jb
461  iarow = mod( iarow+1, nprow )
462  iacol = mod( iacol+1, npcol )
463 *
464 * Loop over remaining block of columns
465 *
466  DO 390 j = jn+1, ja+n-1, desca( nb_ )
467  jb = min( ja+n-j, desca( nb_ ) )
468 *
469  IF( mycol.EQ.iacol ) THEN
470  IF( myrow.EQ.iarow ) THEN
471  IF( udiag ) THEN
472  DO 340 ll = jj, jj + jb -1
473  sum = zero
474  DO 330 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
475  sum = sum + abs( a( ioffa+kk ) )
476  330 CONTINUE
477 * Unit diagonal entry
478  kk = ii+ll-jj
479  IF (kk <= iia+mp-1) THEN
480  sum = sum + one
481  ENDIF
482  ioffa = ioffa + lda
483  work( ll-jja+1 ) = sum
484  340 CONTINUE
485  ELSE
486  DO 360 ll = jj, jj + jb -1
487  sum = zero
488  DO 350 kk = iia, min( ii+ll-jj, iia+mp-1 )
489  sum = sum + abs( a( ioffa+kk ) )
490  350 CONTINUE
491  ioffa = ioffa + lda
492  work( ll-jja+1 ) = sum
493  360 CONTINUE
494  END IF
495  ELSE
496  DO 380 ll = jj, jj + jb -1
497  sum = zero
498  DO 370 kk = iia, min( ii-1, iia+mp-1 )
499  sum = sum + abs( a( ioffa+kk ) )
500  370 CONTINUE
501  ioffa = ioffa + lda
502  work( ll-jja+1 ) = sum
503  380 CONTINUE
504  END IF
505  jj = jj + jb
506  END IF
507 *
508  IF( myrow.EQ.iarow )
509  $ ii = ii + jb
510  iarow = mod( iarow+1, nprow )
511  iacol = mod( iacol+1, npcol )
512 *
513  390 CONTINUE
514 *
515  ELSE
516 *
517 * Lower triangular matrix
518 *
519  ii = iia
520  jj = jja
521  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
522  jb = jn-ja+1
523 *
524  IF( mycol.EQ.iacol ) THEN
525  IF( myrow.EQ.iarow ) THEN
526  IF( udiag ) THEN
527  DO 410 ll = jj, jj + jb -1
528  sum = one
529  DO 400 kk = ii+ll-jj+1, iia+mp-1
530  sum = sum + abs( a( ioffa+kk ) )
531  400 CONTINUE
532  ioffa = ioffa + lda
533  work( ll-jja+1 ) = sum
534  410 CONTINUE
535  ELSE
536  DO 430 ll = jj, jj + jb -1
537  sum = zero
538  DO 420 kk = ii+ll-jj, iia+mp-1
539  sum = sum + abs( a( ioffa+kk ) )
540  420 CONTINUE
541  ioffa = ioffa + lda
542  work( ll-jja+1 ) = sum
543  430 CONTINUE
544  END IF
545  ELSE
546  DO 450 ll = jj, jj + jb -1
547  sum = zero
548  DO 440 kk = ii, iia+mp-1
549  sum = sum + abs( a( ioffa+kk ) )
550  440 CONTINUE
551  ioffa = ioffa + lda
552  work( ll-jja+1 ) = sum
553  450 CONTINUE
554  END IF
555  jj = jj + jb
556  END IF
557 *
558  IF( myrow.EQ.iarow )
559  $ ii = ii + jb
560  iarow = mod( iarow+1, nprow )
561  iacol = mod( iacol+1, npcol )
562 *
563 * Loop over remaining block of columns
564 *
565  DO 520 j = jn+1, ja+n-1, desca( nb_ )
566  jb = min( ja+n-j, desca( nb_ ) )
567 *
568  IF( mycol.EQ.iacol ) THEN
569  IF( myrow.EQ.iarow ) THEN
570  IF( udiag ) THEN
571  DO 470 ll = jj, jj + jb -1
572  sum = one
573  DO 460 kk = ii+ll-jj+1, iia+mp-1
574  sum = sum + abs( a( ioffa+kk ) )
575  460 CONTINUE
576  ioffa = ioffa + lda
577  work( ll-jja+1 ) = sum
578  470 CONTINUE
579  ELSE
580  DO 490 ll = jj, jj + jb -1
581  sum = zero
582  DO 480 kk = ii+ll-jj, iia+mp-1
583  sum = sum + abs( a( ioffa+kk ) )
584  480 CONTINUE
585  ioffa = ioffa + lda
586  work( ll-jja+1 ) = sum
587  490 CONTINUE
588  END IF
589  ELSE
590  DO 510 ll = jj, jj + jb -1
591  sum = zero
592  DO 500 kk = ii, iia+mp-1
593  sum = sum + abs( a( ioffa+kk ) )
594  500 CONTINUE
595  ioffa = ioffa + lda
596  work( ll-jja+1 ) = sum
597  510 CONTINUE
598  END IF
599  jj = jj + jb
600  END IF
601 *
602  IF( myrow.EQ.iarow )
603  $ ii = ii + jb
604  iarow = mod( iarow+1, nprow )
605  iacol = mod( iacol+1, npcol )
606 *
607  520 CONTINUE
608 *
609  END IF
610 *
611 * Find sum of global matrix columns and store on row 0 of
612 * process grid
613 *
614  CALL dgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work, 1,
615  $ 0, mycol )
616 *
617 * Find maximum sum of columns for 1-norm
618 *
619  IF( myrow.EQ.0 ) THEN
620  IF( nq.GT.0 ) THEN
621  VALUE = work( idamax( nq, work, 1 ) )
622  ELSE
623  VALUE = zero
624  END IF
625  CALL dgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, kk, ll,
626  $ -1, 0, 0 )
627  END IF
628 *
629 ************************************************************************
630 * infinity norm
631 *
632  ELSE IF( lsame( norm, 'I' ) ) THEN
633 *
634  IF( lsame( uplo, 'U' ) ) THEN
635  DO 540 kk = iia, iia+mp-1
636  work( kk ) = zero
637  540 CONTINUE
638  ELSE
639  DO 570 kk = iia, iia+mp-1
640  work( kk ) = zero
641  570 CONTINUE
642  END IF
643 *
644  IF( lsame( uplo, 'U' ) ) THEN
645 *
646 * Upper triangular matrix
647 *
648  ii = iia
649  jj = jja
650  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
651  jb = jn-ja+1
652 *
653  IF( mycol.EQ.iacol ) THEN
654  IF( myrow.EQ.iarow ) THEN
655  IF( udiag ) THEN
656  DO 590 ll = jj, jj + jb -1
657  DO 580 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
658  work( kk-iia+1 ) = work( kk-iia+1 ) +
659  $ abs( a( ioffa+kk ) )
660  580 CONTINUE
661 * Unit diagonal entry
662  kk = ii+ll-jj
663  IF (kk <= iia+mp-1) THEN
664  work( kk-iia+1 ) = work( kk-iia+1 ) + one
665  ENDIF
666  ioffa = ioffa + lda
667  590 CONTINUE
668  ELSE
669  DO 610 ll = jj, jj + jb -1
670  DO 600 kk = iia, min( ii+ll-jj, iia+mp-1 )
671  work( kk-iia+1 ) = work( kk-iia+1 ) +
672  $ abs( a( ioffa+kk ) )
673  600 CONTINUE
674  ioffa = ioffa + lda
675  610 CONTINUE
676  END IF
677  ELSE
678  DO 630 ll = jj, jj + jb -1
679  DO 620 kk = iia, min( ii-1, iia+mp-1 )
680  work( kk-iia+1 ) = work( kk-iia+1 ) +
681  $ abs( a( ioffa+kk ) )
682  620 CONTINUE
683  ioffa = ioffa + lda
684  630 CONTINUE
685  END IF
686  jj = jj + jb
687  END IF
688 *
689  IF( myrow.EQ.iarow )
690  $ ii = ii + jb
691  iarow = mod( iarow+1, nprow )
692  iacol = mod( iacol+1, npcol )
693 *
694 * Loop over remaining block of columns
695 *
696  DO 700 j = jn+1, ja+n-1, desca( nb_ )
697  jb = min( ja+n-j, desca( nb_ ) )
698 *
699  IF( mycol.EQ.iacol ) THEN
700  IF( myrow.EQ.iarow ) THEN
701  IF( udiag ) THEN
702  DO 650 ll = jj, jj + jb -1
703  DO 640 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
704  work( kk-iia+1 ) = work( kk-iia+1 ) +
705  $ abs( a( ioffa+kk ) )
706  640 CONTINUE
707 * Unit diagonal entry
708  kk = ii+ll-jj
709  IF (kk <= iia+mp-1) THEN
710  work( kk-iia+1 ) = work( kk-iia+1 ) + one
711  ENDIF
712  ioffa = ioffa + lda
713  650 CONTINUE
714  ELSE
715  DO 670 ll = jj, jj + jb -1
716  DO 660 kk = iia, min( ii+ll-jj, iia+mp-1 )
717  work( kk-iia+1 ) = work( kk-iia+1 ) +
718  $ abs( a( ioffa+kk ) )
719  660 CONTINUE
720  ioffa = ioffa + lda
721  670 CONTINUE
722  END IF
723  ELSE
724  DO 690 ll = jj, jj + jb -1
725  DO 680 kk = iia, min( ii-1, iia+mp-1 )
726  work( kk-iia+1 ) = work( kk-iia+1 ) +
727  $ abs( a( ioffa+kk ) )
728  680 CONTINUE
729  ioffa = ioffa + lda
730  690 CONTINUE
731  END IF
732  jj = jj + jb
733  END IF
734 *
735  IF( myrow.EQ.iarow )
736  $ ii = ii + jb
737  iarow = mod( iarow+1, nprow )
738  iacol = mod( iacol+1, npcol )
739 *
740  700 CONTINUE
741 *
742  ELSE
743 *
744 * Lower triangular matrix
745 *
746  ii = iia
747  jj = jja
748  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
749  jb = jn-ja+1
750 *
751  IF( mycol.EQ.iacol ) THEN
752  IF( myrow.EQ.iarow ) THEN
753  IF( udiag ) THEN
754  DO 720 ll = jj, jj + jb -1
755 * Unit diagonal entry
756  kk = ii+ll-jj
757  work( kk-iia+1 ) = work( kk-iia+1 ) + one
758  DO 710 kk = ii+ll-jj+1, iia+mp-1
759  work( kk-iia+1 ) = work( kk-iia+1 ) +
760  $ abs( a( ioffa+kk ) )
761  710 CONTINUE
762  ioffa = ioffa + lda
763  720 CONTINUE
764  ELSE
765  DO 740 ll = jj, jj + jb -1
766  DO 730 kk = ii+ll-jj, iia+mp-1
767  work( kk-iia+1 ) = work( kk-iia+1 ) +
768  $ abs( a( ioffa+kk ) )
769  730 CONTINUE
770  ioffa = ioffa + lda
771  740 CONTINUE
772  END IF
773  ELSE
774  DO 760 ll = jj, jj + jb -1
775  DO 750 kk = ii, iia+mp-1
776  work( kk-iia+1 ) = work( kk-iia+1 ) +
777  $ abs( a( ioffa+kk ) )
778  750 CONTINUE
779  ioffa = ioffa + lda
780  760 CONTINUE
781  END IF
782  jj = jj + jb
783  END IF
784 *
785  IF( myrow.EQ.iarow )
786  $ ii = ii + jb
787  iarow = mod( iarow+1, nprow )
788  iacol = mod( iacol+1, npcol )
789 *
790 * Loop over remaining block of columns
791 *
792  DO 830 j = jn+1, ja+n-1, desca( nb_ )
793  jb = min( ja+n-j, desca( nb_ ) )
794 *
795  IF( mycol.EQ.iacol ) THEN
796  IF( myrow.EQ.iarow ) THEN
797  IF( udiag ) THEN
798  DO 780 ll = jj, jj + jb -1
799 * Unit diagonal entry
800  kk = ii+ll-jj
801  work( kk-iia+1 ) = work( kk-iia+1 ) + one
802  DO 770 kk = ii+ll-jj+1, iia+mp-1
803  work( kk-iia+1 ) = work( kk-iia+1 ) +
804  $ abs( a( ioffa+kk ) )
805  770 CONTINUE
806  ioffa = ioffa + lda
807  780 CONTINUE
808  ELSE
809  DO 800 ll = jj, jj + jb -1
810  DO 790 kk = ii+ll-jj, iia+mp-1
811  work( kk-iia+1 ) = work( kk-iia+1 ) +
812  $ abs( a( ioffa+kk ) )
813  790 CONTINUE
814  ioffa = ioffa + lda
815  800 CONTINUE
816  END IF
817  ELSE
818  DO 820 ll = jj, jj + jb -1
819  DO 810 kk = ii, iia+mp-1
820  work( kk-iia+1 ) = work( kk-iia+1 ) +
821  $ abs( a( ioffa+kk ) )
822  810 CONTINUE
823  ioffa = ioffa + lda
824  820 CONTINUE
825  END IF
826  jj = jj + jb
827  END IF
828 *
829  IF( myrow.EQ.iarow )
830  $ ii = ii + jb
831  iarow = mod( iarow+1, nprow )
832  iacol = mod( iacol+1, npcol )
833 *
834  830 CONTINUE
835 *
836  END IF
837 *
838 * Find sum of global matrix rows and store on column 0 of
839 * process grid
840 *
841  CALL dgsum2d( ictxt, 'Rowwise', ' ', mp, 1, work, max( 1, mp ),
842  $ myrow, 0 )
843 *
844 * Find maximum sum of rows for Infinity-norm
845 *
846  IF( mycol.EQ.0 ) THEN
847  IF( mp.GT.0 ) THEN
848  VALUE = work( idamax( mp, work, 1 ) )
849  ELSE
850  VALUE = zero
851  END IF
852  CALL dgamx2d( ictxt, 'Columnwise', ' ', 1, 1, VALUE, 1, kk,
853  $ ll, -1, 0, 0 )
854  END IF
855 *
856 ************************************************************************
857 * Frobenius norm
858 * SSQ(1) is scale
859 * SSQ(2) is sum-of-squares
860 *
861  ELSE IF( lsame( norm, 'F' ) .OR. lsame( norm, 'E' ) ) THEN
862 *
863  IF( udiag ) THEN
864  ssq(1) = one
865  ssq(2) = dble( min( m, n ) ) / dble( nprow*npcol )
866  ELSE
867  ssq(1) = zero
868  ssq(2) = one
869  END IF
870 *
871  IF( lsame( uplo, 'U' ) ) THEN
872 *
873 * ***********************
874 * Upper triangular matrix
875 *
876  ii = iia
877  jj = jja
878  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
879  jb = jn-ja+1
880 *
881 * First block column of sub-matrix.
882 *
883  IF( mycol.EQ.iacol ) THEN
884  IF( myrow.EQ.iarow ) THEN
885 * This process has part of current block column,
886 * including diagonal block.
887  IF( udiag ) THEN
888  DO 840 ll = jj, jj + jb -1
889  colssq(1) = zero
890  colssq(2) = one
891  CALL zlassq( min( ii+ll-jj-1, iia+mp-1 )-iia+1,
892  $ a( iia+ioffa ), 1,
893  $ colssq(1), colssq(2) )
894  CALL dcombssq( ssq, colssq )
895  ioffa = ioffa + lda
896  840 CONTINUE
897  ELSE
898  DO 850 ll = jj, jj + jb -1
899  colssq(1) = zero
900  colssq(2) = one
901  CALL zlassq( min( ii+ll-jj, iia+mp-1 )-iia+1,
902  $ a( iia+ioffa ), 1,
903  $ colssq(1), colssq(2) )
904  CALL dcombssq( ssq, colssq )
905  ioffa = ioffa + lda
906  850 CONTINUE
907  END IF
908  ELSE
909 * This rank has part of current block column,
910 * but not diagonal block.
911 * It seems this lassq will be length 0, since ii = iia.
912  DO 860 ll = jj, jj + jb -1
913  colssq(1) = zero
914  colssq(2) = one
915  CALL zlassq( min( ii-1, iia+mp-1 )-iia+1,
916  $ a( iia+ioffa ), 1,
917  $ colssq(1), colssq(2) )
918  CALL dcombssq( ssq, colssq )
919  ioffa = ioffa + lda
920  860 CONTINUE
921  END IF
922  jj = jj + jb
923  END IF
924 *
925 * If this process has part of current block row, advance ii,
926 * then advance iarow, iacol to next diagonal block.
927 *
928  IF( myrow.EQ.iarow )
929  $ ii = ii + jb
930  iarow = mod( iarow+1, nprow )
931  iacol = mod( iacol+1, npcol )
932 *
933 * Loop over remaining block columns
934 *
935  DO 900 j = jn+1, ja+n-1, desca( nb_ )
936  jb = min( ja+n-j, desca( nb_ ) )
937 *
938  IF( mycol.EQ.iacol ) THEN
939  IF( myrow.EQ.iarow ) THEN
940  IF( udiag ) THEN
941  DO 870 ll = jj, jj + jb -1
942  colssq(1) = zero
943  colssq(2) = one
944  CALL zlassq( min(ii+ll-jj-1, iia+mp-1)-iia+1,
945  $ a( iia+ioffa ), 1,
946  $ colssq(1), colssq(2) )
947  CALL dcombssq( ssq, colssq )
948  ioffa = ioffa + lda
949  870 CONTINUE
950  ELSE
951  DO 880 ll = jj, jj + jb -1
952  colssq(1) = zero
953  colssq(2) = one
954  CALL zlassq( min( ii+ll-jj, iia+mp-1 )-iia+1,
955  $ a( iia+ioffa ), 1,
956  $ colssq(1), colssq(2) )
957  CALL dcombssq( ssq, colssq )
958  ioffa = ioffa + lda
959  880 CONTINUE
960  END IF
961  ELSE
962  DO 890 ll = jj, jj + jb -1
963  colssq(1) = zero
964  colssq(2) = one
965  CALL zlassq( min( ii-1, iia+mp-1 )-iia+1,
966  $ a( iia+ioffa ), 1,
967  $ colssq(1), colssq(2) )
968  CALL dcombssq( ssq, colssq )
969  ioffa = ioffa + lda
970  890 CONTINUE
971  END IF
972  jj = jj + jb
973  END IF
974 *
975  IF( myrow.EQ.iarow )
976  $ ii = ii + jb
977  iarow = mod( iarow+1, nprow )
978  iacol = mod( iacol+1, npcol )
979 *
980  900 CONTINUE
981 *
982  ELSE
983 *
984 * ***********************
985 * Lower triangular matrix
986 *
987  ii = iia
988  jj = jja
989  jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
990  jb = jn-ja+1
991 *
992  IF( mycol.EQ.iacol ) THEN
993  IF( myrow.EQ.iarow ) THEN
994  IF( udiag ) THEN
995  DO 910 ll = jj, jj + jb -1
996  colssq(1) = zero
997  colssq(2) = one
998  CALL zlassq( iia+mp-(ii+ll-jj+1),
999  $ a( ii+ll-jj+1+ioffa ), 1,
1000  $ colssq(1), colssq(2) )
1001  CALL dcombssq( ssq, colssq )
1002  ioffa = ioffa + lda
1003  910 CONTINUE
1004  ELSE
1005  DO 920 ll = jj, jj + jb -1
1006  colssq(1) = zero
1007  colssq(2) = one
1008  CALL zlassq( iia+mp-(ii+ll-jj),
1009  $ a( ii+ll-jj+ioffa ), 1,
1010  $ colssq(1), colssq(2) )
1011  CALL dcombssq( ssq, colssq )
1012  ioffa = ioffa + lda
1013  920 CONTINUE
1014  END IF
1015  ELSE
1016  DO 930 ll = jj, jj + jb -1
1017  colssq(1) = zero
1018  colssq(2) = one
1019  CALL zlassq( iia+mp-ii, a( ii+ioffa ), 1,
1020  $ colssq(1), colssq(2) )
1021  CALL dcombssq( ssq, colssq )
1022  ioffa = ioffa + lda
1023  930 CONTINUE
1024  END IF
1025  jj = jj + jb
1026  END IF
1027 *
1028  IF( myrow.EQ.iarow )
1029  $ ii = ii + jb
1030  iarow = mod( iarow+1, nprow )
1031  iacol = mod( iacol+1, npcol )
1032 *
1033 * Loop over remaining block of columns
1034 *
1035  DO 970 j = jn+1, ja+n-1, desca( nb_ )
1036  jb = min( ja+n-j, desca( nb_ ) )
1037 *
1038  IF( mycol.EQ.iacol ) THEN
1039  IF( myrow.EQ.iarow ) THEN
1040  IF( udiag ) THEN
1041  DO 940 ll = jj, jj + jb -1
1042  colssq(1) = zero
1043  colssq(2) = one
1044  CALL zlassq( iia+mp-(ii+ll-jj+1),
1045  $ a( ii+ll-jj+1+ioffa ), 1,
1046  $ colssq(1), colssq(2) )
1047  CALL dcombssq( ssq, colssq )
1048  ioffa = ioffa + lda
1049  940 CONTINUE
1050  ELSE
1051  DO 950 ll = jj, jj + jb -1
1052  colssq(1) = zero
1053  colssq(2) = one
1054  CALL zlassq( iia+mp-(ii+ll-jj),
1055  $ a( ii+ll-jj+ioffa ), 1,
1056  $ colssq(1), colssq(2) )
1057  CALL dcombssq( ssq, colssq )
1058  ioffa = ioffa + lda
1059  950 CONTINUE
1060  END IF
1061  ELSE
1062  DO 960 ll = jj, jj + jb -1
1063  colssq(1) = zero
1064  colssq(2) = one
1065  CALL zlassq( iia+mp-ii, a( ii+ioffa ), 1,
1066  $ colssq(1), colssq(2) )
1067  CALL dcombssq( ssq, colssq )
1068  ioffa = ioffa + lda
1069  960 CONTINUE
1070  END IF
1071  jj = jj + jb
1072  END IF
1073 *
1074  IF( myrow.EQ.iarow )
1075  $ ii = ii + jb
1076  iarow = mod( iarow+1, nprow )
1077  iacol = mod( iacol+1, npcol )
1078 *
1079  970 CONTINUE
1080 *
1081  END IF
1082 *
1083 * ***********************
1084 * Perform the global scaled sum
1085 *
1086  CALL pdtreecomb( ictxt, 'All', 2, ssq, 0, 0, dcombssq )
1087  VALUE = ssq( 1 ) * sqrt( ssq( 2 ) )
1088 *
1089  END IF
1090 *
1091 * Broadcast the result to every process in the grid.
1092 *
1093  IF( myrow.EQ.0 .AND. mycol.EQ.0 ) THEN
1094  CALL dgebs2d( ictxt, 'All', ' ', 1, 1, VALUE, 1 )
1095  ELSE
1096  CALL dgebr2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, 0, 0 )
1097  END IF
1098 *
1099  pzlantr = VALUE
1100 *
1101  RETURN
1102 *
1103 * End of PZLANTR
1104 *
1105  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
pzlantr
double precision function pzlantr(NORM, UPLO, DIAG, M, N, A, IA, JA, DESCA, WORK)
Definition: pzlantr.f:3
min
#define min(A, B)
Definition: pcgemr.c:181
iceil
integer function iceil(INUM, IDENOM)
Definition: iceil.f:2