OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pdlansy.f File Reference

Go to the source code of this file.

Functions/Subroutines

double precision function pdlansy (norm, uplo, n, a, ia, ja, desca, work)

Function/Subroutine Documentation

◆ pdlansy()

double precision function pdlansy ( character norm,
character uplo,
integer n,
double precision, dimension( * ) a,
integer ia,
integer ja,
integer, dimension( * ) desca,
double precision, dimension( * ) work )

Definition at line 1 of file pdlansy.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 NORM, UPLO
12 INTEGER IA, JA, N
13* ..
14* .. Array Arguments ..
15 INTEGER DESCA( * )
16 DOUBLE PRECISION A( * ), WORK( * )
17* ..
18*
19* Purpose
20* =======
21*
22* PDLANSY 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* PDLANSY 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 PDLANSY 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, PDLANSY is set to zero. N >= 0.
113*
114* A (local input) DOUBLE PRECISION 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) DOUBLE PRECISION 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 DOUBLE PRECISION ONE, ZERO
170 parameter( one = 1.0d+0, zero = 0.0d+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 DOUBLE PRECISION SUM, VALUE
178* ..
179* .. Local Arrays ..
180 DOUBLE PRECISION SSQ( 2 ), COLSSQ( 2 )
181* ..
182* .. External Subroutines ..
183 EXTERNAL blacs_gridinfo, daxpy, dcombssq,
184 $ dgamx2d, dgsum2d, dgebr2d,
185 $ dgebs2d, dlassq, pdcol2row,
186 $ pdtreecomb
187* ..
188* .. External Functions ..
189 LOGICAL LSAME
190 INTEGER ICEIL, IDAMAX, NUMROC
191 EXTERNAL iceil, idamax, lsame, numroc
192* ..
193* .. Intrinsic Functions ..
194 INTRINSIC abs, max, min, mod, sqrt
195* ..
196* .. Executable Statements ..
197*
198* Get grid parameters and local indexes.
199*
200 ictxt = desca( ctxt_ )
201 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
202 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol,
203 $ iia, jja, iarow, iacol )
204*
205 iroff = mod( ia-1, desca( mb_ ) )
206 icoff = mod( ja-1, desca( nb_ ) )
207 np = numroc( n+iroff, desca( mb_ ), myrow, iarow, nprow )
208 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
209 icsr = 1
210 irsr = icsr + nq
211 irsc = irsr + nq
212 IF( myrow.EQ.iarow ) THEN
213 irsc0 = irsc + iroff
214 np = np - iroff
215 ELSE
216 irsc0 = irsc
217 END IF
218 IF( mycol.EQ.iacol ) THEN
219 icsr0 = icsr + icoff
220 irsr0 = irsr + icoff
221 nq = nq - icoff
222 ELSE
223 icsr0 = icsr
224 irsr0 = irsr
225 END IF
226 in = min( iceil( ia, desca( mb_ ) ) * desca( mb_ ), ia+n-1 )
227 lda = desca( lld_ )
228*
229* If the matrix is symmetric, we address only a triangular portion
230* of the matrix. A sum of row (column) i of the complete matrix
231* can be obtained by adding along row i and column i of the the
232* triangular matrix, stopping/starting at the diagonal, which is
233* the point of reflection. The pictures below demonstrate this.
234* In the following code, the row sums created by --- rows below are
235* refered to as ROWSUMS, and the column sums shown by | are refered
236* to as COLSUMS. Infinity-norm = 1-norm = ROWSUMS+COLSUMS.
237*
238* UPLO = 'U' UPLO = 'L'
239* ____i______ ___________
240* |\ | | |\ |
241* | \ | | | \ |
242* | \ | | | \ |
243* | \|------| i i|---\ |
244* | \ | | |\ |
245* | \ | | | \ |
246* | \ | | | \ |
247* | \ | | | \ |
248* | \ | | | \ |
249* | \ | | | \ |
250* |__________\| |___|______\|
251* i
252*
253* II, JJ : local indices into array A
254* ICURROW : process row containing diagonal block
255* ICURCOL : process column containing diagonal block
256* IRSC0 : pointer to part of work used to store the ROWSUMS while
257* they are stored along a process column
258* IRSR0 : pointer to part of work used to store the ROWSUMS after
259* they have been transposed to be along a process row
260*
261 ii = iia
262 jj = jja
263*
264 IF( n.EQ.0 ) THEN
265*
266 VALUE = zero
267*
268************************************************************************
269* max norm
270*
271 ELSE IF( lsame( norm, 'M' ) ) THEN
272*
273* Find max(abs(A(i,j))).
274*
275 VALUE = zero
276*
277 IF( lsame( uplo, 'U' ) ) THEN
278*
279* Handle first block separately
280*
281 ib = in-ia+1
282*
283* Find COLMAXS
284*
285 IF( mycol.EQ.iacol ) THEN
286 DO 20 k = (jj-1)*lda, (jj+ib-2)*lda, lda
287 IF( ii.GT.iia ) THEN
288 DO 10 ll = iia, ii-1
289 VALUE = max( VALUE, abs( a( ll+k ) ) )
290 10 CONTINUE
291 END IF
292 IF( myrow.EQ.iarow )
293 $ ii = ii + 1
294 20 CONTINUE
295*
296* Reset local indices so we can find ROWMAXS
297*
298 IF( myrow.EQ.iarow )
299 $ ii = ii - ib
300*
301 END IF
302*
303* Find ROWMAXS
304*
305 IF( myrow.EQ.iarow ) THEN
306 DO 40 k = ii, ii+ib-1
307 IF( jj.LE.jja+nq-1 ) THEN
308 DO 30 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
309 VALUE = max( VALUE, abs( a( k+ll ) ) )
310 30 CONTINUE
311 END IF
312 IF( mycol.EQ.iacol )
313 $ jj = jj + 1
314 40 CONTINUE
315 ii = ii + ib
316 ELSE IF( mycol.EQ.iacol ) THEN
317 jj = jj + ib
318 END IF
319*
320 icurrow = mod( iarow+1, nprow )
321 icurcol = mod( iacol+1, npcol )
322*
323* Loop over the remaining rows/columns of the matrix.
324*
325 DO 90 i = in+1, ia+n-1, desca( mb_ )
326 ib = min( desca( mb_ ), ia+n-i )
327*
328* Find COLMAXS
329*
330 IF( mycol.EQ.icurcol ) THEN
331 DO 60 k = (jj-1)*lda, (jj+ib-2)*lda, lda
332 IF( ii.GT.iia ) THEN
333 DO 50 ll = iia, ii-1
334 VALUE = max( VALUE, abs( a( ll+k ) ) )
335 50 CONTINUE
336 END IF
337 IF( myrow.EQ.icurrow )
338 $ ii = ii + 1
339 60 CONTINUE
340*
341* Reset local indices so we can find ROWMAXS
342*
343 IF( myrow.EQ.icurrow )
344 $ ii = ii - ib
345 END IF
346*
347* Find ROWMAXS
348*
349 IF( myrow.EQ.icurrow ) THEN
350 DO 80 k = ii, ii+ib-1
351 IF( jj.LE.jja+nq-1 ) THEN
352 DO 70 ll = (jj-1)*lda, (jja+nq-2)*lda, lda
353 VALUE = max( VALUE, abs( a( k+ll ) ) )
354 70 CONTINUE
355 END IF
356 IF( mycol.EQ.icurcol )
357 $ jj = jj + 1
358 80 CONTINUE
359 ii = ii + ib
360 ELSE IF( mycol.EQ.icurcol ) THEN
361 jj = jj + ib
362 END IF
363 icurrow = mod( icurrow+1, nprow )
364 icurcol = mod( icurcol+1, npcol )
365 90 CONTINUE
366*
367 ELSE
368*
369* Handle first block separately
370*
371 ib = in-ia+1
372*
373* Find COLMAXS
374*
375 IF( mycol.EQ.iacol ) THEN
376 DO 110 k = (jj-1)*lda, (jj+ib-2)*lda, lda
377 IF( ii.LE.iia+np-1 ) THEN
378 DO 100 ll = ii, iia+np-1
379 VALUE = max( VALUE, abs( a( ll+k ) ) )
380 100 CONTINUE
381 END IF
382 IF( myrow.EQ.iarow )
383 $ ii = ii + 1
384 110 CONTINUE
385*
386* Reset local indices so we can find ROWMAXS
387*
388 IF( myrow.EQ.iarow )
389 $ ii = ii - ib
390 END IF
391*
392* Find ROWMAXS
393*
394 IF( myrow.EQ.iarow ) THEN
395 DO 130 k = 0, ib-1
396 IF( jj.GT.jja ) THEN
397 DO 120 ll = (jja-1)*lda, (jj-2)*lda, lda
398 VALUE = max( VALUE, abs( a( ii+ll ) ) )
399 120 CONTINUE
400 END IF
401 ii = ii + 1
402 IF( mycol.EQ.iacol )
403 $ jj = jj + 1
404 130 CONTINUE
405 ELSE IF( mycol.EQ.iacol ) THEN
406 jj = jj + ib
407 END IF
408*
409 icurrow = mod( iarow+1, nprow )
410 icurcol = mod( iacol+1, npcol )
411*
412* Loop over rows/columns of global matrix.
413*
414 DO 180 i = in+1, ia+n-1, desca( mb_ )
415 ib = min( desca( mb_ ), ia+n-i )
416*
417* Find COLMAXS
418*
419 IF( mycol.EQ.icurcol ) THEN
420 DO 150 k = (jj-1)*lda, (jj+ib-2)*lda, lda
421 IF( ii.LE.iia+np-1 ) THEN
422 DO 140 ll = ii, iia+np-1
423 VALUE = max( VALUE, abs( a( ll+k ) ) )
424 140 CONTINUE
425 END IF
426 IF( myrow.EQ.icurrow )
427 $ ii = ii + 1
428 150 CONTINUE
429*
430* Reset local indices so we can find ROWMAXS
431*
432 IF( myrow.EQ.icurrow )
433 $ ii = ii - ib
434 END IF
435*
436* Find ROWMAXS
437*
438 IF( myrow.EQ.icurrow ) THEN
439 DO 170 k = 0, ib-1
440 IF( jj.GT.jja ) THEN
441 DO 160 ll = (jja-1)*lda, (jj-2)*lda, lda
442 VALUE = max( VALUE, abs( a( ii+ll ) ) )
443 160 CONTINUE
444 END IF
445 ii = ii + 1
446 IF( mycol.EQ.icurcol )
447 $ jj = jj + 1
448 170 CONTINUE
449 ELSE IF( mycol.EQ.icurcol ) THEN
450 jj = jj + ib
451 END IF
452 icurrow = mod( icurrow+1, nprow )
453 icurcol = mod( icurcol+1, npcol )
454*
455 180 CONTINUE
456*
457 END IF
458*
459* Gather the result on process (IAROW,IACOL).
460*
461 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, i, k, -1,
462 $ iarow, iacol )
463*
464************************************************************************
465* one or inf norm
466*
467 ELSE IF( lsame( norm, 'I' ) .OR. lsame( norm, 'O' ) .OR.
468 $ norm.EQ.'1' ) THEN
469*
470* Find normI( sub( A ) ) ( = norm1( sub( A ) ), since sub( A ) is
471* symmetric).
472*
473 IF( LSAME( UPLO, 'u' ) ) THEN
474*
475* Handle first block separately
476*
477 IB = IN-IA+1
478*
479* Find COLSUMS
480*
481.EQ. IF( MYCOLIACOL ) THEN
482 IOFFA = ( JJ - 1 ) * LDA
483 DO 200 K = 0, IB-1
484 SUM = ZERO
485.GT. IF( IIIIA ) THEN
486 DO 190 LL = IIA, II-1
487 SUM = SUM + ABS( A( LL+IOFFA ) )
488 190 CONTINUE
489 END IF
490 IOFFA = IOFFA + LDA
491 WORK( JJ+K-JJA+ICSR0 ) = SUM
492.EQ. IF( MYROWIAROW )
493 $ II = II + 1
494 200 CONTINUE
495*
496* Reset local indices so we can find ROWSUMS
497*
498.EQ. IF( MYROWIAROW )
499 $ II = II - IB
500*
501 END IF
502*
503* Find ROWSUMS
504*
505.EQ. IF( MYROWIAROW ) THEN
506 DO 220 K = II, II+IB-1
507 SUM = ZERO
508.GT. IF( JJA+NQJJ ) THEN
509 DO 210 LL = (JJ-1)*LDA, (JJA+NQ-2)*LDA, LDA
510 SUM = SUM + ABS( A( K+LL ) )
511 210 CONTINUE
512 END IF
513 WORK( K-IIA+IRSC0 ) = SUM
514.EQ. IF( MYCOLIACOL )
515 $ JJ = JJ + 1
516 220 CONTINUE
517 II = II + IB
518.EQ. ELSE IF( MYCOLIACOL ) THEN
519 JJ = JJ + IB
520 END IF
521*
522 ICURROW = MOD( IAROW+1, NPROW )
523 ICURCOL = MOD( IACOL+1, NPCOL )
524*
525* Loop over remaining rows/columns of global matrix.
526*
527 DO 270 I = IN+1, IA+N-1, DESCA( MB_ )
528 IB = MIN( DESCA( MB_ ), IA+N-I )
529*
530* Find COLSUMS
531*
532.EQ. IF( MYCOLICURCOL ) THEN
533 IOFFA = ( JJ - 1 ) * LDA
534 DO 240 K = 0, IB-1
535 SUM = ZERO
536.GT. IF( IIIIA ) THEN
537 DO 230 LL = IIA, II-1
538 SUM = SUM + ABS( A( IOFFA+LL ) )
539 230 CONTINUE
540 END IF
541 IOFFA = IOFFA + LDA
542 WORK( JJ+K-JJA+ICSR0 ) = SUM
543.EQ. IF( MYROWICURROW )
544 $ II = II + 1
545 240 CONTINUE
546*
547* Reset local indices so we can find ROWSUMS
548*
549.EQ. IF( MYROWICURROW )
550 $ II = II - IB
551*
552 END IF
553*
554* Find ROWSUMS
555*
556.EQ. IF( MYROWICURROW ) THEN
557 DO 260 K = II, II+IB-1
558 SUM = ZERO
559.GT. IF( JJA+NQJJ ) THEN
560 DO 250 LL = (JJ-1)*LDA, (JJA+NQ-2)*LDA, LDA
561 SUM = SUM + ABS( A( K+LL ) )
562 250 CONTINUE
563 END IF
564 WORK( K-IIA+IRSC0 ) = SUM
565.EQ. IF( MYCOLICURCOL )
566 $ JJ = JJ + 1
567 260 CONTINUE
568 II = II + IB
569.EQ. ELSE IF( MYCOLICURCOL ) THEN
570 JJ = JJ + IB
571 END IF
572*
573 ICURROW = MOD( ICURROW+1, NPROW )
574 ICURCOL = MOD( ICURCOL+1, NPCOL )
575*
576 270 CONTINUE
577*
578 ELSE
579*
580* Handle first block separately
581*
582 IB = IN-IA+1
583*
584* Find COLSUMS
585*
586.EQ. IF( MYCOLIACOL ) THEN
587 IOFFA = (JJ-1)*LDA
588 DO 290 K = 0, IB-1
589 SUM = ZERO
590.GT. IF( IIA+NPII ) THEN
591 DO 280 LL = II, IIA+NP-1
592 SUM = SUM + ABS( A( IOFFA+LL ) )
593 280 CONTINUE
594 END IF
595 IOFFA = IOFFA + LDA
596 WORK( JJ+K-JJA+ICSR0 ) = SUM
597.EQ. IF( MYROWIAROW )
598 $ II = II + 1
599 290 CONTINUE
600*
601* Reset local indices so we can find ROWSUMS
602*
603.EQ. IF( MYROWIAROW )
604 $ II = II - IB
605*
606 END IF
607*
608* Find ROWSUMS
609*
610.EQ. IF( MYROWIAROW ) THEN
611 DO 310 K = II, II+IB-1
612 SUM = ZERO
613.GT. IF( JJJJA ) THEN
614 DO 300 LL = (JJA-1)*LDA, (JJ-2)*LDA, LDA
615 SUM = SUM + ABS( A( K+LL ) )
616 300 CONTINUE
617 END IF
618 WORK( K-IIA+IRSC0 ) = SUM
619.EQ. IF( MYCOLIACOL )
620 $ JJ = JJ + 1
621 310 CONTINUE
622 II = II + IB
623.EQ. ELSE IF( MYCOLIACOL ) THEN
624 JJ = JJ + IB
625 END IF
626*
627 ICURROW = MOD( IAROW+1, NPROW )
628 ICURCOL = MOD( IACOL+1, NPCOL )
629*
630* Loop over rows/columns of global matrix.
631*
632 DO 360 I = IN+1, IA+N-1, DESCA( MB_ )
633 IB = MIN( DESCA( MB_ ), IA+N-I )
634*
635* Find COLSUMS
636*
637.EQ. IF( MYCOLICURCOL ) THEN
638 IOFFA = ( JJ - 1 ) * LDA
639 DO 330 K = 0, IB-1
640 SUM = ZERO
641.GT. IF( IIA+NPII ) THEN
642 DO 320 LL = II, IIA+NP-1
643 SUM = SUM + ABS( A( LL+IOFFA ) )
644 320 CONTINUE
645 END IF
646 IOFFA = IOFFA + LDA
647 WORK( JJ+K-JJA+ICSR0 ) = SUM
648.EQ. IF( MYROWICURROW )
649 $ II = II + 1
650 330 CONTINUE
651*
652* Reset local indices so we can find ROWSUMS
653*
654.EQ. IF( MYROWICURROW )
655 $ II = II - IB
656*
657 END IF
658*
659* Find ROWSUMS
660*
661.EQ. IF( MYROWICURROW ) THEN
662 DO 350 K = II, II+IB-1
663 SUM = ZERO
664.GT. IF( JJJJA ) THEN
665 DO 340 LL = (JJA-1)*LDA, (JJ-2)*LDA, LDA
666 SUM = SUM + ABS( A( K+LL ) )
667 340 CONTINUE
668 END IF
669 WORK(K-IIA+IRSC0) = SUM
670.EQ. IF( MYCOLICURCOL )
671 $ JJ = JJ + 1
672 350 CONTINUE
673 II = II + IB
674.EQ. ELSE IF( MYCOLICURCOL ) THEN
675 JJ = JJ + IB
676 END IF
677*
678 ICURROW = MOD( ICURROW+1, NPROW )
679 ICURCOL = MOD( ICURCOL+1, NPCOL )
680*
681 360 CONTINUE
682 END IF
683*
684* After calls to DGSUM2D, process row 0 will have global
685* COLSUMS and process column 0 will have global ROWSUMS.
686* Transpose ROWSUMS and add to COLSUMS to get global row/column
687* sum, the max of which is the infinity or 1 norm.
688*
689.EQ. IF( MYCOLIACOL )
690 $ NQ = NQ + ICOFF
691 CALL DGSUM2D( ICTXT, 'columnwise', ' ', 1, NQ, WORK( ICSR ), 1,
692 $ IAROW, MYCOL )
693.EQ. IF( MYROWIAROW )
694 $ NP = NP + IROFF
695 CALL DGSUM2D( ICTXT, 'rowwise', ' ', NP, 1, WORK( IRSC ),
696 $ MAX( 1, NP ), MYROW, IACOL )
697*
698 CALL PDCOL2ROW( ICTXT, N, 1, DESCA( MB_ ), WORK( IRSC ),
699 $ MAX( 1, NP ), WORK( IRSR ), MAX( 1, NQ ),
700 $ IAROW, IACOL, IAROW, IACOL, WORK( IRSC+NP ) )
701*
702.EQ. IF( MYROWIAROW ) THEN
703.EQ. IF( MYCOLIACOL )
704 $ NQ = NQ - ICOFF
705 CALL DAXPY( NQ, ONE, WORK( IRSR0 ), 1, WORK( ICSR0 ), 1 )
706.LT. IF( NQ1 ) THEN
707 VALUE = ZERO
708 ELSE
709 VALUE = WORK( IDAMAX( NQ, WORK( ICSR0 ), 1 ) )
710 END IF
711 CALL DGAMX2D( ICTXT, 'rowwise', ' ', 1, 1, VALUE, 1, I, K,
712 $ -1, IAROW, IACOL )
713 END IF
714*
715************************************************************************
716* Frobenius norm
717* SSQ(1) is scale
718* SSQ(2) is sum-of-squares
719*
720 ELSE IF( LSAME( NORM, 'f.OR.' ) LSAME( NORM, 'e' ) ) THEN
721*
722* Find normF( sub( A ) ).
723*
724 SSQ(1) = ZERO
725 SSQ(2) = ONE
726*
727* Add off-diagonal entries, first
728*
729 IF( LSAME( UPLO, 'u' ) ) THEN
730*
731* Handle first block separately
732*
733 IB = IN-IA+1
734*
735.EQ. IF( MYCOLIACOL ) THEN
736 DO 370 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
737 COLSSQ(1) = ZERO
738 COLSSQ(2) = ONE
739 CALL DLASSQ( II-IIA, A( IIA+K ), 1,
740 $ COLSSQ(1), COLSSQ(2) )
741.EQ. IF( MYROWIAROW )
742 $ II = II + 1
743 CALL DLASSQ( II-IIA, A( IIA+K ), 1,
744 $ COLSSQ(1), COLSSQ(2) )
745 CALL DCOMBSSQ( SSQ, COLSSQ )
746 370 CONTINUE
747*
748 JJ = JJ + IB
749.EQ. ELSE IF( MYROWIAROW ) THEN
750 II = II + IB
751 END IF
752*
753 ICURROW = MOD( IAROW+1, NPROW )
754 ICURCOL = MOD( IACOL+1, NPCOL )
755*
756* Loop over rows/columns of global matrix.
757*
758 DO 390 I = IN+1, IA+N-1, DESCA( MB_ )
759 IB = MIN( DESCA( MB_ ), IA+N-I )
760*
761.EQ. IF( MYCOLICURCOL ) THEN
762 DO 380 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
763 COLSSQ(1) = ZERO
764 COLSSQ(2) = ONE
765 CALL DLASSQ( II-IIA, A( IIA+K ), 1,
766 $ COLSSQ(1), COLSSQ(2) )
767.EQ. IF( MYROWICURROW )
768 $ II = II + 1
769 CALL DLASSQ( II-IIA, A (IIA+K ), 1,
770 $ COLSSQ(1), COLSSQ(2) )
771 CALL DCOMBSSQ( SSQ, COLSSQ )
772 380 CONTINUE
773*
774 JJ = JJ + IB
775.EQ. ELSE IF( MYROWICURROW ) THEN
776 II = II + IB
777 END IF
778*
779 ICURROW = MOD( ICURROW+1, NPROW )
780 ICURCOL = MOD( ICURCOL+1, NPCOL )
781*
782 390 CONTINUE
783*
784 ELSE
785*
786* Handle first block separately
787*
788 IB = IN-IA+1
789*
790.EQ. IF( MYCOLIACOL ) THEN
791 DO 400 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
792 COLSSQ(1) = ZERO
793 COLSSQ(2) = ONE
794 CALL DLASSQ( IIA+NP-II, A( II+K ), 1,
795 $ COLSSQ(1), COLSSQ(2) )
796.EQ. IF( MYROWIAROW )
797 $ II = II + 1
798 CALL DLASSQ( IIA+NP-II, A( II+K ), 1,
799 $ COLSSQ(1), COLSSQ(2) )
800 CALL DCOMBSSQ( SSQ, COLSSQ )
801 400 CONTINUE
802*
803 JJ = JJ + IB
804.EQ. ELSE IF( MYROWIAROW ) THEN
805 II = II + IB
806 END IF
807*
808 ICURROW = MOD( IAROW+1, NPROW )
809 ICURCOL = MOD( IACOL+1, NPCOL )
810*
811* Loop over rows/columns of global matrix.
812*
813 DO 420 I = IN+1, IA+N-1, DESCA( MB_ )
814 IB = MIN( DESCA( MB_ ), IA+N-I )
815*
816.EQ. IF( MYCOLICURCOL ) THEN
817 DO 410 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
818 COLSSQ(1) = ZERO
819 COLSSQ(2) = ONE
820 CALL DLASSQ( IIA+NP-II, A( II+K ), 1,
821 $ COLSSQ(1), COLSSQ(2) )
822.EQ. IF( MYROWICURROW )
823 $ II = II + 1
824 CALL DLASSQ( IIA+NP-II, A( II+K ), 1,
825 $ COLSSQ(1), COLSSQ(2) )
826 CALL DCOMBSSQ( SSQ, COLSSQ )
827 410 CONTINUE
828*
829 JJ = JJ + IB
830.EQ. ELSE IF( MYROWICURROW ) THEN
831 II = II + IB
832 END IF
833*
834 ICURROW = MOD( ICURROW+1, NPROW )
835 ICURCOL = MOD( ICURCOL+1, NPCOL )
836*
837 420 CONTINUE
838*
839 END IF
840*
841* Perform the global scaled sum
842*
843 CALL PDTREECOMB( ICTXT, 'all', 2, SSQ, IAROW, IACOL,
844 $ DCOMBSSQ )
845 VALUE = SSQ( 1 ) * SQRT( SSQ( 2 ) )
846*
847 END IF
848*
849* Broadcast the result to the other processes
850*
851.EQ..AND..EQ. IF( MYROWIAROW MYCOLIACOL ) THEN
852 CALL DGEBS2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1 )
853 ELSE
854 CALL DGEBR2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1, IAROW,
855 $ IACOL )
856 END IF
857*
858 PDLANSY = VALUE
859*
860 RETURN
861*
862* End of PDLANSY
863*
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
dlassq
subroutine dlassq(n, x, incx, scl, sumsq)
DLASSQ updates a sum of squares represented in scaled form.
Definition dlassq.f90:137
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
idamax
integer function idamax(n, dx, incx)
IDAMAX
Definition idamax.f:71
daxpy
subroutine daxpy(n, da, dx, incx, dy, incy)
DAXPY
Definition daxpy.f:89
iceil
integer function iceil(inum, idenom)
Definition iceil.f:2
min
#define min(a, b)
Definition macros.h:20
max
#define max(a, b)
Definition macros.h:21
dgebs2d
subroutine dgebs2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1082
dgebr2d
subroutine dgebr2d(contxt, scope, top, m, n, a, lda)
Definition mpi.f:1123
infog2l
subroutine infog2l(grindx, gcindx, desc, nprow, npcol, myrow, mycol, lrindx, lcindx, rsrc, csrc)
Definition mpi.f:937
blacs_gridinfo
subroutine blacs_gridinfo(cntxt, nprow, npcol, myrow, mycol)
Definition mpi.f:754
numroc
integer function numroc(n, nb, iproc, isrcproc, nprocs)
Definition mpi.f:786
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