SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
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
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
pdlantr
double precision function pdlantr(norm, uplo, diag, m, n, a, ia, ja, desca, work)
Definition pdlantr.f:3
dcombssq
subroutine dcombssq(v1, v2)
Definition pdtreecomb.f:259
pdtreecomb
subroutine pdtreecomb(ictxt, scope, n, mine, rdest0, cdest0, subptr)
Definition pdtreecomb.f:3
lsame
logical function lsame(ca, cb)
Definition tools.f:1724