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