OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
pzlantr.f
Go to the documentation of this file.
1 DOUBLE PRECISION FUNCTION pzlantr( NORM, UPLO, DIAG, M, N, A,
2 $ IA, JA, DESCA, WORK )
3 IMPLICIT NONE
4*
5* -- ScaLAPACK auxiliary routine (version 1.7) --
6* University of Tennessee, Knoxville, Oak Ridge National Laboratory,
7* and University of California, Berkeley.
8* May 1, 1997
9*
10* .. Scalar Arguments ..
11 CHARACTER diag, norm, uplo
12 INTEGER ia, ja, m, n
13* ..
14* .. Array Arguments ..
15 INTEGER desca( * )
16 DOUBLE PRECISION work( * )
17 COMPLEX*16 a( * )
18* ..
19*
20* Purpose
21* =======
22*
23* PZLANTR 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* trapezoidal or triangular distributed matrix sub( A ) denoting
26* A(IA:IA+M-1, JA:JA+N-1).
27*
28* PZLANTR returns the value
29*
30* ( max(abs(A(i,j))), NORM = 'M' or 'm' with ia <= i <= ia+m-1,
31* ( and ja <= j <= ja+n-1,
32* (
33* ( norm1( sub( A ) ), NORM = '1', 'O' or 'o'
34* (
35* ( normI( sub( A ) ), NORM = 'I' or 'i'
36* (
37* ( normF( sub( A ) ), NORM = 'F', 'f', 'E' or 'e'
38*
39* where norm1 denotes the one norm of a matrix (maximum column sum),
40* normI denotes the infinity norm of a matrix (maximum row sum) and
41* normF denotes the Frobenius norm of a matrix (square root of sum of
42* squares). Note that max(abs(A(i,j))) is not a matrix norm.
43*
44* Notes
45* =====
46*
47* Each global data object is described by an associated description
48* vector. This vector stores the information required to establish
49* the mapping between an object element and its corresponding process
50* and memory location.
51*
52* Let A be a generic term for any 2D block cyclicly distributed array.
53* Such a global array has an associated description vector DESCA.
54* In the following comments, the character _ should be read as
55* "of the global array".
56*
57* NOTATION STORED IN EXPLANATION
58* --------------- -------------- --------------------------------------
59* DTYPE_A(global) DESCA( DTYPE_ )The descriptor type. In this case,
60* DTYPE_A = 1.
61* CTXT_A (global) DESCA( CTXT_ ) The BLACS context handle, indicating
62* the BLACS process grid A is distribu-
63* ted over. The context itself is glo-
64* bal, but the handle (the integer
65* value) may vary.
66* M_A (global) DESCA( M_ ) The number of rows in the global
67* array A.
68* N_A (global) DESCA( N_ ) The number of columns in the global
69* array A.
70* MB_A (global) DESCA( MB_ ) The blocking factor used to distribute
71* the rows of the array.
72* NB_A (global) DESCA( NB_ ) The blocking factor used to distribute
73* the columns of the array.
74* RSRC_A (global) DESCA( RSRC_ ) The process row over which the first
75* row of the array A is distributed.
76* CSRC_A (global) DESCA( CSRC_ ) The process column over which the
77* first column of the array A is
78* distributed.
79* LLD_A (local) DESCA( LLD_ ) The leading dimension of the local
80* array. LLD_A >= MAX(1,LOCr(M_A)).
81*
82* Let K be the number of rows or columns of a distributed matrix,
83* and assume that its process grid has dimension p x q.
84* LOCr( K ) denotes the number of elements of K that a process
85* would receive if K were distributed over the p processes of its
86* process column.
87* Similarly, LOCc( K ) denotes the number of elements of K that a
88* process would receive if K were distributed over the q processes of
89* its process row.
90* The values of LOCr() and LOCc() may be determined via a call to the
91* ScaLAPACK tool function, NUMROC:
92* LOCr( M ) = NUMROC( M, MB_A, MYROW, RSRC_A, NPROW ),
93* LOCc( N ) = NUMROC( N, NB_A, MYCOL, CSRC_A, NPCOL ).
94* An upper bound for these quantities may be computed by:
95* LOCr( M ) <= ceil( ceil(M/MB_A)/NPROW )*MB_A
96* LOCc( N ) <= ceil( ceil(N/NB_A)/NPCOL )*NB_A
97*
98* Arguments
99* =========
100*
101* NORM (global input) CHARACTER
102* Specifies the value to be returned in PZLANTR as described
103* above.
104*
105* UPLO (global input) CHARACTER
106* Specifies whether the matrix sub( A ) is upper or lower
107* trapezoidal.
108* = 'U': Upper trapezoidal
109* = 'L': Lower trapezoidal
110* Note that sub( A ) is triangular instead of trapezoidal
111* if M = N.
112*
113* DIAG (global input) CHARACTER
114* Specifies whether or not the distributed matrix sub( A ) has
115* unit diagonal.
116* = 'N': Non-unit diagonal
117* = 'U': Unit diagonal
118*
119* M (global input) INTEGER
120* The number of rows to be operated on i.e the number of rows
121* of the distributed submatrix sub( A ). When M = 0, PZLANTR is
122* set to zero. M >= 0.
123*
124* N (global input) INTEGER
125* The number of columns to be operated on i.e the number of
126* columns of the distributed submatrix sub( A ). When N = 0,
127* PZLANTR is set to zero. N >= 0.
128*
129* A (local input) COMPLEX*16 pointer into the local memory
130* to an array of dimension (LLD_A, LOCc(JA+N-1) ) containing
131* the local pieces of sub( A ).
132*
133* IA (global input) INTEGER
134* The row index in the global array A indicating the first
135* row of sub( A ).
136*
137* JA (global input) INTEGER
138* The column index in the global array A indicating the
139* first column of sub( A ).
140*
141* DESCA (global and local input) INTEGER array of dimension DLEN_.
142* The array descriptor for the distributed matrix A.
143*
144* WORK (local workspace) DOUBLE PRECISION array dimension (LWORK)
145* LWORK >= 0 if NORM = 'M' or 'm' (not referenced),
146* Nq0 if NORM = '1', 'O' or 'o',
147* Mp0 if NORM = 'I' or 'i',
148* 0 if NORM = 'F', 'f', 'E' or 'e' (not referenced),
149* where
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* Mp0 = NUMROC( M+IROFFA, MB_A, MYROW, IAROW, NPROW ),
155* Nq0 = NUMROC( N+ICOFFA, NB_A, MYCOL, IACOL, NPCOL ),
156*
157* INDXG2P and NUMROC are ScaLAPACK tool functions; MYROW,
158* MYCOL, NPROW and NPCOL can be determined by calling the
159* 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 LOGICAL udiag
174 INTEGER iacol, iarow, ictxt, ii, iia, icoff, ioffa,
175 $ iroff, j, jb, jj, jja, jn, kk, lda, ll, mp,
176 $ 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, dcombssq, dgebr2d,
184 $ dgebs2d, dgamx2d, dgsum2d, infog2l,
185 $ pdtreecomb, zlassq
186* ..
187* .. External Functions ..
188 LOGICAL lsame
189 INTEGER iceil, idamax, numroc
190 EXTERNAL lsame, iceil, idamax, numroc
191* ..
192* .. Intrinsic Functions ..
193 INTRINSIC abs, dble, max, min, mod, sqrt
194* ..
195* .. Executable Statements ..
196*
197* Get grid parameters
198*
199 ictxt = desca( ctxt_ )
200 CALL blacs_gridinfo( ictxt, nprow, npcol, myrow, mycol )
201*
202 udiag = lsame( diag, 'U' )
203 CALL infog2l( ia, ja, desca, nprow, npcol, myrow, mycol, iia, jja,
204 $ iarow, iacol )
205 iroff = mod( ia-1, desca( mb_ ) )
206 icoff = mod( ja-1, desca( nb_ ) )
207 mp = numroc( m+iroff, desca( mb_ ), myrow, iarow, nprow )
208 nq = numroc( n+icoff, desca( nb_ ), mycol, iacol, npcol )
209 IF( myrow.EQ.iarow )
210 $ mp = mp - iroff
211 IF( mycol.EQ.iacol )
212 $ nq = nq - icoff
213 lda = desca( lld_ )
214 ioffa = ( jja - 1 ) * lda
215*
216 IF( min( m, n ).EQ.0 ) THEN
217*
218 VALUE = zero
219*
220************************************************************************
221* max norm
222*
223 ELSE IF( lsame( norm, 'M' ) ) THEN
224*
225* Find max(abs(A(i,j))).
226*
227 IF( udiag ) THEN
228 VALUE = one
229 ELSE
230 VALUE = zero
231 END IF
232*
233 IF( lsame( uplo, 'U' ) ) THEN
234*
235* Upper triangular matrix
236*
237 ii = iia
238 jj = jja
239 jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
240 jb = jn-ja+1
241*
242 IF( mycol.EQ.iacol ) THEN
243 IF( myrow.EQ.iarow ) THEN
244 IF( udiag ) THEN
245 DO 20 ll = jj, jj + jb -1
246 DO 10 kk = iia, min(ii+ll-jj-1,iia+mp-1)
247 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
248 10 CONTINUE
249 ioffa = ioffa + lda
250 20 CONTINUE
251 ELSE
252 DO 40 ll = jj, jj + jb -1
253 DO 30 kk = iia, min( ii+ll-jj, iia+mp-1 )
254 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
255 30 CONTINUE
256 ioffa = ioffa + lda
257 40 CONTINUE
258 END IF
259 ELSE
260 DO 60 ll = jj, jj + jb -1
261 DO 50 kk = iia, min( ii-1, iia+mp-1 )
262 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
263 50 CONTINUE
264 ioffa = ioffa + lda
265 60 CONTINUE
266 END IF
267 jj = jj + jb
268 END IF
269*
270 IF( myrow.EQ.iarow )
271 $ ii = ii + jb
272 iarow = mod( iarow+1, nprow )
273 iacol = mod( iacol+1, npcol )
274*
275* Loop over remaining block of columns
276*
277 DO 130 j = jn+1, ja+n-1, desca( nb_ )
278 jb = min( ja+n-j, desca( nb_ ) )
279*
280 IF( mycol.EQ.iacol ) THEN
281 IF( myrow.EQ.iarow ) THEN
282 IF( udiag ) THEN
283 DO 80 ll = jj, jj + jb -1
284 DO 70 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
285 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
286 70 CONTINUE
287 ioffa = ioffa + lda
288 80 CONTINUE
289 ELSE
290 DO 100 ll = jj, jj + jb -1
291 DO 90 kk = iia, min( ii+ll-jj, iia+mp-1 )
292 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
293 90 CONTINUE
294 ioffa = ioffa + lda
295 100 CONTINUE
296 END IF
297 ELSE
298 DO 120 ll = jj, jj + jb -1
299 DO 110 kk = iia, min( ii-1, iia+mp-1 )
300 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
301 110 CONTINUE
302 ioffa = ioffa + lda
303 120 CONTINUE
304 END IF
305 jj = jj + jb
306 END IF
307*
308 IF( myrow.EQ.iarow )
309 $ ii = ii + jb
310 iarow = mod( iarow+1, nprow )
311 iacol = mod( iacol+1, npcol )
312*
313 130 CONTINUE
314*
315 ELSE
316*
317* Lower triangular matrix
318*
319 ii = iia
320 jj = jja
321 jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
322 jb = jn-ja+1
323*
324 IF( mycol.EQ.iacol ) THEN
325 IF( myrow.EQ.iarow ) THEN
326 IF( udiag ) THEN
327 DO 150 ll = jj, jj + jb -1
328 DO 140 kk = ii+ll-jj+1, iia+mp-1
329 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
330 140 CONTINUE
331 ioffa = ioffa + lda
332 150 CONTINUE
333 ELSE
334 DO 170 ll = jj, jj + jb -1
335 DO 160 kk = ii+ll-jj, iia+mp-1
336 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
337 160 CONTINUE
338 ioffa = ioffa + lda
339 170 CONTINUE
340 END IF
341 ELSE
342 DO 190 ll = jj, jj + jb -1
343 DO 180 kk = ii, iia+mp-1
344 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
345 180 CONTINUE
346 ioffa = ioffa + lda
347 190 CONTINUE
348 END IF
349 jj = jj + jb
350 END IF
351*
352 IF( myrow.EQ.iarow )
353 $ ii = ii + jb
354 iarow = mod( iarow+1, nprow )
355 iacol = mod( iacol+1, npcol )
356*
357* Loop over remaining block of columns
358*
359 DO 260 j = jn+1, ja+n-1, desca( nb_ )
360 jb = min( ja+n-j, desca( nb_ ) )
361*
362 IF( mycol.EQ.iacol ) THEN
363 IF( myrow.EQ.iarow ) THEN
364 IF( udiag ) THEN
365 DO 210 ll = jj, jj + jb -1
366 DO 200 kk = ii+ll-jj+1, iia+mp-1
367 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
368 200 CONTINUE
369 ioffa = ioffa + lda
370 210 CONTINUE
371 ELSE
372 DO 230 ll = jj, jj + jb -1
373 DO 220 kk = ii+ll-jj, iia+mp-1
374 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
375 220 CONTINUE
376 ioffa = ioffa + lda
377 230 CONTINUE
378 END IF
379 ELSE
380 DO 250 ll = jj, jj + jb -1
381 DO 240 kk = ii, iia+mp-1
382 VALUE = max( VALUE, abs( a( ioffa+kk ) ) )
383 240 CONTINUE
384 ioffa = ioffa + lda
385 250 CONTINUE
386 END IF
387 jj = jj + jb
388 END IF
389*
390 IF( myrow.EQ.iarow )
391 $ ii = ii + jb
392 iarow = mod( iarow+1, nprow )
393 iacol = mod( iacol+1, npcol )
394*
395 260 CONTINUE
396*
397 END IF
398*
399* Gather the intermediate results to process (0,0).
400*
401 CALL dgamx2d( ictxt, 'All', ' ', 1, 1, VALUE, 1, kk, ll, -1,
402 $ 0, 0 )
403*
404************************************************************************
405* one norm
406*
407 ELSE IF( lsame( norm, 'O' ) .OR. norm.EQ.'1' ) THEN
408*
409 VALUE = zero
410*
411 IF( lsame( uplo, 'U' ) ) THEN
412*
413* Upper triangular matrix
414*
415 ii = iia
416 jj = jja
417 jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
418 jb = jn-ja+1
419*
420 IF( mycol.EQ.iacol ) THEN
421 IF( myrow.EQ.iarow ) THEN
422 IF( udiag ) THEN
423 DO 280 ll = jj, jj + jb -1
424 sum = zero
425 DO 270 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
426 sum = sum + abs( a( ioffa+kk ) )
427 270 CONTINUE
428* Unit diagonal entry
429 kk = ii+ll-jj
430 IF (kk <= iia+mp-1) THEN
431 sum = sum + one
432 ENDIF
433 ioffa = ioffa + lda
434 work( ll-jja+1 ) = sum
435 280 CONTINUE
436 ELSE
437 DO 300 ll = jj, jj + jb -1
438 sum = zero
439 DO 290 kk = iia, min( ii+ll-jj, iia+mp-1 )
440 sum = sum + abs( a( ioffa+kk ) )
441 290 CONTINUE
442 ioffa = ioffa + lda
443 work( ll-jja+1 ) = sum
444 300 CONTINUE
445 END IF
446 ELSE
447 DO 320 ll = jj, jj + jb -1
448 sum = zero
449 DO 310 kk = iia, min( ii-1, iia+mp-1 )
450 sum = sum + abs( a( ioffa+kk ) )
451 310 CONTINUE
452 ioffa = ioffa + lda
453 work( ll-jja+1 ) = sum
454 320 CONTINUE
455 END IF
456 jj = jj + jb
457 END IF
458*
459 IF( myrow.EQ.iarow )
460 $ ii = ii + jb
461 iarow = mod( iarow+1, nprow )
462 iacol = mod( iacol+1, npcol )
463*
464* Loop over remaining block of columns
465*
466 DO 390 j = jn+1, ja+n-1, desca( nb_ )
467 jb = min( ja+n-j, desca( nb_ ) )
468*
469 IF( mycol.EQ.iacol ) THEN
470 IF( myrow.EQ.iarow ) THEN
471 IF( udiag ) THEN
472 DO 340 ll = jj, jj + jb -1
473 sum = zero
474 DO 330 kk = iia, min( ii+ll-jj-1, iia+mp-1 )
475 sum = sum + abs( a( ioffa+kk ) )
476 330 CONTINUE
477* Unit diagonal entry
478 kk = ii+ll-jj
479 IF (kk <= iia+mp-1) THEN
480 sum = sum + one
481 ENDIF
482 ioffa = ioffa + lda
483 work( ll-jja+1 ) = sum
484 340 CONTINUE
485 ELSE
486 DO 360 ll = jj, jj + jb -1
487 sum = zero
488 DO 350 kk = iia, min( ii+ll-jj, iia+mp-1 )
489 sum = sum + abs( a( ioffa+kk ) )
490 350 CONTINUE
491 ioffa = ioffa + lda
492 work( ll-jja+1 ) = sum
493 360 CONTINUE
494 END IF
495 ELSE
496 DO 380 ll = jj, jj + jb -1
497 sum = zero
498 DO 370 kk = iia, min( ii-1, iia+mp-1 )
499 sum = sum + abs( a( ioffa+kk ) )
500 370 CONTINUE
501 ioffa = ioffa + lda
502 work( ll-jja+1 ) = sum
503 380 CONTINUE
504 END IF
505 jj = jj + jb
506 END IF
507*
508 IF( myrow.EQ.iarow )
509 $ ii = ii + jb
510 iarow = mod( iarow+1, nprow )
511 iacol = mod( iacol+1, npcol )
512*
513 390 CONTINUE
514*
515 ELSE
516*
517* Lower triangular matrix
518*
519 ii = iia
520 jj = jja
521 jn = min( iceil( ja, desca( nb_ ) ) * desca( nb_ ), ja+n-1 )
522 jb = jn-ja+1
523*
524 IF( mycol.EQ.iacol ) THEN
525 IF( myrow.EQ.iarow ) THEN
526 IF( udiag ) THEN
527 DO 410 ll = jj, jj + jb -1
528 sum = one
529 DO 400 kk = ii+ll-jj+1, iia+mp-1
530 sum = sum + abs( a( ioffa+kk ) )
531 400 CONTINUE
532 ioffa = ioffa + lda
533 work( ll-jja+1 ) = sum
534 410 CONTINUE
535 ELSE
536 DO 430 ll = jj, jj + jb -1
537 sum = zero
538 DO 420 kk = ii+ll-jj, iia+mp-1
539 sum = sum + abs( a( ioffa+kk ) )
540 420 CONTINUE
541 ioffa = ioffa + lda
542 work( ll-jja+1 ) = sum
543 430 CONTINUE
544 END IF
545 ELSE
546 DO 450 ll = jj, jj + jb -1
547 sum = zero
548 DO 440 kk = ii, iia+mp-1
549 sum = sum + abs( a( ioffa+kk ) )
550 440 CONTINUE
551 ioffa = ioffa + lda
552 work( ll-jja+1 ) = sum
553 450 CONTINUE
554 END IF
555 jj = jj + jb
556 END IF
557*
558 IF( myrow.EQ.iarow )
559 $ ii = ii + jb
560 iarow = mod( iarow+1, nprow )
561 iacol = mod( iacol+1, npcol )
562*
563* Loop over remaining block of columns
564*
565 DO 520 j = jn+1, ja+n-1, desca( nb_ )
566 jb = min( ja+n-j, desca( nb_ ) )
567*
568 IF( mycol.EQ.iacol ) THEN
569 IF( myrow.EQ.iarow ) THEN
570 IF( udiag ) THEN
571 DO 470 ll = jj, jj + jb -1
572 sum = one
573 DO 460 kk = ii+ll-jj+1, iia+mp-1
574 sum = sum + abs( a( ioffa+kk ) )
575 460 CONTINUE
576 ioffa = ioffa + lda
577 work( ll-jja+1 ) = sum
578 470 CONTINUE
579 ELSE
580 DO 490 ll = jj, jj + jb -1
581 sum = zero
582 DO 480 kk = ii+ll-jj, iia+mp-1
583 sum = sum + abs( a( ioffa+kk ) )
584 480 CONTINUE
585 ioffa = ioffa + lda
586 work( ll-jja+1 ) = sum
587 490 CONTINUE
588 END IF
589 ELSE
590 DO 510 ll = jj, jj + jb -1
591 sum = zero
592 DO 500 kk = ii, iia+mp-1
593 sum = sum + abs( a( ioffa+kk ) )
594 500 CONTINUE
595 ioffa = ioffa + lda
596 work( ll-jja+1 ) = sum
597 510 CONTINUE
598 END IF
599 jj = jj + jb
600 END IF
601*
602 IF( myrow.EQ.iarow )
603 $ ii = ii + jb
604 iarow = mod( iarow+1, nprow )
605 iacol = mod( iacol+1, npcol )
606*
607 520 CONTINUE
608*
609 END IF
610*
611* Find sum of global matrix columns and store on row 0 of
612* process grid
613*
614 CALL dgsum2d( ictxt, 'Columnwise', ' ', 1, nq, work, 1,
615 $ 0, mycol )
616*
617* Find maximum sum of columns for 1-norm
618*
619 IF( myrow.EQ.0 ) THEN
620 IF( nq.GT.0 ) THEN
621 VALUE = work( idamax( nq, work, 1 ) )
622 ELSE
623 VALUE = zero
624 END IF
625 CALL dgamx2d( ictxt, 'Rowwise', ' ', 1, 1, VALUE, 1, KK, LL,
626 $ -1, 0, 0 )
627 END IF
628*
629************************************************************************
630* infinity norm
631*
632 ELSE IF( LSAME( NORM, 'i' ) ) THEN
633*
634 IF( LSAME( UPLO, 'u' ) ) THEN
635 DO 540 KK = IIA, IIA+MP-1
636 WORK( KK ) = ZERO
637 540 CONTINUE
638 ELSE
639 DO 570 KK = IIA, IIA+MP-1
640 WORK( KK ) = ZERO
641 570 CONTINUE
642 END IF
643*
644 IF( LSAME( UPLO, 'u' ) ) THEN
645*
646* Upper triangular matrix
647*
648 II = IIA
649 JJ = JJA
650 JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
651 JB = JN-JA+1
652*
653.EQ. IF( MYCOLIACOL ) THEN
654.EQ. IF( MYROWIAROW ) THEN
655 IF( UDIAG ) THEN
656 DO 590 LL = JJ, JJ + JB -1
657 DO 580 KK = IIA, MIN( II+LL-JJ-1, IIA+MP-1 )
658 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
659 $ ABS( A( IOFFA+KK ) )
660 580 CONTINUE
661* Unit diagonal entry
662 KK = II+LL-JJ
663 IF (KK <= IIA+MP-1) THEN
664 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) + ONE
665 ENDIF
666 IOFFA = IOFFA + LDA
667 590 CONTINUE
668 ELSE
669 DO 610 LL = JJ, JJ + JB -1
670 DO 600 KK = IIA, MIN( II+LL-JJ, IIA+MP-1 )
671 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
672 $ ABS( A( IOFFA+KK ) )
673 600 CONTINUE
674 IOFFA = IOFFA + LDA
675 610 CONTINUE
676 END IF
677 ELSE
678 DO 630 LL = JJ, JJ + JB -1
679 DO 620 KK = IIA, MIN( II-1, IIA+MP-1 )
680 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
681 $ ABS( A( IOFFA+KK ) )
682 620 CONTINUE
683 IOFFA = IOFFA + LDA
684 630 CONTINUE
685 END IF
686 JJ = JJ + JB
687 END IF
688*
689.EQ. IF( MYROWIAROW )
690 $ II = II + JB
691 IAROW = MOD( IAROW+1, NPROW )
692 IACOL = MOD( IACOL+1, NPCOL )
693*
694* Loop over remaining block of columns
695*
696 DO 700 J = JN+1, JA+N-1, DESCA( NB_ )
697 JB = MIN( JA+N-J, DESCA( NB_ ) )
698*
699.EQ. IF( MYCOLIACOL ) THEN
700.EQ. IF( MYROWIAROW ) THEN
701 IF( UDIAG ) THEN
702 DO 650 LL = JJ, JJ + JB -1
703 DO 640 KK = IIA, MIN( II+LL-JJ-1, IIA+MP-1 )
704 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
705 $ ABS( A( IOFFA+KK ) )
706 640 CONTINUE
707* Unit diagonal entry
708 KK = II+LL-JJ
709 IF (KK <= IIA+MP-1) THEN
710 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) + ONE
711 ENDIF
712 IOFFA = IOFFA + LDA
713 650 CONTINUE
714 ELSE
715 DO 670 LL = JJ, JJ + JB -1
716 DO 660 KK = IIA, MIN( II+LL-JJ, IIA+MP-1 )
717 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
718 $ ABS( A( IOFFA+KK ) )
719 660 CONTINUE
720 IOFFA = IOFFA + LDA
721 670 CONTINUE
722 END IF
723 ELSE
724 DO 690 LL = JJ, JJ + JB -1
725 DO 680 KK = IIA, MIN( II-1, IIA+MP-1 )
726 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
727 $ ABS( A( IOFFA+KK ) )
728 680 CONTINUE
729 IOFFA = IOFFA + LDA
730 690 CONTINUE
731 END IF
732 JJ = JJ + JB
733 END IF
734*
735.EQ. IF( MYROWIAROW )
736 $ II = II + JB
737 IAROW = MOD( IAROW+1, NPROW )
738 IACOL = MOD( IACOL+1, NPCOL )
739*
740 700 CONTINUE
741*
742 ELSE
743*
744* Lower triangular matrix
745*
746 II = IIA
747 JJ = JJA
748 JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
749 JB = JN-JA+1
750*
751.EQ. IF( MYCOLIACOL ) THEN
752.EQ. IF( MYROWIAROW ) THEN
753 IF( UDIAG ) THEN
754 DO 720 LL = JJ, JJ + JB -1
755* Unit diagonal entry
756 KK = II+LL-JJ
757 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) + ONE
758 DO 710 KK = II+LL-JJ+1, IIA+MP-1
759 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
760 $ ABS( A( IOFFA+KK ) )
761 710 CONTINUE
762 IOFFA = IOFFA + LDA
763 720 CONTINUE
764 ELSE
765 DO 740 LL = JJ, JJ + JB -1
766 DO 730 KK = II+LL-JJ, IIA+MP-1
767 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
768 $ ABS( A( IOFFA+KK ) )
769 730 CONTINUE
770 IOFFA = IOFFA + LDA
771 740 CONTINUE
772 END IF
773 ELSE
774 DO 760 LL = JJ, JJ + JB -1
775 DO 750 KK = II, IIA+MP-1
776 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
777 $ ABS( A( IOFFA+KK ) )
778 750 CONTINUE
779 IOFFA = IOFFA + LDA
780 760 CONTINUE
781 END IF
782 JJ = JJ + JB
783 END IF
784*
785.EQ. IF( MYROWIAROW )
786 $ II = II + JB
787 IAROW = MOD( IAROW+1, NPROW )
788 IACOL = MOD( IACOL+1, NPCOL )
789*
790* Loop over remaining block of columns
791*
792 DO 830 J = JN+1, JA+N-1, DESCA( NB_ )
793 JB = MIN( JA+N-J, DESCA( NB_ ) )
794*
795.EQ. IF( MYCOLIACOL ) THEN
796.EQ. IF( MYROWIAROW ) THEN
797 IF( UDIAG ) THEN
798 DO 780 LL = JJ, JJ + JB -1
799* Unit diagonal entry
800 KK = II+LL-JJ
801 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) + ONE
802 DO 770 KK = II+LL-JJ+1, IIA+MP-1
803 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
804 $ ABS( A( IOFFA+KK ) )
805 770 CONTINUE
806 IOFFA = IOFFA + LDA
807 780 CONTINUE
808 ELSE
809 DO 800 LL = JJ, JJ + JB -1
810 DO 790 KK = II+LL-JJ, IIA+MP-1
811 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
812 $ ABS( A( IOFFA+KK ) )
813 790 CONTINUE
814 IOFFA = IOFFA + LDA
815 800 CONTINUE
816 END IF
817 ELSE
818 DO 820 LL = JJ, JJ + JB -1
819 DO 810 KK = II, IIA+MP-1
820 WORK( KK-IIA+1 ) = WORK( KK-IIA+1 ) +
821 $ ABS( A( IOFFA+KK ) )
822 810 CONTINUE
823 IOFFA = IOFFA + LDA
824 820 CONTINUE
825 END IF
826 JJ = JJ + JB
827 END IF
828*
829.EQ. IF( MYROWIAROW )
830 $ II = II + JB
831 IAROW = MOD( IAROW+1, NPROW )
832 IACOL = MOD( IACOL+1, NPCOL )
833*
834 830 CONTINUE
835*
836 END IF
837*
838* Find sum of global matrix rows and store on column 0 of
839* process grid
840*
841 CALL DGSUM2D( ICTXT, 'rowwise', ' ', MP, 1, WORK, MAX( 1, MP ),
842 $ MYROW, 0 )
843*
844* Find maximum sum of rows for Infinity-norm
845*
846.EQ. IF( MYCOL0 ) THEN
847.GT. IF( MP0 ) THEN
848 VALUE = WORK( IDAMAX( MP, WORK, 1 ) )
849 ELSE
850 VALUE = ZERO
851 END IF
852 CALL DGAMX2D( ICTXT, 'columnwise', ' ', 1, 1, VALUE, 1, KK,
853 $ LL, -1, 0, 0 )
854 END IF
855*
856************************************************************************
857* Frobenius norm
858* SSQ(1) is scale
859* SSQ(2) is sum-of-squares
860*
861 ELSE IF( LSAME( NORM, 'f.OR.' ) LSAME( NORM, 'e' ) ) THEN
862*
863 IF( UDIAG ) THEN
864 SSQ(1) = ONE
865 SSQ(2) = DBLE( MIN( M, N ) ) / DBLE( NPROW*NPCOL )
866 ELSE
867 SSQ(1) = ZERO
868 SSQ(2) = ONE
869 END IF
870*
871 IF( LSAME( UPLO, 'u' ) ) THEN
872*
873* ***********************
874* Upper triangular matrix
875*
876 II = IIA
877 JJ = JJA
878 JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
879 JB = JN-JA+1
880*
881* First block column of sub-matrix.
882*
883.EQ. IF( MYCOLIACOL ) THEN
884.EQ. IF( MYROWIAROW ) THEN
885* This process has part of current block column,
886* including diagonal block.
887 IF( UDIAG ) THEN
888 DO 840 LL = JJ, JJ + JB -1
889 COLSSQ(1) = ZERO
890 COLSSQ(2) = ONE
891 CALL ZLASSQ( MIN( II+LL-JJ-1, IIA+MP-1 )-IIA+1,
892 $ A( IIA+IOFFA ), 1,
893 $ COLSSQ(1), COLSSQ(2) )
894 CALL DCOMBSSQ( SSQ, COLSSQ )
895 IOFFA = IOFFA + LDA
896 840 CONTINUE
897 ELSE
898 DO 850 LL = JJ, JJ + JB -1
899 COLSSQ(1) = ZERO
900 COLSSQ(2) = ONE
901 CALL ZLASSQ( MIN( II+LL-JJ, IIA+MP-1 )-IIA+1,
902 $ A( IIA+IOFFA ), 1,
903 $ COLSSQ(1), COLSSQ(2) )
904 CALL DCOMBSSQ( SSQ, COLSSQ )
905 IOFFA = IOFFA + LDA
906 850 CONTINUE
907 END IF
908 ELSE
909* This rank has part of current block column,
910* but not diagonal block.
911* It seems this lassq will be length 0, since ii = iia.
912 DO 860 LL = JJ, JJ + JB -1
913 COLSSQ(1) = ZERO
914 COLSSQ(2) = ONE
915 CALL ZLASSQ( MIN( II-1, IIA+MP-1 )-IIA+1,
916 $ A( IIA+IOFFA ), 1,
917 $ COLSSQ(1), COLSSQ(2) )
918 CALL DCOMBSSQ( SSQ, COLSSQ )
919 IOFFA = IOFFA + LDA
920 860 CONTINUE
921 END IF
922 JJ = JJ + JB
923 END IF
924*
925* If this process has part of current block row, advance ii,
926* then advance iarow, iacol to next diagonal block.
927*
928.EQ. IF( MYROWIAROW )
929 $ II = II + JB
930 IAROW = MOD( IAROW+1, NPROW )
931 IACOL = MOD( IACOL+1, NPCOL )
932*
933* Loop over remaining block columns
934*
935 DO 900 J = JN+1, JA+N-1, DESCA( NB_ )
936 JB = MIN( JA+N-J, DESCA( NB_ ) )
937*
938.EQ. IF( MYCOLIACOL ) THEN
939.EQ. IF( MYROWIAROW ) THEN
940 IF( UDIAG ) THEN
941 DO 870 LL = JJ, JJ + JB -1
942 COLSSQ(1) = ZERO
943 COLSSQ(2) = ONE
944 CALL ZLASSQ( MIN(II+LL-JJ-1, IIA+MP-1)-IIA+1,
945 $ A( IIA+IOFFA ), 1,
946 $ COLSSQ(1), COLSSQ(2) )
947 CALL DCOMBSSQ( SSQ, COLSSQ )
948 IOFFA = IOFFA + LDA
949 870 CONTINUE
950 ELSE
951 DO 880 LL = JJ, JJ + JB -1
952 COLSSQ(1) = ZERO
953 COLSSQ(2) = ONE
954 CALL ZLASSQ( MIN( II+LL-JJ, IIA+MP-1 )-IIA+1,
955 $ A( IIA+IOFFA ), 1,
956 $ COLSSQ(1), COLSSQ(2) )
957 CALL DCOMBSSQ( SSQ, COLSSQ )
958 IOFFA = IOFFA + LDA
959 880 CONTINUE
960 END IF
961 ELSE
962 DO 890 LL = JJ, JJ + JB -1
963 COLSSQ(1) = ZERO
964 COLSSQ(2) = ONE
965 CALL ZLASSQ( MIN( II-1, IIA+MP-1 )-IIA+1,
966 $ A( IIA+IOFFA ), 1,
967 $ COLSSQ(1), COLSSQ(2) )
968 CALL DCOMBSSQ( SSQ, COLSSQ )
969 IOFFA = IOFFA + LDA
970 890 CONTINUE
971 END IF
972 JJ = JJ + JB
973 END IF
974*
975.EQ. IF( MYROWIAROW )
976 $ II = II + JB
977 IAROW = MOD( IAROW+1, NPROW )
978 IACOL = MOD( IACOL+1, NPCOL )
979*
980 900 CONTINUE
981*
982 ELSE
983*
984* ***********************
985* Lower triangular matrix
986*
987 II = IIA
988 JJ = JJA
989 JN = MIN( ICEIL( JA, DESCA( NB_ ) ) * DESCA( NB_ ), JA+N-1 )
990 JB = JN-JA+1
991*
992.EQ. IF( MYCOLIACOL ) THEN
993.EQ. IF( MYROWIAROW ) THEN
994 IF( UDIAG ) THEN
995 DO 910 LL = JJ, JJ + JB -1
996 COLSSQ(1) = ZERO
997 COLSSQ(2) = ONE
998 CALL ZLASSQ( IIA+MP-(II+LL-JJ+1),
999 $ A( II+LL-JJ+1+IOFFA ), 1,
1000 $ COLSSQ(1), COLSSQ(2) )
1001 CALL DCOMBSSQ( SSQ, COLSSQ )
1002 IOFFA = IOFFA + LDA
1003 910 CONTINUE
1004 ELSE
1005 DO 920 LL = JJ, JJ + JB -1
1006 COLSSQ(1) = ZERO
1007 COLSSQ(2) = ONE
1008 CALL ZLASSQ( IIA+MP-(II+LL-JJ),
1009 $ A( II+LL-JJ+IOFFA ), 1,
1010 $ COLSSQ(1), COLSSQ(2) )
1011 CALL DCOMBSSQ( SSQ, COLSSQ )
1012 IOFFA = IOFFA + LDA
1013 920 CONTINUE
1014 END IF
1015 ELSE
1016 DO 930 LL = JJ, JJ + JB -1
1017 COLSSQ(1) = ZERO
1018 COLSSQ(2) = ONE
1019 CALL ZLASSQ( IIA+MP-II, A( II+IOFFA ), 1,
1020 $ COLSSQ(1), COLSSQ(2) )
1021 CALL DCOMBSSQ( SSQ, COLSSQ )
1022 IOFFA = IOFFA + LDA
1023 930 CONTINUE
1024 END IF
1025 JJ = JJ + JB
1026 END IF
1027*
1028.EQ. IF( MYROWIAROW )
1029 $ II = II + JB
1030 IAROW = MOD( IAROW+1, NPROW )
1031 IACOL = MOD( IACOL+1, NPCOL )
1032*
1033* Loop over remaining block of columns
1034*
1035 DO 970 J = JN+1, JA+N-1, DESCA( NB_ )
1036 JB = MIN( JA+N-J, DESCA( NB_ ) )
1037*
1038.EQ. IF( MYCOLIACOL ) THEN
1039.EQ. IF( MYROWIAROW ) THEN
1040 IF( UDIAG ) THEN
1041 DO 940 LL = JJ, JJ + JB -1
1042 COLSSQ(1) = ZERO
1043 COLSSQ(2) = ONE
1044 CALL ZLASSQ( IIA+MP-(II+LL-JJ+1),
1045 $ A( II+LL-JJ+1+IOFFA ), 1,
1046 $ COLSSQ(1), COLSSQ(2) )
1047 CALL DCOMBSSQ( SSQ, COLSSQ )
1048 IOFFA = IOFFA + LDA
1049 940 CONTINUE
1050 ELSE
1051 DO 950 LL = JJ, JJ + JB -1
1052 COLSSQ(1) = ZERO
1053 COLSSQ(2) = ONE
1054 CALL ZLASSQ( IIA+MP-(II+LL-JJ),
1055 $ A( II+LL-JJ+IOFFA ), 1,
1056 $ COLSSQ(1), COLSSQ(2) )
1057 CALL DCOMBSSQ( SSQ, COLSSQ )
1058 IOFFA = IOFFA + LDA
1059 950 CONTINUE
1060 END IF
1061 ELSE
1062 DO 960 LL = JJ, JJ + JB -1
1063 COLSSQ(1) = ZERO
1064 COLSSQ(2) = ONE
1065 CALL ZLASSQ( IIA+MP-II, A( II+IOFFA ), 1,
1066 $ COLSSQ(1), COLSSQ(2) )
1067 CALL DCOMBSSQ( SSQ, COLSSQ )
1068 IOFFA = IOFFA + LDA
1069 960 CONTINUE
1070 END IF
1071 JJ = JJ + JB
1072 END IF
1073*
1074.EQ. IF( MYROWIAROW )
1075 $ II = II + JB
1076 IAROW = MOD( IAROW+1, NPROW )
1077 IACOL = MOD( IACOL+1, NPCOL )
1078*
1079 970 CONTINUE
1080*
1081 END IF
1082*
1083* ***********************
1084* Perform the global scaled sum
1085*
1086 CALL PDTREECOMB( ICTXT, 'all', 2, SSQ, 0, 0, DCOMBSSQ )
1087 VALUE = SSQ( 1 ) * SQRT( SSQ( 2 ) )
1088*
1089 END IF
1090*
1091* Broadcast the result to every process in the grid.
1092*
1093.EQ..AND..EQ. IF( MYROW0 MYCOL0 ) THEN
1094 CALL DGEBS2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1 )
1095 ELSE
1096 CALL DGEBR2D( ICTXT, 'all', ' ', 1, 1, VALUE, 1, 0, 0 )
1097 END IF
1098*
1099 PZLANTR = VALUE
1100*
1101 RETURN
1102*
1103* End of PZLANTR
1104*
1105 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
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
dcombssq
subroutine dcombssq(v1, v2)
Definition pdtreecomb.f:259
pdtreecomb
subroutine pdtreecomb(ictxt, scope, n, mine, rdest0, cdest0, subptr)
Definition pdtreecomb.f:3
pzlantr
double precision function pzlantr(norm, uplo, diag, m, n, a, ia, ja, desca, work)
Definition pzlantr.f:3