SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
pzlansy.f
Go to the documentation of this file.
1 DOUBLE PRECISION FUNCTION pzlansy( NORM, UPLO, N, A, IA, JA,
2 $ DESCA, WORK )
3 IMPLICIT NONE
4*
5* -- ScaLAPACK auxiliary routine (version 1.7) --
6* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7* and University of California, Berkeley.
8* May 1, 1997
9*
10* .. Scalar Arguments ..
11 CHARACTER norm, uplo
12 INTEGER ia, ja, n
13* ..
14* .. Array Arguments ..
15 INTEGER desca( * )
16 DOUBLE PRECISION work( * )
17 COMPLEX*16 a( * )
18* ..
19*
20* Purpose
21* =======
22*
23* PZLANSY 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* real symmetric distributed matrix sub( A ) = A(IA:IA+N-1,JA:JA+N-1).
26*
27* PZLANSY returns the value
28*
29* ( max(abs(A(i,j))), NORM = 'M' or 'm' with IA <= i <= IA+N-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 PZLANSY as described
102* above.
103*
104* UPLO (global input) CHARACTER
105* Specifies whether the upper or lower triangular part of the
106* symmetric matrix sub( A ) is to be referenced.
107* = 'U': Upper triangular part of sub( A ) is referenced,
108* = 'L': Lower triangular part of sub( A ) is referenced.
109*
110* N (global input) INTEGER
111* The number of rows and columns to be operated on i.e the
112* number of rows and columns of the distributed submatrix
113* sub( A ). When N = 0, PZLANSY is set to zero. N >= 0.
114*
115* A (local input) COMPLEX*16 pointer into the local memory
116* to an array of dimension (LLD_A, LOCc(JA+N-1)) containing the
117* local pieces of the symmetric distributed matrix sub( A ).
118* If UPLO = 'U', the leading N-by-N upper triangular part of
119* sub( A ) contains the upper triangular matrix which norm is
120* to be computed, and the strictly lower triangular part of
121* this matrix is not referenced. If UPLO = 'L', the leading
122* N-by-N lower triangular part of sub( A ) contains the lower
123* triangular matrix which norm is to be computed, and the
124* strictly upper triangular part of sub( A ) is not referenced.
125*
126* IA (global input) INTEGER
127* The row index in the global array A indicating the first
128* row of sub( A ).
129*
130* JA (global input) INTEGER
131* The column index in the global array A indicating the
132* first column of sub( A ).
133*
134* DESCA (global and local input) INTEGER array of dimension DLEN_.
135* The array descriptor for the distributed matrix A.
136*
137* WORK (local workspace) DOUBLE PRECISION array dimension (LWORK)
138* LWORK >= 0 if NORM = 'M' or 'm' (not referenced),
139* 2*Nq0+Np0+LDW if NORM = '1', 'O', 'o', 'I' or 'i',
140* where LDW is given by:
141* IF( NPROW.NE.NPCOL ) THEN
142* LDW = MB_A*CEIL(CEIL(Np0/MB_A)/(LCM/NPROW))
143* ELSE
144* LDW = 0
145* END IF
146* 0 if NORM = 'F', 'f', 'E' or 'e' (not referenced),
147*
148* where LCM is the least common multiple of NPROW and NPCOL
149* LCM = ILCM( NPROW, NPCOL ) and CEIL denotes the ceiling
150* operation (ICEIL).
151*
152* IROFFA = MOD( IA-1, MB_A ), ICOFFA = MOD( JA-1, NB_A ),
153* IAROW = INDXG2P( IA, MB_A, MYROW, RSRC_A, NPROW ),
154* IACOL = INDXG2P( JA, NB_A, MYCOL, CSRC_A, NPCOL ),
155* Np0 = NUMROC( N+IROFFA, MB_A, MYROW, IAROW, NPROW ),
156* Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
157*
158* ICEIL, ILCM, INDXG2P and NUMROC are ScaLAPACK tool functions;
159* MYROW, MYCOL, NPROW and NPCOL can be determined by calling
160* the subroutine BLACS_GRIDINFO.
161*
162* =====================================================================
163*
164* .. Parameters ..
165 INTEGER block_cyclic_2d, csrc_, ctxt_, dlen_, dtype_,
166 $ lld_, mb_, m_, nb_, n_, rsrc_
167 parameter( block_cyclic_2d = 1, dlen_ = 9, dtype_ = 1,
168 $ ctxt_ = 2, m_ = 3, n_ = 4, mb_ = 5, nb_ = 6,
169 $ rsrc_ = 7, csrc_ = 8, lld_ = 9 )
170 DOUBLE PRECISION one, zero
171 parameter( one = 1.0d+0, zero = 0.0d+0 )
172* ..
173* .. Local Scalars ..
174 INTEGER i, iarow, iacol, ib, icoff, ictxt, icurcol,
175 $ icurrow, ii, iia, in, iroff, icsr, icsr0,
176 $ ioffa, irsc, irsc0, irsr, irsr0, jj, jja, k,
177 $ lda, ll, mycol, myrow, np, npcol, nprow, nq
178 DOUBLE PRECISION sum, value
179* ..
180* .. Local Arrays ..
181 DOUBLE PRECISION ssq( 2 ), colssq( 2 )
182* ..
183* .. External Subroutines ..
184 EXTERNAL blacs_gridinfo, daxpy, dcombssq,
185 $ dgamx2d, dgsum2d, dgebr2d,
186 $ dgebs2d, pdcol2row, pdtreecomb,
187 $ zlassq
188* ..
189* .. External Functions ..
190 LOGICAL lsame
191 INTEGER iceil, idamax, numroc
192 EXTERNAL iceil, idamax, lsame, numroc
193* ..
194* .. Intrinsic Functions ..
195 INTRINSIC abs, max, min, mod, sqrt
196* ..
197* .. Executable Statements ..
198*
199* Get grid parameters and local indexes.
200*
201 ictxt = desca( ctxt_ )
202 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
203 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol,
204 $ iia, jja, iarow, iacol )
205*
206 iroff = mod( ia-1, desca( mb_ ) )
207 icoff = mod( ja-1, desca( nb_ ) )
208 np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
209 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
210 icsr = 1
211 irsr = icsr + nq
212 irsc = irsr + nq
213 IF( myrow.EQ.iarow ) THEN
214 irsc0 = irsc + iroff
215 np = np - iroff
216 ELSE
217 irsc0 = irsc
218 END IF
219 IF( mycol.EQ.iacol ) THEN
220 icsr0 = icsr + icoff
221 irsr0 = irsr + icoff
222 nq = nq - icoff
223 ELSE
224 icsr0 = icsr
225 irsr0 = irsr
226 END IF
227 in = min( iceil( ia, desca( mb_ ) ) * desca( mb_ ), ia+n-1 )
228 lda = desca( lld_ )
229*
230* If the matrix is symmetric, we address only a triangular portion
231* of the matrix. A sum of row (column) i of the complete matrix
232* can be obtained by adding along row i and column i of the the
233* triangular matrix, stopping/starting at the diagonal, which is
234* the point of reflection. The pictures below demonstrate this.
235* In the following code, the row sums created by --- rows below are
236* refered to as ROWSUMS, and the column sums shown by | are refered
237* to as COLSUMS. Infinity-norm = 1-norm = ROWSUMS+COLSUMS.
238*
239* UPLO = 'U' UPLO = 'L'
240* ____i______ ___________
241* |\ | | |\ |
242* | \ | | | \ |
243* | \ | | | \ |
244* | \|------| i i|---\ |
245* | \ | | |\ |
246* | \ | | | \ |
247* | \ | | | \ |
248* | \ | | | \ |
249* | \ | | | \ |
250* | \ | | | \ |
251* |__________\| |___|______\|
252* i
253*
254* II, JJ : local indices into array A
255* ICURROW : process row containing diagonal block
256* ICURCOL : process column containing diagonal block
257* IRSC0 : pointer to part of work used to store the ROWSUMS while
258* they are stored along a process column
259* IRSR0 : pointer to part of work used to store the ROWSUMS after
260* they have been transposed to be along a process row
261*
262 ii = iia
263 jj = jja
264*
265 IF( n.EQ.0 ) THEN
266*
267 VALUE = zero
268*
269************************************************************************
270* max norm
271*
272 ELSE IF( lsame( norm, 'M' ) ) THEN
273*
274* Find max(abs(A(i,j))).
275*
276 VALUE = zero
277*
278 IF( lsame( uplo, 'U' ) ) THEN
279*
280* Handle first block separately
281*
282 ib = in-ia+1
283*
284* Find COLMAXS
285*
286 IF( mycol.EQ.iacol ) THEN
287 DO 20 k = (jj-1)*lda, (jj+ib-2)*lda, lda
288 IF( ii.GT.iia ) THEN
289 DO 10 ll = iia, ii-1
290 VALUE = max( VALUE, abs( a( ll+k ) ) )
291 10 CONTINUE
292 END IF
293 IF( myrow.EQ.iarow )
294 $ ii = ii + 1
295 20 CONTINUE
296*
297* Reset local indices so we can find ROWMAXS
298*
299 IF( myrow.EQ.iarow )
300 $ ii = ii - ib
301*
302 END IF
303*
304* Find ROWMAXS
305*
306 IF( myrow.EQ.iarow ) THEN
307 DO 40 k = ii, ii+ib-1
308 IF( jj.LE.jja+nq-1 ) THEN
309 DO 30 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
310 VALUE = max( VALUE, abs( a( k+ll ) ) )
311 30 CONTINUE
312 END IF
313 IF( mycol.EQ.iacol )
314 $ jj = jj + 1
315 40 CONTINUE
316 ii = ii + ib
317 ELSE IF( mycol.EQ.iacol ) THEN
318 jj = jj + ib
319 END IF
320*
321 icurrow = mod( iarow+1, nprow )
322 icurcol = mod( iacol+1, npcol )
323*
324* Loop over the remaining rows/columns of the matrix.
325*
326 DO 90 i = in+1, ia+n-1, desca( mb_ )
327 ib = min( desca( mb_ ), ia+n-i )
328*
329* Find COLMAXS
330*
331 IF( mycol.EQ.icurcol ) THEN
332 DO 60 k = (jj-1)*lda, (jj+ib-2)*lda, lda
333 IF( ii.GT.iia ) THEN
334 DO 50 ll = iia, ii-1
335 VALUE = max( VALUE, abs( a( ll+k ) ) )
336 50 CONTINUE
337 END IF
338 IF( myrow.EQ.icurrow )
339 $ ii = ii + 1
340 60 CONTINUE
341*
342* Reset local indices so we can find ROWMAXS
343*
344 IF( myrow.EQ.icurrow )
345 $ ii = ii - ib
346 END IF
347*
348* Find ROWMAXS
349*
350 IF( myrow.EQ.icurrow ) THEN
351 DO 80 k = ii, ii+ib-1
352 IF( jj.LE.jja+nq-1 ) THEN
353 DO 70 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
354 VALUE = max( VALUE, abs( a( k+ll ) ) )
355 70 CONTINUE
356 END IF
357 IF( mycol.EQ.icurcol )
358 $ jj = jj + 1
359 80 CONTINUE
360 ii = ii + ib
361 ELSE IF( mycol.EQ.icurcol ) THEN
362 jj = jj + ib
363 END IF
364 icurrow = mod( icurrow+1, nprow )
365 icurcol = mod( icurcol+1, npcol )
366 90 CONTINUE
367*
368 ELSE
369*
370* Handle first block separately
371*
372 ib = in-ia+1
373*
374* Find COLMAXS
375*
376 IF( mycol.EQ.iacol ) THEN
377 DO 110 k = (jj-1)*lda, (jj+ib-2)*lda, lda
378 IF( ii.LE.iia+np-1 ) THEN
379 DO 100 ll = ii, iia+np-1
380 VALUE = max( VALUE, abs( a( ll+k ) ) )
381 100 CONTINUE
382 END IF
383 IF( myrow.EQ.iarow )
384 $ ii = ii + 1
385 110 CONTINUE
386*
387* Reset local indices so we can find ROWMAXS
388*
389 IF( myrow.EQ.iarow )
390 $ ii = ii - ib
391 END IF
392*
393* Find ROWMAXS
394*
395 IF( myrow.EQ.iarow ) THEN
396 DO 130 k = 0, ib-1
397 IF( jj.GT.jja ) THEN
398 DO 120 ll = (jja-1)*lda, (jj-2)*lda, lda
399 VALUE = max( VALUE, abs( a( ii+ll ) ) )
400 120 CONTINUE
401 END IF
402 ii = ii + 1
403 IF( mycol.EQ.iacol )
404 $ jj = jj + 1
405 130 CONTINUE
406 ELSE IF( mycol.EQ.iacol ) THEN
407 jj = jj + ib
408 END IF
409*
410 icurrow = mod( iarow+1, nprow )
411 icurcol = mod( iacol+1, npcol )
412*
413* Loop over rows/columns of global matrix.
414*
415 DO 180 i = in+1, ia+n-1, desca( mb_ )
416 ib = min( desca( mb_ ), ia+n-i )
417*
418* Find COLMAXS
419*
420 IF( mycol.EQ.icurcol ) THEN
421 DO 150 k = (jj-1)*lda, (jj+ib-2)*lda, lda
422 IF( ii.LE.iia+np-1 ) THEN
423 DO 140 ll = ii, iia+np-1
424 VALUE = max( VALUE, abs( a( ll+k ) ) )
425 140 CONTINUE
426 END IF
427 IF( myrow.EQ.icurrow )
428 $ ii = ii + 1
429 150 CONTINUE
430*
431* Reset local indices so we can find ROWMAXS
432*
433 IF( myrow.EQ.icurrow )
434 $ ii = ii - ib
435 END IF
436*
437* Find ROWMAXS
438*
439 IF( myrow.EQ.icurrow ) THEN
440 DO 170 k = 0, ib-1
441 IF( jj.GT.jja ) THEN
442 DO 160 ll = (jja-1)*lda, (jj-2)*lda, lda
443 VALUE = max( VALUE, abs( a( ii+ll ) ) )
444 160 CONTINUE
445 END IF
446 ii = ii + 1
447 IF( mycol.EQ.icurcol )
448 $ jj = jj + 1
449 170 CONTINUE
450 ELSE IF( mycol.EQ.icurcol ) THEN
451 jj = jj + ib
452 END IF
453 icurrow = mod( icurrow+1, nprow )
454 icurcol = mod( icurcol+1, npcol )
455*
456 180 CONTINUE
457*
458 END IF
459*
460* Gather the result on process (IAROW,IACOL).
461*
462 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, i, k, -1,
463 $ iarow, iacol )
464*
465************************************************************************
466* one or inf norm
467*
468 ELSE IF( lsame( norm, 'I' ) .OR. lsame( norm, 'O' ) .OR.
469 $ norm.EQ.'1' ) THEN
470*
471* Find normI( sub( A ) ) ( = norm1( sub( A ) ), since sub( A ) is
472* symmetric).
473*
474 IF( lsame( uplo, 'U' ) ) THEN
475*
476* Handle first block separately
477*
478 ib = in-ia+1
479*
480* Find COLSUMS
481*
482 IF( mycol.EQ.iacol ) THEN
483 ioffa = ( jj - 1 ) * lda
484 DO 200 k = 0, ib-1
485 sum = zero
486 IF( ii.GT.iia ) THEN
487 DO 190 ll = iia, ii-1
488 sum = sum + abs( a( ll+ioffa ) )
489 190 CONTINUE
490 END IF
491 ioffa = ioffa + lda
492 work( jj+k-jja+icsr0 ) = sum
493 IF( myrow.EQ.iarow )
494 $ ii = ii + 1
495 200 CONTINUE
496*
497* Reset local indices so we can find ROWSUMS
498*
499 IF( myrow.EQ.iarow )
500 $ ii = ii - ib
501*
502 END IF
503*
504* Find ROWSUMS
505*
506 IF( myrow.EQ.iarow ) THEN
507 DO 220 k = ii, ii+ib-1
508 sum = zero
509 IF( jja+nq.GT.jj ) THEN
510 DO 210 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
511 sum = sum + abs( a( k+ll ) )
512 210 CONTINUE
513 END IF
514 work( k-iia+irsc0 ) = sum
515 IF( mycol.EQ.iacol )
516 $ jj = jj + 1
517 220 CONTINUE
518 ii = ii + ib
519 ELSE IF( mycol.EQ.iacol ) THEN
520 jj = jj + ib
521 END IF
522*
523 icurrow = mod( iarow+1, nprow )
524 icurcol = mod( iacol+1, npcol )
525*
526* Loop over remaining rows/columns of global matrix.
527*
528 DO 270 i = in+1, ia+n-1, desca( mb_ )
529 ib = min( desca( mb_ ), ia+n-i )
530*
531* Find COLSUMS
532*
533 IF( mycol.EQ.icurcol ) THEN
534 ioffa = ( jj - 1 ) * lda
535 DO 240 k = 0, ib-1
536 sum = zero
537 IF( ii.GT.iia ) THEN
538 DO 230 ll = iia, ii-1
539 sum = sum + abs( a( ioffa+ll ) )
540 230 CONTINUE
541 END IF
542 ioffa = ioffa + lda
543 work( jj+k-jja+icsr0 ) = sum
544 IF( myrow.EQ.icurrow )
545 $ ii = ii + 1
546 240 CONTINUE
547*
548* Reset local indices so we can find ROWSUMS
549*
550 IF( myrow.EQ.icurrow )
551 $ ii = ii - ib
552*
553 END IF
554*
555* Find ROWSUMS
556*
557 IF( myrow.EQ.icurrow ) THEN
558 DO 260 k = ii, ii+ib-1
559 sum = zero
560 IF( jja+nq.GT.jj ) THEN
561 DO 250 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
562 sum = sum + abs( a( k+ll ) )
563 250 CONTINUE
564 END IF
565 work( k-iia+irsc0 ) = sum
566 IF( mycol.EQ.icurcol )
567 $ jj = jj + 1
568 260 CONTINUE
569 ii = ii + ib
570 ELSE IF( mycol.EQ.icurcol ) THEN
571 jj = jj + ib
572 END IF
573*
574 icurrow = mod( icurrow+1, nprow )
575 icurcol = mod( icurcol+1, npcol )
576*
577 270 CONTINUE
578*
579 ELSE
580*
581* Handle first block separately
582*
583 ib = in-ia+1
584*
585* Find COLSUMS
586*
587 IF( mycol.EQ.iacol ) THEN
588 ioffa = (jj-1)*lda
589 DO 290 k = 0, ib-1
590 sum = zero
591 IF( iia+np.GT.ii ) THEN
592 DO 280 ll = ii, iia+np-1
593 sum = sum + abs( a( ioffa+ll ) )
594 280 CONTINUE
595 END IF
596 ioffa = ioffa + lda
597 work( jj+k-jja+icsr0 ) = sum
598 IF( myrow.EQ.iarow )
599 $ ii = ii + 1
600 290 CONTINUE
601*
602* Reset local indices so we can find ROWSUMS
603*
604 IF( myrow.EQ.iarow )
605 $ ii = ii - ib
606*
607 END IF
608*
609* Find ROWSUMS
610*
611 IF( myrow.EQ.iarow ) THEN
612 DO 310 k = ii, ii+ib-1
613 sum = zero
614 IF( jj.GT.jja ) THEN
615 DO 300 ll = (jja-1)*lda, (jj-2)*lda, lda
616 sum = sum + abs( a( k+ll ) )
617 300 CONTINUE
618 END IF
619 work( k-iia+irsc0 ) = sum
620 IF( mycol.EQ.iacol )
621 $ jj = jj + 1
622 310 CONTINUE
623 ii = ii + ib
624 ELSE IF( mycol.EQ.iacol ) THEN
625 jj = jj + ib
626 END IF
627*
628 icurrow = mod( iarow+1, nprow )
629 icurcol = mod( iacol+1, npcol )
630*
631* Loop over rows/columns of global matrix.
632*
633 DO 360 i = in+1, ia+n-1, desca( mb_ )
634 ib = min( desca( mb_ ), ia+n-i )
635*
636* Find COLSUMS
637*
638 IF( mycol.EQ.icurcol ) THEN
639 ioffa = ( jj - 1 ) * lda
640 DO 330 k = 0, ib-1
641 sum = zero
642 IF( iia+np.GT.ii ) THEN
643 DO 320 ll = ii, iia+np-1
644 sum = sum + abs( a( ll+ioffa ) )
645 320 CONTINUE
646 END IF
647 ioffa = ioffa + lda
648 work( jj+k-jja+icsr0 ) = sum
649 IF( myrow.EQ.icurrow )
650 $ ii = ii + 1
651 330 CONTINUE
652*
653* Reset local indices so we can find ROWSUMS
654*
655 IF( myrow.EQ.icurrow )
656 $ ii = ii - ib
657*
658 END IF
659*
660* Find ROWSUMS
661*
662 IF( myrow.EQ.icurrow ) THEN
663 DO 350 k = ii, ii+ib-1
664 sum = zero
665 IF( jj.GT.jja ) THEN
666 DO 340 ll = (jja-1)*lda, (jj-2)*lda, lda
667 sum = sum + abs( a( k+ll ) )
668 340 CONTINUE
669 END IF
670 work(k-iia+irsc0) = sum
671 IF( mycol.EQ.icurcol )
672 $ jj = jj + 1
673 350 CONTINUE
674 ii = ii + ib
675 ELSE IF( mycol.EQ.icurcol ) THEN
676 jj = jj + ib
677 END IF
678*
679 icurrow = mod( icurrow+1, nprow )
680 icurcol = mod( icurcol+1, npcol )
681*
682 360 CONTINUE
683 END IF
684*
685* After calls to DGSUM2D, process row 0 will have global
686* COLSUMS and process column 0 will have global ROWSUMS.
687* Transpose ROWSUMS and add to COLSUMS to get global row/column
688* sum, the max of which is the infinity or 1 norm.
689*
690 IF( mycol.EQ.iacol )
691 $ nq = nq + icoff
692 CALL dgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work( icsr ), 1,
693 $ iarow, mycol )
694 IF( myrow.EQ.iarow )
695 $ np = np + iroff
696 CALL dgsum2d( ictxt, 'Rowwise', ' ', np, 1, work( irsc ),
697 $ max( 1, np ), myrow, iacol )
698*
699 CALL pdcol2row( ictxt, n, 1, desca( mb_ ), work( irsc ),
700 $ max( 1, np ), work( irsr ), max( 1, nq ),
701 $ iarow, iacol, iarow, iacol, work( irsc+np ) )
702*
703 IF( myrow.EQ.iarow ) THEN
704 IF( mycol.EQ.iacol )
705 $ nq = nq - icoff
706 CALL daxpy( nq, one, work( irsr0 ), 1, work( icsr0 ), 1 )
707 IF( nq.LT.1 ) THEN
708 VALUE = zero
709 ELSE
710 VALUE = work( idamax( nq, work( icsr0 ), 1 ) )
711 END IF
712 CALL dgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, i, k,
713 $ -1, iarow, iacol )
714 END IF
715*
716************************************************************************
717* Frobenius norm
718* SSQ(1) is scale
719* SSQ(2) is sum-of-squares
720*
721 ELSE IF( lsame( norm, 'F' ) .OR. lsame( norm, 'E' ) ) THEN
722*
723* Find normF( sub( A ) ).
724*
725 ssq(1) = zero
726 ssq(2) = one
727*
728* Add off-diagonal entries, first
729*
730 IF( lsame( uplo, 'U' ) ) THEN
731*
732* Handle first block separately
733*
734 ib = in-ia+1
735*
736 IF( mycol.EQ.iacol ) THEN
737 DO 370 k = (jj-1)*lda, (jj+ib-2)*lda, lda
738 colssq(1) = zero
739 colssq(2) = one
740 CALL zlassq( ii-iia, a( iia+k ), 1,
741 $ colssq(1), colssq(2) )
742 IF( myrow.EQ.iarow )
743 $ ii = ii + 1
744 CALL zlassq( ii-iia, a( iia+k ), 1,
745 $ colssq(1), colssq(2) )
746 CALL dcombssq( ssq, colssq )
747 370 CONTINUE
748*
749 jj = jj + ib
750 ELSE IF( myrow.EQ.iarow ) THEN
751 ii = ii + ib
752 END IF
753*
754 icurrow = mod( iarow+1, nprow )
755 icurcol = mod( iacol+1, npcol )
756*
757* Loop over rows/columns of global matrix.
758*
759 DO 390 i = in+1, ia+n-1, desca( mb_ )
760 ib = min( desca( mb_ ), ia+n-i )
761*
762 IF( mycol.EQ.icurcol ) THEN
763 DO 380 k = (jj-1)*lda, (jj+ib-2)*lda, lda
764 colssq(1) = zero
765 colssq(2) = one
766 CALL zlassq( ii-iia, a( iia+k ), 1,
767 $ colssq(1), colssq(2) )
768 IF( myrow.EQ.icurrow )
769 $ ii = ii + 1
770 CALL zlassq( ii-iia, a(iia+k ), 1,
771 $ colssq(1), colssq(2) )
772 CALL dcombssq( ssq, colssq )
773 380 CONTINUE
774*
775 jj = jj + ib
776 ELSE IF( myrow.EQ.icurrow ) THEN
777 ii = ii + ib
778 END IF
779*
780 icurrow = mod( icurrow+1, nprow )
781 icurcol = mod( icurcol+1, npcol )
782*
783 390 CONTINUE
784*
785 ELSE
786*
787* Handle first block separately
788*
789 ib = in-ia+1
790*
791 IF( mycol.EQ.iacol ) THEN
792 DO 400 k = (jj-1)*lda, (jj+ib-2)*lda, lda
793 colssq(1) = zero
794 colssq(2) = one
795 CALL zlassq( iia+np-ii, a( ii+k ), 1,
796 $ colssq(1), colssq(2) )
797 IF( myrow.EQ.iarow )
798 $ ii = ii + 1
799 CALL zlassq( iia+np-ii, a( ii+k ), 1,
800 $ colssq(1), colssq(2) )
801 CALL dcombssq( ssq, colssq )
802 400 CONTINUE
803*
804 jj = jj + ib
805 ELSE IF( myrow.EQ.iarow ) THEN
806 ii = ii + ib
807 END IF
808*
809 icurrow = mod( iarow+1, nprow )
810 icurcol = mod( iacol+1, npcol )
811*
812* Loop over rows/columns of global matrix.
813*
814 DO 420 i = in+1, ia+n-1, desca( mb_ )
815 ib = min( desca( mb_ ), ia+n-i )
816*
817 IF( mycol.EQ.icurcol ) THEN
818 DO 410 k = (jj-1)*lda, (jj+ib-2)*lda, lda
819 colssq(1) = zero
820 colssq(2) = one
821 CALL zlassq( iia+np-ii, a( ii+k ), 1,
822 $ colssq(1), colssq(2) )
823 IF( myrow.EQ.icurrow )
824 $ ii = ii + 1
825 CALL zlassq( iia+np-ii, a( ii+k ), 1,
826 $ colssq(1), colssq(2) )
827 CALL dcombssq( ssq, colssq )
828 410 CONTINUE
829*
830 jj = jj + ib
831 ELSE IF( myrow.EQ.icurrow ) THEN
832 ii = ii + ib
833 END IF
834*
835 icurrow = mod( icurrow+1, nprow )
836 icurcol = mod( icurcol+1, npcol )
837*
838 420 CONTINUE
839*
840 END IF
841*
842* Perform the global scaled sum
843*
844 CALL pdtreecomb( ictxt, 'All', 2, ssq, iarow, iacol,
845 $ dcombssq )
846 VALUE = ssq( 1 ) * sqrt( ssq( 2 ) )
847*
848 END IF
849*
850* Broadcast the result to the other processes
851*
852 IF( myrow.EQ.iarow .AND. mycol.EQ.iacol ) THEN
853 CALL dgebs2d( ictxt, 'All', ' ', 1, 1, VALUE, 1 )
854 ELSE
855 CALL dgebr2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, iarow,
856 $ iacol )
857 END IF
858*
859 pzlansy = VALUE
860*
861 RETURN
862*
863* End of PZLANSY
864*
865 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
pdcol2row
subroutine pdcol2row(ictxt, m, n, nb, vs, ldvs, vd, ldvd, rsrc, csrc, rdest, cdest, work)
Definition pdcol2row.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
pzlansy
double precision function pzlansy(norm, uplo, n, a, ia, ja, desca, work)
Definition pzlansy.f:3
lsame
logical function lsame(ca, cb)
Definition tools.f:1724