OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pzlanhe.f
Go to the documentation of this file.
1 DOUBLE PRECISION FUNCTION pzlanhe( NORM, UPLO, N, A, IA, JA,
2 $ DESCA, WORK )
3*
4* -- ScaLAPACK auxiliary routine (version 1.7) --
5* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
6* and University of California, Berkeley.
7* May 1, 1997
8*
9* .. Scalar Arguments ..
10 CHARACTER norm, uplo
11 INTEGER ia, ja, n
12* ..
13* .. Array Arguments ..
14 INTEGER desca( * )
15 DOUBLE PRECISION WORK( * )
16 COMPLEX*16 a( * )
17* ..
18*
19* Purpose
20* =======
21*
22* PZLANHE 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* complex hermitian distributed matrix sub(A) = A(IA:IA+N-1,JA:JA+N-1).
25*
26* PZLANHE 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 PZLANHE as described
101* above.
102*
103* UPLO (global input) CHARACTER
104* Specifies whether the upper or lower triangular part of the
105* hermitian 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, PZLANHE is set to zero. N >= 0.
113*
114* A (local input) COMPLEX*16 pointer into the local memory
115* to an array of dimension (LLD_A, LOCc(JA+N-1)) containing the
116* local pieces of the hermitian 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 absa, scale, sum, value
178* ..
179* .. Local Arrays ..
180 DOUBLE PRECISION rwork( 2 )
181* ..
182* .. External Subroutines ..
183 EXTERNAL blacs_gridinfo, daxpy, dcombssq,
184 $ dgamx2d, dgsum2d, dgebr2d,
185 $ dgebs2d, pdcol2row, pdtreecomb,
186 $ zlassq
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, dble, 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 Hermitian, 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 ELSE IF( lsame( norm, 'M' ) ) THEN
269*
270* Find max(abs(A(i,j))).
271*
272 VALUE = zero
273*
274 IF( lsame( uplo, 'u' ) ) THEN
275*
276* Handle first block separately
277*
278 IB = IN-IA+1
279*
280* Find COLMAXS
281*
282.EQ. IF( MYCOLIACOL ) THEN
283 DO 20 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
284.GT. IF( IIIIA ) THEN
285 DO 10 LL = IIA, II-1
286 VALUE = MAX( VALUE, ABS( A( LL+K ) ) )
287 10 CONTINUE
288 END IF
289.EQ. IF( MYROWIAROW )
290 $ II = II + 1
291 20 CONTINUE
292*
293* Reset local indices so we can find ROWMAXS
294*
295.EQ. IF( MYROWIAROW )
296 $ II = II - IB
297*
298 END IF
299*
300* Find ROWMAXS
301*
302.EQ. IF( MYROWIAROW ) THEN
303 DO 40 K = II, II+IB-1
304.EQ. IF( MYCOLIACOL ) THEN
305.LE. IF( JJJJA+NQ-1 ) THEN
306 VALUE = MAX( VALUE,
307 $ ABS( DBLE( A( K+(JJ-1)*LDA ) ) ) )
308 DO 30 LL = JJ*LDA, (JJA+NQ-2)*LDA, LDA
309 VALUE = MAX( VALUE, ABS( A( K+LL ) ) )
310 30 CONTINUE
311 END IF
312 ELSE
313.LE. IF( JJJJA+NQ-1 ) THEN
314 DO 35 LL = (JJ-1)*LDA, (JJA+NQ-2)*LDA, LDA
315 VALUE = MAX( VALUE, ABS( A( K+LL ) ) )
316 35 CONTINUE
317 END IF
318 END IF
319.EQ. IF( MYCOLIACOL )
320 $ JJ = JJ + 1
321 40 CONTINUE
322 II = II + IB
323.EQ. ELSE IF( MYCOLIACOL ) THEN
324 JJ = JJ + IB
325 END IF
326*
327 ICURROW = MOD( IAROW+1, NPROW )
328 ICURCOL = MOD( IACOL+1, NPCOL )
329*
330* Loop over the remaining rows/columns of the matrix.
331*
332 DO 90 I = IN+1, IA+N-1, DESCA( MB_ )
333 IB = MIN( DESCA( MB_ ), IA+N-I )
334*
335* Find COLMAXS
336*
337.EQ. IF( MYCOLICURCOL ) THEN
338 DO 60 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
339.GT. IF( IIIIA ) THEN
340 DO 50 LL = IIA, II-1
341 VALUE = MAX( VALUE, ABS( A( LL+K ) ) )
342 50 CONTINUE
343 END IF
344.EQ. IF( MYROWICURROW )
345 $ II = II + 1
346 60 CONTINUE
347*
348* Reset local indices so we can find ROWMAXS
349*
350.EQ. IF( MYROWICURROW )
351 $ II = II - IB
352 END IF
353*
354* Find ROWMAXS
355*
356.EQ. IF( MYROWICURROW ) THEN
357 DO 80 K = II, II+IB-1
358.EQ. IF( MYCOLICURCOL ) THEN
359.LE. IF( JJJJA+NQ-1 ) THEN
360 VALUE = MAX( VALUE,
361 $ ABS( DBLE( A( K+(JJ-1)*LDA ) ) ) )
362 DO 70 LL = JJ*LDA, (JJA+NQ-2)*LDA, LDA
363 VALUE = MAX( VALUE, ABS( A( K+LL ) ) )
364 70 CONTINUE
365 END IF
366 ELSE
367.LE. IF( JJJJA+NQ-1 ) THEN
368 DO 75 LL = (JJ-1)*LDA, (JJA+NQ-2)*LDA, LDA
369 VALUE = MAX( VALUE, ABS( A( K+LL ) ) )
370 75 CONTINUE
371 END IF
372 END IF
373.EQ. IF( MYCOLICURCOL )
374 $ JJ = JJ + 1
375 80 CONTINUE
376 II = II + IB
377.EQ. ELSE IF( MYCOLICURCOL ) THEN
378 JJ = JJ + IB
379 END IF
380 ICURROW = MOD( ICURROW+1, NPROW )
381 ICURCOL = MOD( ICURCOL+1, NPCOL )
382 90 CONTINUE
383*
384 ELSE
385*
386* Handle first block separately
387*
388 IB = IN-IA+1
389*
390* Find COLMAXS
391*
392.EQ. IF( MYCOLIACOL ) THEN
393 DO 110 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
394.EQ. IF( MYROWIAROW ) THEN
395.LE. IF( IIIIA+NP-1 ) THEN
396 VALUE = MAX( VALUE, ABS( DBLE( A( II+K ) ) ) )
397 DO 100 LL = II+1, IIA+NP-1
398 VALUE = MAX( VALUE, ABS( A( LL+K ) ) )
399 100 CONTINUE
400 END IF
401 ELSE
402.LE. IF( IIIIA+NP-1 ) THEN
403 DO 105 LL = II, IIA+NP-1
404 VALUE = MAX( VALUE, ABS( A( LL+K ) ) )
405 105 CONTINUE
406 END IF
407 END IF
408.EQ. IF( MYROWIAROW )
409 $ II = II + 1
410 110 CONTINUE
411*
412* Reset local indices so we can find ROWMAXS
413*
414.EQ. IF( MYROWIAROW )
415 $ II = II - IB
416 END IF
417*
418* Find ROWMAXS
419*
420.EQ. IF( MYROWIAROW ) THEN
421 DO 130 K = 0, IB-1
422.GT. IF( JJJJA ) THEN
423 DO 120 LL = (JJA-1)*LDA, (JJ-2)*LDA, LDA
424 VALUE = MAX( VALUE, ABS( A( II+LL ) ) )
425 120 CONTINUE
426 END IF
427 II = II + 1
428.EQ. IF( MYCOLIACOL )
429 $ JJ = JJ + 1
430 130 CONTINUE
431.EQ. ELSE IF( MYCOLIACOL ) THEN
432 JJ = JJ + IB
433 END IF
434*
435 ICURROW = MOD( IAROW+1, NPROW )
436 ICURCOL = MOD( IACOL+1, NPCOL )
437*
438* Loop over rows/columns of global matrix.
439*
440 DO 180 I = IN+1, IA+N-1, DESCA( MB_ )
441 IB = MIN( DESCA( MB_ ), IA+N-I )
442*
443* Find COLMAXS
444*
445.EQ. IF( MYCOLICURCOL ) THEN
446 DO 150 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
447.EQ. IF( MYROWICURROW ) THEN
448.LE. IF( IIIIA+NP-1 ) THEN
449 VALUE = MAX( VALUE,
450 $ ABS( DBLE( A( II+K ) ) ) )
451 DO 140 LL = II+1, IIA+NP-1
452 VALUE = MAX( VALUE, ABS( A( LL+K ) ) )
453 140 CONTINUE
454 END IF
455 ELSE
456.LE. IF( IIIIA+NP-1 ) THEN
457 DO 145 LL = II, IIA+NP-1
458 VALUE = MAX( VALUE, ABS( A( LL+K ) ) )
459 145 CONTINUE
460 END IF
461 END IF
462.EQ. IF( MYROWICURROW )
463 $ II = II + 1
464 150 CONTINUE
465*
466* Reset local indices so we can find ROWMAXS
467*
468.EQ. IF( MYROWICURROW )
469 $ II = II - IB
470 END IF
471*
472* Find ROWMAXS
473*
474.EQ. IF( MYROWICURROW ) THEN
475 DO 170 K = 0, IB-1
476.GT. IF( JJJJA ) THEN
477 DO 160 LL = (JJA-1)*LDA, (JJ-2)*LDA, LDA
478 VALUE = MAX( VALUE, ABS( A( II+LL ) ) )
479 160 CONTINUE
480 END IF
481 II = II + 1
482.EQ. IF( MYCOLICURCOL )
483 $ JJ = JJ + 1
484 170 CONTINUE
485.EQ. ELSE IF( MYCOLICURCOL ) THEN
486 JJ = JJ + IB
487 END IF
488 ICURROW = MOD( ICURROW+1, NPROW )
489 ICURCOL = MOD( ICURCOL+1, NPCOL )
490*
491 180 CONTINUE
492*
493 END IF
494*
495* Gather the result on process (IAROW,IACOL).
496*
497 CALL DGAMX2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1, I, K, -1,
498 $ IAROW, IACOL )
499*
500 ELSE IF( LSAME( NORM, 'i.OR.' ) LSAME( NORM, 'o.OR.' )
501.EQ. $ NORM'1' ) THEN
502*
503* Find normI( sub( A ) ) ( = norm1( sub( A ) ), since sub( A ) is
504* hermitian).
505*
506 IF( LSAME( UPLO, 'u' ) ) THEN
507*
508* Handle first block separately
509*
510 IB = IN-IA+1
511*
512* Find COLSUMS
513*
514.EQ. IF( MYCOLIACOL ) THEN
515 IOFFA = ( JJ - 1 ) * LDA
516 DO 200 K = 0, IB-1
517 SUM = ZERO
518.GT. IF( IIIIA ) THEN
519 DO 190 LL = IIA, II-1
520 SUM = SUM + ABS( A( LL+IOFFA ) )
521 190 CONTINUE
522 END IF
523 IOFFA = IOFFA + LDA
524 WORK( JJ+K-JJA+ICSR0 ) = SUM
525.EQ. IF( MYROWIAROW )
526 $ II = II + 1
527 200 CONTINUE
528*
529* Reset local indices so we can find ROWSUMS
530*
531.EQ. IF( MYROWIAROW )
532 $ II = II - IB
533*
534 END IF
535*
536* Find ROWSUMS
537*
538.EQ. IF( MYROWIAROW ) THEN
539 DO 220 K = II, II+IB-1
540 SUM = ZERO
541.EQ. IF( MYCOLIACOL ) THEN
542.GT. IF( JJA+NQJJ ) THEN
543 SUM = ABS( DBLE( A( K+(JJ-1)*LDA ) ) )
544 DO 210 LL = JJ*LDA, (JJA+NQ-2)*LDA, LDA
545 SUM = SUM + ABS( A( K+LL ) )
546 210 CONTINUE
547 END IF
548 ELSE
549.GT. IF( JJA+NQJJ ) THEN
550 DO 215 LL = (JJ-1)*LDA, (JJA+NQ-2)*LDA, LDA
551 SUM = SUM + ABS( A( K+LL ) )
552 215 CONTINUE
553 END IF
554 END IF
555 WORK( K-IIA+IRSC0 ) = SUM
556.EQ. IF( MYCOLIACOL )
557 $ JJ = JJ + 1
558 220 CONTINUE
559 II = II + IB
560.EQ. ELSE IF( MYCOLIACOL ) THEN
561 JJ = JJ + IB
562 END IF
563*
564 ICURROW = MOD( IAROW+1, NPROW )
565 ICURCOL = MOD( IACOL+1, NPCOL )
566*
567* Loop over remaining rows/columns of global matrix.
568*
569 DO 270 I = IN+1, IA+N-1, DESCA( MB_ )
570 IB = MIN( DESCA( MB_ ), IA+N-I )
571*
572* Find COLSUMS
573*
574.EQ. IF( MYCOLICURCOL ) THEN
575 IOFFA = ( JJ - 1 ) * LDA
576 DO 240 K = 0, IB-1
577 SUM = ZERO
578.GT. IF( IIIIA ) THEN
579 DO 230 LL = IIA, II-1
580 SUM = SUM + ABS( A( IOFFA+LL ) )
581 230 CONTINUE
582 END IF
583 IOFFA = IOFFA + LDA
584 WORK( JJ+K-JJA+ICSR0 ) = SUM
585.EQ. IF( MYROWICURROW )
586 $ II = II + 1
587 240 CONTINUE
588*
589* Reset local indices so we can find ROWSUMS
590*
591.EQ. IF( MYROWICURROW )
592 $ II = II - IB
593*
594 END IF
595*
596* Find ROWSUMS
597*
598.EQ. IF( MYROWICURROW ) THEN
599 DO 260 K = II, II+IB-1
600 SUM = ZERO
601.EQ. IF( MYCOLICURCOL ) THEN
602.GT. IF( JJA+NQJJ ) THEN
603 SUM = ABS( DBLE( A( K+(JJ-1)*LDA ) ) )
604 DO 250 LL = JJ*LDA, (JJA+NQ-2)*LDA, LDA
605 SUM = SUM + ABS( A( K+LL ) )
606 250 CONTINUE
607 END IF
608 ELSE
609.GT. IF( JJA+NQJJ ) THEN
610 DO 255 LL = (JJ-1)*LDA, (JJA+NQ-2)*LDA, LDA
611 SUM = SUM + ABS( A( K+LL ) )
612 255 CONTINUE
613 END IF
614 END IF
615 WORK( K-IIA+IRSC0 ) = SUM
616.EQ. IF( MYCOLICURCOL )
617 $ JJ = JJ + 1
618 260 CONTINUE
619 II = II + IB
620.EQ. ELSE IF( MYCOLICURCOL ) THEN
621 JJ = JJ + IB
622 END IF
623*
624 ICURROW = MOD( ICURROW+1, NPROW )
625 ICURCOL = MOD( ICURCOL+1, NPCOL )
626*
627 270 CONTINUE
628*
629 ELSE
630*
631* Handle first block separately
632*
633 IB = IN-IA+1
634*
635* Find COLSUMS
636*
637.EQ. IF( MYCOLIACOL ) THEN
638 IOFFA = (JJ-1)*LDA
639 DO 290 K = 0, IB-1
640 SUM = ZERO
641.EQ. IF( MYROWIAROW ) THEN
642.GT. IF( IIA+NPII ) THEN
643 SUM = ABS( DBLE( A( IOFFA+II ) ) )
644 DO 280 LL = II+1, IIA+NP-1
645 SUM = SUM + ABS( A( IOFFA+LL ) )
646 280 CONTINUE
647 END IF
648 ELSE
649 DO 285 LL = II, IIA+NP-1
650 SUM = SUM + ABS( A( IOFFA+LL ) )
651 285 CONTINUE
652 END IF
653 IOFFA = IOFFA + LDA
654 WORK( JJ+K-JJA+ICSR0 ) = SUM
655.EQ. IF( MYROWIAROW )
656 $ II = II + 1
657 290 CONTINUE
658*
659* Reset local indices so we can find ROWSUMS
660*
661.EQ. IF( MYROWIAROW )
662 $ II = II - IB
663*
664 END IF
665*
666* Find ROWSUMS
667*
668.EQ. IF( MYROWIAROW ) THEN
669 DO 310 K = II, II+IB-1
670 SUM = ZERO
671.GT. IF( JJJJA ) THEN
672 DO 300 LL = (JJA-1)*LDA, (JJ-2)*LDA, LDA
673 SUM = SUM + ABS( A( K+LL ) )
674 300 CONTINUE
675 END IF
676 WORK( K-IIA+IRSC0 ) = SUM
677.EQ. IF( MYCOLIACOL )
678 $ JJ = JJ + 1
679 310 CONTINUE
680 II = II + IB
681.EQ. ELSE IF( MYCOLIACOL ) THEN
682 JJ = JJ + IB
683 END IF
684*
685 ICURROW = MOD( IAROW+1, NPROW )
686 ICURCOL = MOD( IACOL+1, NPCOL )
687*
688* Loop over rows/columns of global matrix.
689*
690 DO 360 I = IN+1, IA+N-1, DESCA( MB_ )
691 IB = MIN( DESCA( MB_ ), IA+N-I )
692*
693* Find COLSUMS
694*
695.EQ. IF( MYCOLICURCOL ) THEN
696 IOFFA = ( JJ - 1 ) * LDA
697 DO 330 K = 0, IB-1
698 SUM = ZERO
699.EQ. IF( MYROWICURROW ) THEN
700.GT. IF( IIA+NPII ) THEN
701 SUM = ABS( DBLE( A( II+IOFFA ) ) )
702 DO 320 LL = II+1, IIA+NP-1
703 SUM = SUM + ABS( A( LL+IOFFA ) )
704 320 CONTINUE
705.EQ. ELSE IF( IIIIA+NP-1 ) THEN
706 SUM = ABS( DBLE( A( II+IOFFA ) ) )
707 END IF
708 ELSE
709 DO 325 LL = II, IIA+NP-1
710 SUM = SUM + ABS( A( LL+IOFFA ) )
711 325 CONTINUE
712 END IF
713 IOFFA = IOFFA + LDA
714 WORK( JJ+K-JJA+ICSR0 ) = SUM
715.EQ. IF( MYROWICURROW )
716 $ II = II + 1
717 330 CONTINUE
718*
719* Reset local indices so we can find ROWSUMS
720*
721.EQ. IF( MYROWICURROW )
722 $ II = II - IB
723*
724 END IF
725*
726* Find ROWSUMS
727*
728.EQ. IF( MYROWICURROW ) THEN
729 DO 350 K = II, II+IB-1
730 SUM = ZERO
731.GT. IF( JJJJA ) THEN
732 DO 340 LL = (JJA-1)*LDA, (JJ-2)*LDA, LDA
733 SUM = SUM + ABS( A( K+LL ) )
734 340 CONTINUE
735 END IF
736 WORK(K-IIA+IRSC0) = SUM
737.EQ. IF( MYCOLICURCOL )
738 $ JJ = JJ + 1
739 350 CONTINUE
740 II = II + IB
741.EQ. ELSE IF( MYCOLICURCOL ) THEN
742 JJ = JJ + IB
743 END IF
744*
745 ICURROW = MOD( ICURROW+1, NPROW )
746 ICURCOL = MOD( ICURCOL+1, NPCOL )
747*
748 360 CONTINUE
749 END IF
750*
751* After calls to DGSUM2D, process row 0 will have global
752* COLSUMS and process column 0 will have global ROWSUMS.
753* Transpose ROWSUMS and add to COLSUMS to get global row/column
754* sum, the max of which is the infinity or 1 norm.
755*
756.EQ. IF( MYCOLIACOL )
757 $ NQ = NQ + ICOFF
758 CALL DGSUM2D( ICTXT, 'columnwise', ' ', 1, NQ, WORK( ICSR ), 1,
759 $ IAROW, MYCOL )
760.EQ. IF( MYROWIAROW )
761 $ NP = NP + IROFF
762 CALL DGSUM2D( ICTXT, 'rowwise', ' ', NP, 1, WORK( IRSC ),
763 $ MAX( 1, NP ), MYROW, IACOL )
764*
765 CALL PDCOL2ROW( ICTXT, N, 1, DESCA( MB_ ), WORK( IRSC ),
766 $ MAX( 1, NP ), WORK( IRSR ), MAX( 1, NQ ),
767 $ IAROW, IACOL, IAROW, IACOL, WORK( IRSC+NP ) )
768*
769.EQ. IF( MYROWIAROW ) THEN
770.EQ. IF( MYCOLIACOL )
771 $ NQ = NQ - ICOFF
772 CALL DAXPY( NQ, ONE, WORK( IRSR0 ), 1, WORK( ICSR0 ), 1 )
773.LT. IF( NQ1 ) THEN
774 VALUE = ZERO
775 ELSE
776 VALUE = WORK( IDAMAX( NQ, WORK( ICSR0 ), 1 ) )
777 END IF
778 CALL DGAMX2D( ICTXT, 'rowwise', ' ', 1, 1, VALUE, 1, I, K,
779 $ -1, IAROW, IACOL )
780 END IF
781*
782 ELSE IF( LSAME( NORM, 'f.OR.' ) LSAME( NORM, 'e' ) ) THEN
783*
784* Find normF( sub( A ) ).
785*
786 SCALE = ZERO
787 SUM = ONE
788*
789* Add off-diagonal entries, first
790*
791 IF( LSAME( UPLO, 'u' ) ) THEN
792*
793* Handle first block separately
794*
795 IB = IN-IA+1
796*
797.EQ. IF( MYCOLIACOL ) THEN
798 DO 370 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
799 CALL ZLASSQ( II-IIA, A( IIA+K ), 1, SCALE, SUM )
800 CALL ZLASSQ( II-IIA, A( IIA+K ), 1, SCALE, SUM )
801.EQ. IF( MYROWIAROW ) THEN
802.NE. IF( DBLE( A( II+K ) )ZERO ) THEN
803 ABSA = ABS( DBLE( A( II+K ) ) )
804.LT. IF( SCALEABSA ) THEN
805 SUM = ONE + SUM * ( SCALE / ABSA )**2
806 SCALE = ABSA
807 ELSE
808 SUM = SUM + ( ABSA / SCALE )**2
809 END IF
810 END IF
811 II = II + 1
812 END IF
813 370 CONTINUE
814*
815 JJ = JJ + IB
816.EQ. ELSE IF( MYROWIAROW ) THEN
817 II = II + IB
818 END IF
819*
820 ICURROW = MOD( IAROW+1, NPROW )
821 ICURCOL = MOD( IACOL+1, NPCOL )
822*
823* Loop over rows/columns of global matrix.
824*
825 DO 390 I = IN+1, IA+N-1, DESCA( MB_ )
826 IB = MIN( DESCA( MB_ ), IA+N-I )
827*
828.EQ. IF( MYCOLICURCOL ) THEN
829 DO 380 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
830 CALL ZLASSQ( II-IIA, A( IIA+K ), 1, SCALE, SUM )
831 CALL ZLASSQ( II-IIA, A( IIA+K ), 1, SCALE, SUM )
832.EQ. IF( MYROWICURROW ) THEN
833.NE. IF( DBLE( A( II+K ) )ZERO ) THEN
834 ABSA = ABS( DBLE( A( II+K ) ) )
835.LT. IF( SCALEABSA ) THEN
836 SUM = ONE + SUM * ( SCALE / ABSA )**2
837 SCALE = ABSA
838 ELSE
839 SUM = SUM + ( ABSA / SCALE )**2
840 END IF
841 END IF
842 II = II + 1
843 END IF
844 380 CONTINUE
845*
846 JJ = JJ + IB
847.EQ. ELSE IF( MYROWICURROW ) THEN
848 II = II + IB
849 END IF
850*
851 ICURROW = MOD( ICURROW+1, NPROW )
852 ICURCOL = MOD( ICURCOL+1, NPCOL )
853*
854 390 CONTINUE
855*
856 ELSE
857*
858* Handle first block separately
859*
860 IB = IN-IA+1
861*
862.EQ. IF( MYCOLIACOL ) THEN
863 DO 400 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
864.EQ. IF( MYROWIAROW ) THEN
865.NE. IF( DBLE( A( II+K ) )ZERO ) THEN
866 ABSA = ABS( DBLE( A( II+K ) ) )
867.LT. IF( SCALEABSA ) THEN
868 SUM = ONE + SUM * ( SCALE / ABSA )**2
869 SCALE = ABSA
870 ELSE
871 SUM = SUM + ( ABSA / SCALE )**2
872 END IF
873 END IF
874 II = II + 1
875 END IF
876 CALL ZLASSQ( IIA+NP-II, A( II+K ), 1, SCALE, SUM )
877 CALL ZLASSQ( IIA+NP-II, A( II+K ), 1, SCALE, SUM )
878 400 CONTINUE
879*
880 JJ = JJ + IB
881.EQ. ELSE IF( MYROWIAROW ) THEN
882 II = II + IB
883 END IF
884*
885 ICURROW = MOD( IAROW+1, NPROW )
886 ICURCOL = MOD( IACOL+1, NPCOL )
887*
888* Loop over rows/columns of global matrix.
889*
890 DO 420 I = IN+1, IA+N-1, DESCA( MB_ )
891 IB = MIN( DESCA( MB_ ), IA+N-I )
892*
893.EQ. IF( MYCOLICURCOL ) THEN
894 DO 410 K = (JJ-1)*LDA, (JJ+IB-2)*LDA, LDA
895.EQ. IF( MYROWICURROW ) THEN
896.NE. IF( DBLE( A( II+K ) )ZERO ) THEN
897 ABSA = ABS( DBLE( A( II+K ) ) )
898.LT. IF( SCALEABSA ) THEN
899 SUM = ONE + SUM * ( SCALE / ABSA )**2
900 SCALE = ABSA
901 ELSE
902 SUM = SUM + ( ABSA / SCALE )**2
903 END IF
904 END IF
905 II = II + 1
906 END IF
907 CALL ZLASSQ( IIA+NP-II, A( II+K ), 1, SCALE, SUM )
908 CALL ZLASSQ( IIA+NP-II, A( II+K ), 1, SCALE, SUM )
909 410 CONTINUE
910*
911 JJ = JJ + IB
912.EQ. ELSE IF( MYROWICURROW ) THEN
913 II = II + IB
914 END IF
915*
916 ICURROW = MOD( ICURROW+1, NPROW )
917 ICURCOL = MOD( ICURCOL+1, NPCOL )
918*
919 420 CONTINUE
920*
921 END IF
922*
923* Perform the global scaled sum
924*
925 RWORK( 1 ) = SCALE
926 RWORK( 2 ) = SUM
927*
928 CALL PDTREECOMB( ICTXT, 'all', 2, RWORK, IAROW, IACOL,
929 $ DCOMBSSQ )
930 VALUE = RWORK( 1 ) * SQRT( RWORK( 2 ) )
931*
932 END IF
933*
934* Broadcast the result to the other processes
935*
936.EQ..AND..EQ. IF( MYROWIAROW MYCOLIACOL ) THEN
937 CALL DGEBS2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1 )
938 ELSE
939 CALL DGEBR2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1, IAROW,
940 $ IACOL )
941 END IF
942*
943 PZLANHE = VALUE
944*
945 RETURN
946*
947* End of PZLANHE
948*
949 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
zlassq
subroutine zlassq(n, x, incx, scl, sumsq)
ZLASSQ updates a sum of squares represented in scaled form.
Definition zlassq.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
pzlanhe
double precision function pzlanhe(norm, uplo, n, a, ia, ja, desca, work)
Definition pzlanhe.f:3