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