SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ pzlantr()

double precision function pzlantr ( character  norm,
character  uplo,
character  diag,
integer  m,
integer  n,
complex*16, dimension( * )  a,
integer  ia,
integer  ja,
integer, dimension( * )  desca,
double precision, dimension( * )  work 
)

Definition at line 1 of file pzlantr.f.

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*
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition infog2l.f:3
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition numroc.f:2
max
#define max(A, B)
Definition pcgemr.c:180
min
#define min(A, B)
Definition pcgemr.c:181
dcombssq
subroutine dcombssq(v1, v2)
Definition pdtreecomb.f:259
pdtreecomb
subroutine pdtreecomb(ictxt, scope, n, mine, rdest0, cdest0, subptr)
Definition pdtreecomb.f:3
pzlantr
double precision function pzlantr(norm, uplo, diag, m, n, a, ia, ja, desca, work)
Definition pzlantr.f:3
lsame
logical function lsame(ca, cb)
Definition tools.f:1724
Here is the call graph for this function:
Here is the caller graph for this function: