OpenRadioss 2025.1.11
OpenRadioss project
Loading...
Searching...
No Matches
slansf.f
Go to the documentation of this file.
1*> \brief \b SLANSF
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download SLANSF + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slansf.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slansf.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slansf.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* REAL FUNCTION SLANSF( NORM, TRANSR, UPLO, N, A, WORK )
22*
23* .. Scalar Arguments ..
24* CHARACTER NORM, TRANSR, UPLO
25* INTEGER N
26* ..
27* .. Array Arguments ..
28* REAL A( 0: * ), WORK( 0: * )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> SLANSF returns the value of the one norm, or the Frobenius norm, or
38*> the infinity norm, or the element of largest absolute value of a
39*> real symmetric matrix A in RFP format.
40*> \endverbatim
41*>
42*> \return SLANSF
43*> \verbatim
44*>
45*> SLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
46*> (
47*> ( norm1(A), NORM = '1', 'O' or 'o'
48*> (
49*> ( normI(A), NORM = 'I' or 'i'
50*> (
51*> ( normF(A), NORM = 'F', 'f', 'E' or 'e'
52*>
53*> where norm1 denotes the one norm of a matrix (maximum column sum),
54*> normI denotes the infinity norm of a matrix (maximum row sum) and
55*> normF denotes the Frobenius norm of a matrix (square root of sum of
56*> squares). Note that max(abs(A(i,j))) is not a matrix norm.
57*> \endverbatim
58*
59* Arguments:
60* ==========
61*
62*> \param[in] NORM
63*> \verbatim
64*> NORM is CHARACTER*1
65*> Specifies the value to be returned in SLANSF as described
66*> above.
67*> \endverbatim
68*>
69*> \param[in] TRANSR
70*> \verbatim
71*> TRANSR is CHARACTER*1
72*> Specifies whether the RFP format of A is normal or
73*> transposed format.
74*> = 'N': RFP format is Normal;
75*> = 'T': RFP format is Transpose.
76*> \endverbatim
77*>
78*> \param[in] UPLO
79*> \verbatim
80*> UPLO is CHARACTER*1
81*> On entry, UPLO specifies whether the RFP matrix A came from
82*> an upper or lower triangular matrix as follows:
83*> = 'U': RFP A came from an upper triangular matrix;
84*> = 'L': RFP A came from a lower triangular matrix.
85*> \endverbatim
86*>
87*> \param[in] N
88*> \verbatim
89*> N is INTEGER
90*> The order of the matrix A. N >= 0. When N = 0, SLANSF is
91*> set to zero.
92*> \endverbatim
93*>
94*> \param[in] A
95*> \verbatim
96*> A is REAL array, dimension ( N*(N+1)/2 );
97*> On entry, the upper (if UPLO = 'U') or lower (if UPLO = 'L')
98*> part of the symmetric matrix A stored in RFP format. See the
99*> "Notes" below for more details.
100*> Unchanged on exit.
101*> \endverbatim
102*>
103*> \param[out] WORK
104*> \verbatim
105*> WORK is REAL array, dimension (MAX(1,LWORK)),
106*> where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
107*> WORK is not referenced.
108*> \endverbatim
109*
110* Authors:
111* ========
112*
113*> \author Univ. of Tennessee
114*> \author Univ. of California Berkeley
115*> \author Univ. of Colorado Denver
116*> \author NAG Ltd.
117*
118*> \ingroup realOTHERcomputational
119*
120*> \par Further Details:
121* =====================
122*>
123*> \verbatim
124*>
125*> We first consider Rectangular Full Packed (RFP) Format when N is
126*> even. We give an example where N = 6.
127*>
128*> AP is Upper AP is Lower
129*>
130*> 00 01 02 03 04 05 00
131*> 11 12 13 14 15 10 11
132*> 22 23 24 25 20 21 22
133*> 33 34 35 30 31 32 33
134*> 44 45 40 41 42 43 44
135*> 55 50 51 52 53 54 55
136*>
137*>
138*> Let TRANSR = 'N'. RFP holds AP as follows:
139*> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
140*> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
141*> the transpose of the first three columns of AP upper.
142*> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
143*> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
144*> the transpose of the last three columns of AP lower.
145*> This covers the case N even and TRANSR = 'N'.
146*>
147*> RFP A RFP A
148*>
149*> 03 04 05 33 43 53
150*> 13 14 15 00 44 54
151*> 23 24 25 10 11 55
152*> 33 34 35 20 21 22
153*> 00 44 45 30 31 32
154*> 01 11 55 40 41 42
155*> 02 12 22 50 51 52
156*>
157*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
158*> transpose of RFP A above. One therefore gets:
159*>
160*>
161*> RFP A RFP A
162*>
163*> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
164*> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
165*> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
166*>
167*>
168*> We then consider Rectangular Full Packed (RFP) Format when N is
169*> odd. We give an example where N = 5.
170*>
171*> AP is Upper AP is Lower
172*>
173*> 00 01 02 03 04 00
174*> 11 12 13 14 10 11
175*> 22 23 24 20 21 22
176*> 33 34 30 31 32 33
177*> 44 40 41 42 43 44
178*>
179*>
180*> Let TRANSR = 'N'. RFP holds AP as follows:
181*> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
182*> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
183*> the transpose of the first two columns of AP upper.
184*> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
185*> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
186*> the transpose of the last two columns of AP lower.
187*> This covers the case N odd and TRANSR = 'N'.
188*>
189*> RFP A RFP A
190*>
191*> 02 03 04 00 33 43
192*> 12 13 14 10 11 44
193*> 22 23 24 20 21 22
194*> 00 33 34 30 31 32
195*> 01 11 44 40 41 42
196*>
197*> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
198*> transpose of RFP A above. One therefore gets:
199*>
200*> RFP A RFP A
201*>
202*> 02 12 22 00 01 00 10 20 30 40 50
203*> 03 13 23 33 11 33 11 21 31 41 51
204*> 04 14 24 34 44 43 44 22 32 42 52
205*> \endverbatim
206*
207* =====================================================================
208 REAL function slansf( norm, transr, uplo, n, a, work )
209*
210* -- LAPACK computational routine --
211* -- LAPACK is a software package provided by Univ. of Tennessee, --
212* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
213*
214* .. Scalar Arguments ..
215 CHARACTER norm, TRANSR, uplo
216 INTEGER n
217* ..
218* .. Array Arguments ..
219 REAL a( 0: * ), work( 0: * )
220* ..
221*
222* =====================================================================
223*
224* ..
225* .. Parameters ..
226 REAL one, zero
227 parameter( one = 1.0e+0, zero = 0.0e+0 )
228* ..
229* .. Local Scalars ..
230 INTEGER i, j, ifm, ilu, noe, n1, k, l, lda
231 REAL scale, s, VALUE, aa, temp
232* ..
233* .. External Functions ..
234 LOGICAL lsame, sisnan
235 EXTERNAL lsame, sisnan
236* ..
237* .. External Subroutines ..
238 EXTERNAL slassq
239* ..
240* .. Intrinsic Functions ..
241 INTRINSIC abs, sqrt
242* ..
243* .. Executable Statements ..
244*
245 IF( n.EQ.0 ) THEN
246 slansf = zero
247 RETURN
248 ELSE IF( n.EQ.1 ) THEN
249 slansf = abs( a(0) )
250 RETURN
251 END IF
252*
253* set noe = 1 if n is odd. if n is even set noe=0
254*
255 noe = 1
256 IF( mod( n, 2 ).EQ.0 )
257 $ noe = 0
258*
259* set ifm = 0 when form='T or 't' and 1 otherwise
260*
261 ifm = 1
262 IF( lsame( transr, 'T' ) )
263 $ ifm = 0
264*
265* set ilu = 0 when uplo='U or 'u' and 1 otherwise
266*
267 ilu = 1
268 IF( lsame( uplo, 'U' ) )
269 $ ilu = 0
270*
271* set lda = (n+1)/2 when ifm = 0
272* set lda = n when ifm = 1 and noe = 1
273* set lda = n+1 when ifm = 1 and noe = 0
274*
275 IF( ifm.EQ.1 ) THEN
276 IF( noe.EQ.1 ) THEN
277 lda = n
278 ELSE
279* noe=0
280 lda = n + 1
281 END IF
282 ELSE
283* ifm=0
284 lda = ( n+1 ) / 2
285 END IF
286*
287 IF( lsame( norm, 'M' ) ) THEN
288*
289* Find max(abs(A(i,j))).
290*
291 k = ( n+1 ) / 2
292 VALUE = zero
293 IF( noe.EQ.1 ) THEN
294* n is odd
295 IF( ifm.EQ.1 ) THEN
296* A is n by k
297 DO j = 0, k - 1
298 DO i = 0, n - 1
299 temp = abs( a( i+j*lda ) )
300 IF( VALUE .LT. temp .OR. sisnan( temp ) )
301 $ VALUE = temp
302 END DO
303 END DO
304 ELSE
305* xpose case; A is k by n
306 DO j = 0, n - 1
307 DO i = 0, k - 1
308 temp = abs( a( i+j*lda ) )
309 IF( VALUE .LT. temp .OR. sisnan( temp ) )
310 $ VALUE = temp
311 END DO
312 END DO
313 END IF
314 ELSE
315* n is even
316 IF( ifm.EQ.1 ) THEN
317* A is n+1 by k
318 DO j = 0, k - 1
319 DO i = 0, n
320 temp = abs( a( i+j*lda ) )
321 IF( VALUE .LT. temp .OR. sisnan( temp ) )
322 $ VALUE = temp
323 END DO
324 END DO
325 ELSE
326* xpose case; A is k by n+1
327 DO j = 0, n
328 DO i = 0, k - 1
329 temp = abs( a( i+j*lda ) )
330 IF( VALUE .LT. temp .OR. sisnan( temp ) )
331 $ VALUE = temp
332 END DO
333 END DO
334 END IF
335 END IF
336 ELSE IF( ( lsame( norm, 'I' ) ) .OR. ( lsame( norm, 'o.OR.' ) )
337.EQ. $ ( NORM'1' ) ) THEN
338*
339* Find normI(A) ( = norm1(A), since A is symmetric).
340*
341.EQ. IF( IFM1 ) THEN
342 K = N / 2
343.EQ. IF( NOE1 ) THEN
344* n is odd
345.EQ. IF( ILU0 ) THEN
346 DO I = 0, K - 1
347 WORK( I ) = ZERO
348 END DO
349 DO J = 0, K
350 S = ZERO
351 DO I = 0, K + J - 1
352 AA = ABS( A( I+J*LDA ) )
353* -> A(i,j+k)
354 S = S + AA
355 WORK( I ) = WORK( I ) + AA
356 END DO
357 AA = ABS( A( I+J*LDA ) )
358* -> A(j+k,j+k)
359 WORK( J+K ) = S + AA
360.EQ. IF( IK+K )
361 $ GO TO 10
362 I = I + 1
363 AA = ABS( A( I+J*LDA ) )
364* -> A(j,j)
365 WORK( J ) = WORK( J ) + AA
366 S = ZERO
367 DO L = J + 1, K - 1
368 I = I + 1
369 AA = ABS( A( I+J*LDA ) )
370* -> A(l,j)
371 S = S + AA
372 WORK( L ) = WORK( L ) + AA
373 END DO
374 WORK( J ) = WORK( J ) + S
375 END DO
376 10 CONTINUE
377 VALUE = WORK( 0 )
378 DO I = 1, N-1
379 TEMP = WORK( I )
380.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
381 $ VALUE = TEMP
382 END DO
383 ELSE
384* ilu = 1
385 K = K + 1
386* k=(n+1)/2 for n odd and ilu=1
387 DO I = K, N - 1
388 WORK( I ) = ZERO
389 END DO
390 DO J = K - 1, 0, -1
391 S = ZERO
392 DO I = 0, J - 2
393 AA = ABS( A( I+J*LDA ) )
394* -> A(j+k,i+k)
395 S = S + AA
396 WORK( I+K ) = WORK( I+K ) + AA
397 END DO
398.GT. IF( J0 ) THEN
399 AA = ABS( A( I+J*LDA ) )
400* -> A(j+k,j+k)
401 S = S + AA
402 WORK( I+K ) = WORK( I+K ) + S
403* i=j
404 I = I + 1
405 END IF
406 AA = ABS( A( I+J*LDA ) )
407* -> A(j,j)
408 WORK( J ) = AA
409 S = ZERO
410 DO L = J + 1, N - 1
411 I = I + 1
412 AA = ABS( A( I+J*LDA ) )
413* -> A(l,j)
414 S = S + AA
415 WORK( L ) = WORK( L ) + AA
416 END DO
417 WORK( J ) = WORK( J ) + S
418 END DO
419 VALUE = WORK( 0 )
420 DO I = 1, N-1
421 TEMP = WORK( I )
422.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
423 $ VALUE = TEMP
424 END DO
425 END IF
426 ELSE
427* n is even
428.EQ. IF( ILU0 ) THEN
429 DO I = 0, K - 1
430 WORK( I ) = ZERO
431 END DO
432 DO J = 0, K - 1
433 S = ZERO
434 DO I = 0, K + J - 1
435 AA = ABS( A( I+J*LDA ) )
436* -> A(i,j+k)
437 S = S + AA
438 WORK( I ) = WORK( I ) + AA
439 END DO
440 AA = ABS( A( I+J*LDA ) )
441* -> A(j+k,j+k)
442 WORK( J+K ) = S + AA
443 I = I + 1
444 AA = ABS( A( I+J*LDA ) )
445* -> A(j,j)
446 WORK( J ) = WORK( J ) + AA
447 S = ZERO
448 DO L = J + 1, K - 1
449 I = I + 1
450 AA = ABS( A( I+J*LDA ) )
451* -> A(l,j)
452 S = S + AA
453 WORK( L ) = WORK( L ) + AA
454 END DO
455 WORK( J ) = WORK( J ) + S
456 END DO
457 VALUE = WORK( 0 )
458 DO I = 1, N-1
459 TEMP = WORK( I )
460.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
461 $ VALUE = TEMP
462 END DO
463 ELSE
464* ilu = 1
465 DO I = K, N - 1
466 WORK( I ) = ZERO
467 END DO
468 DO J = K - 1, 0, -1
469 S = ZERO
470 DO I = 0, J - 1
471 AA = ABS( A( I+J*LDA ) )
472* -> A(j+k,i+k)
473 S = S + AA
474 WORK( I+K ) = WORK( I+K ) + AA
475 END DO
476 AA = ABS( A( I+J*LDA ) )
477* -> A(j+k,j+k)
478 S = S + AA
479 WORK( I+K ) = WORK( I+K ) + S
480* i=j
481 I = I + 1
482 AA = ABS( A( I+J*LDA ) )
483* -> A(j,j)
484 WORK( J ) = AA
485 S = ZERO
486 DO L = J + 1, N - 1
487 I = I + 1
488 AA = ABS( A( I+J*LDA ) )
489* -> A(l,j)
490 S = S + AA
491 WORK( L ) = WORK( L ) + AA
492 END DO
493 WORK( J ) = WORK( J ) + S
494 END DO
495 VALUE = WORK( 0 )
496 DO I = 1, N-1
497 TEMP = WORK( I )
498.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
499 $ VALUE = TEMP
500 END DO
501 END IF
502 END IF
503 ELSE
504* ifm=0
505 K = N / 2
506.EQ. IF( NOE1 ) THEN
507* n is odd
508.EQ. IF( ILU0 ) THEN
509 N1 = K
510* n/2
511 K = K + 1
512* k is the row size and lda
513 DO I = N1, N - 1
514 WORK( I ) = ZERO
515 END DO
516 DO J = 0, N1 - 1
517 S = ZERO
518 DO I = 0, K - 1
519 AA = ABS( A( I+J*LDA ) )
520* A(j,n1+i)
521 WORK( I+N1 ) = WORK( I+N1 ) + AA
522 S = S + AA
523 END DO
524 WORK( J ) = S
525 END DO
526* j=n1=k-1 is special
527 S = ABS( A( 0+J*LDA ) )
528* A(k-1,k-1)
529 DO I = 1, K - 1
530 AA = ABS( A( I+J*LDA ) )
531* A(k-1,i+n1)
532 WORK( I+N1 ) = WORK( I+N1 ) + AA
533 S = S + AA
534 END DO
535 WORK( J ) = WORK( J ) + S
536 DO J = K, N - 1
537 S = ZERO
538 DO I = 0, J - K - 1
539 AA = ABS( A( I+J*LDA ) )
540* A(i,j-k)
541 WORK( I ) = WORK( I ) + AA
542 S = S + AA
543 END DO
544* i=j-k
545 AA = ABS( A( I+J*LDA ) )
546* A(j-k,j-k)
547 S = S + AA
548 WORK( J-K ) = WORK( J-K ) + S
549 I = I + 1
550 S = ABS( A( I+J*LDA ) )
551* A(j,j)
552 DO L = J + 1, N - 1
553 I = I + 1
554 AA = ABS( A( I+J*LDA ) )
555* A(j,l)
556 WORK( L ) = WORK( L ) + AA
557 S = S + AA
558 END DO
559 WORK( J ) = WORK( J ) + S
560 END DO
561 VALUE = WORK( 0 )
562 DO I = 1, N-1
563 TEMP = WORK( I )
564.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
565 $ VALUE = TEMP
566 END DO
567 ELSE
568* ilu=1
569 K = K + 1
570* k=(n+1)/2 for n odd and ilu=1
571 DO I = K, N - 1
572 WORK( I ) = ZERO
573 END DO
574 DO J = 0, K - 2
575* process
576 S = ZERO
577 DO I = 0, J - 1
578 AA = ABS( A( I+J*LDA ) )
579* A(j,i)
580 WORK( I ) = WORK( I ) + AA
581 S = S + AA
582 END DO
583 AA = ABS( A( I+J*LDA ) )
584* i=j so process of A(j,j)
585 S = S + AA
586 WORK( J ) = S
587* is initialised here
588 I = I + 1
589* i=j process A(j+k,j+k)
590 AA = ABS( A( I+J*LDA ) )
591 S = AA
592 DO L = K + J + 1, N - 1
593 I = I + 1
594 AA = ABS( A( I+J*LDA ) )
595* A(l,k+j)
596 S = S + AA
597 WORK( L ) = WORK( L ) + AA
598 END DO
599 WORK( K+J ) = WORK( K+J ) + S
600 END DO
601* j=k-1 is special :process col A(k-1,0:k-1)
602 S = ZERO
603 DO I = 0, K - 2
604 AA = ABS( A( I+J*LDA ) )
605* A(k,i)
606 WORK( I ) = WORK( I ) + AA
607 S = S + AA
608 END DO
609* i=k-1
610 AA = ABS( A( I+J*LDA ) )
611* A(k-1,k-1)
612 S = S + AA
613 WORK( I ) = S
614* done with col j=k+1
615 DO J = K, N - 1
616* process col j of A = A(j,0:k-1)
617 S = ZERO
618 DO I = 0, K - 1
619 AA = ABS( A( I+J*LDA ) )
620* A(j,i)
621 WORK( I ) = WORK( I ) + AA
622 S = S + AA
623 END DO
624 WORK( J ) = WORK( J ) + S
625 END DO
626 VALUE = WORK( 0 )
627 DO I = 1, N-1
628 TEMP = WORK( I )
629.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
630 $ VALUE = TEMP
631 END DO
632 END IF
633 ELSE
634* n is even
635.EQ. IF( ILU0 ) THEN
636 DO I = K, N - 1
637 WORK( I ) = ZERO
638 END DO
639 DO J = 0, K - 1
640 S = ZERO
641 DO I = 0, K - 1
642 AA = ABS( A( I+J*LDA ) )
643* A(j,i+k)
644 WORK( I+K ) = WORK( I+K ) + AA
645 S = S + AA
646 END DO
647 WORK( J ) = S
648 END DO
649* j=k
650 AA = ABS( A( 0+J*LDA ) )
651* A(k,k)
652 S = AA
653 DO I = 1, K - 1
654 AA = ABS( A( I+J*LDA ) )
655* A(k,k+i)
656 WORK( I+K ) = WORK( I+K ) + AA
657 S = S + AA
658 END DO
659 WORK( J ) = WORK( J ) + S
660 DO J = K + 1, N - 1
661 S = ZERO
662 DO I = 0, J - 2 - K
663 AA = ABS( A( I+J*LDA ) )
664* A(i,j-k-1)
665 WORK( I ) = WORK( I ) + AA
666 S = S + AA
667 END DO
668* i=j-1-k
669 AA = ABS( A( I+J*LDA ) )
670* A(j-k-1,j-k-1)
671 S = S + AA
672 WORK( J-K-1 ) = WORK( J-K-1 ) + S
673 I = I + 1
674 AA = ABS( A( I+J*LDA ) )
675* A(j,j)
676 S = AA
677 DO L = J + 1, N - 1
678 I = I + 1
679 AA = ABS( A( I+J*LDA ) )
680* A(j,l)
681 WORK( L ) = WORK( L ) + AA
682 S = S + AA
683 END DO
684 WORK( J ) = WORK( J ) + S
685 END DO
686* j=n
687 S = ZERO
688 DO I = 0, K - 2
689 AA = ABS( A( I+J*LDA ) )
690* A(i,k-1)
691 WORK( I ) = WORK( I ) + AA
692 S = S + AA
693 END DO
694* i=k-1
695 AA = ABS( A( I+J*LDA ) )
696* A(k-1,k-1)
697 S = S + AA
698 WORK( I ) = WORK( I ) + S
699 VALUE = WORK ( 0 )
700 DO I = 1, N-1
701 TEMP = WORK( I )
702.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
703 $ VALUE = TEMP
704 END DO
705 ELSE
706* ilu=1
707 DO I = K, N - 1
708 WORK( I ) = ZERO
709 END DO
710* j=0 is special :process col A(k:n-1,k)
711 S = ABS( A( 0 ) )
712* A(k,k)
713 DO I = 1, K - 1
714 AA = ABS( A( I ) )
715* A(k+i,k)
716 WORK( I+K ) = WORK( I+K ) + AA
717 S = S + AA
718 END DO
719 WORK( K ) = WORK( K ) + S
720 DO J = 1, K - 1
721* process
722 S = ZERO
723 DO I = 0, J - 2
724 AA = ABS( A( I+J*LDA ) )
725* A(j-1,i)
726 WORK( I ) = WORK( I ) + AA
727 S = S + AA
728 END DO
729 AA = ABS( A( I+J*LDA ) )
730* i=j-1 so process of A(j-1,j-1)
731 S = S + AA
732 WORK( J-1 ) = S
733* is initialised here
734 I = I + 1
735* i=j process A(j+k,j+k)
736 AA = ABS( A( I+J*LDA ) )
737 S = AA
738 DO L = K + J + 1, N - 1
739 I = I + 1
740 AA = ABS( A( I+J*LDA ) )
741* A(l,k+j)
742 S = S + AA
743 WORK( L ) = WORK( L ) + AA
744 END DO
745 WORK( K+J ) = WORK( K+J ) + S
746 END DO
747* j=k is special :process col A(k,0:k-1)
748 S = ZERO
749 DO I = 0, K - 2
750 AA = ABS( A( I+J*LDA ) )
751* A(k,i)
752 WORK( I ) = WORK( I ) + AA
753 S = S + AA
754 END DO
755* i=k-1
756 AA = ABS( A( I+J*LDA ) )
757* A(k-1,k-1)
758 S = S + AA
759 WORK( I ) = S
760* done with col j=k+1
761 DO J = K + 1, N
762* process col j-1 of A = A(j-1,0:k-1)
763 S = ZERO
764 DO I = 0, K - 1
765 AA = ABS( A( I+J*LDA ) )
766* A(j-1,i)
767 WORK( I ) = WORK( I ) + AA
768 S = S + AA
769 END DO
770 WORK( J-1 ) = WORK( J-1 ) + S
771 END DO
772 VALUE = WORK( 0 )
773 DO I = 1, N-1
774 TEMP = WORK( I )
775.LT..OR. IF( VALUE TEMP SISNAN( TEMP ) )
776 $ VALUE = TEMP
777 END DO
778 END IF
779 END IF
780 END IF
781 ELSE IF( ( LSAME( NORM, 'f.OR.' ) ) ( LSAME( NORM, 'e' ) ) ) THEN
782*
783* Find normF(A).
784*
785 K = ( N+1 ) / 2
786 SCALE = ZERO
787 S = ONE
788.EQ. IF( NOE1 ) THEN
789* n is odd
790.EQ. IF( IFM1 ) THEN
791* A is normal
792.EQ. IF( ILU0 ) THEN
793* A is upper
794 DO J = 0, K - 3
795 CALL SLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S )
796* L at A(k,0)
797 END DO
798 DO J = 0, K - 1
799 CALL SLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S )
800* trap U at A(0,0)
801 END DO
802 S = S + S
803* double s for the off diagonal elements
804 CALL SLASSQ( K-1, A( K ), LDA+1, SCALE, S )
805* tri L at A(k,0)
806 CALL SLASSQ( K, A( K-1 ), LDA+1, SCALE, S )
807* tri U at A(k-1,0)
808 ELSE
809* ilu=1 & A is lower
810 DO J = 0, K - 1
811 CALL SLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S )
812* trap L at A(0,0)
813 END DO
814 DO J = 0, K - 2
815 CALL SLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S )
816* U at A(0,1)
817 END DO
818 S = S + S
819* double s for the off diagonal elements
820 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
821* tri L at A(0,0)
822 CALL SLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S )
823* tri U at A(0,1)
824 END IF
825 ELSE
826* A is xpose
827.EQ. IF( ILU0 ) THEN
828* A**T is upper
829 DO J = 1, K - 2
830 CALL SLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )
831* U at A(0,k)
832 END DO
833 DO J = 0, K - 2
834 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
835* k by k-1 rect. at A(0,0)
836 END DO
837 DO J = 0, K - 2
838 CALL SLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,
839 $ SCALE, S )
840* L at A(0,k-1)
841 END DO
842 S = S + S
843* double s for the off diagonal elements
844 CALL SLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S )
845* tri U at A(0,k)
846 CALL SLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )
847* tri L at A(0,k-1)
848 ELSE
849* A**T is lower
850 DO J = 1, K - 1
851 CALL SLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
852* U at A(0,0)
853 END DO
854 DO J = K, N - 1
855 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
856* k by k-1 rect. at A(0,k)
857 END DO
858 DO J = 0, K - 3
859 CALL SLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S )
860* L at A(1,0)
861 END DO
862 S = S + S
863* double s for the off diagonal elements
864 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
865* tri U at A(0,0)
866 CALL SLASSQ( K-1, A( 1 ), LDA+1, SCALE, S )
867* tri L at A(1,0)
868 END IF
869 END IF
870 ELSE
871* n is even
872.EQ. IF( IFM1 ) THEN
873* A is normal
874.EQ. IF( ILU0 ) THEN
875* A is upper
876 DO J = 0, K - 2
877 CALL SLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S )
878* L at A(k+1,0)
879 END DO
880 DO J = 0, K - 1
881 CALL SLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S )
882* trap U at A(0,0)
883 END DO
884 S = S + S
885* double s for the off diagonal elements
886 CALL SLASSQ( K, A( K+1 ), LDA+1, SCALE, S )
887* tri L at A(k+1,0)
888 CALL SLASSQ( K, A( K ), LDA+1, SCALE, S )
889* tri U at A(k,0)
890 ELSE
891* ilu=1 & A is lower
892 DO J = 0, K - 1
893 CALL SLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S )
894* trap L at A(1,0)
895 END DO
896 DO J = 1, K - 1
897 CALL SLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
898* U at A(0,0)
899 END DO
900 S = S + S
901* double s for the off diagonal elements
902 CALL SLASSQ( K, A( 1 ), LDA+1, SCALE, S )
903* tri L at A(1,0)
904 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
905* tri U at A(0,0)
906 END IF
907 ELSE
908* A is xpose
909.EQ. IF( ILU0 ) THEN
910* A**T is upper
911 DO J = 1, K - 1
912 CALL SLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )
913* U at A(0,k+1)
914 END DO
915 DO J = 0, K - 1
916 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
917* k by k rect. at A(0,0)
918 END DO
919 DO J = 0, K - 2
920 CALL SLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,
921 $ S )
922* L at A(0,k)
923 END DO
924 S = S + S
925* double s for the off diagonal elements
926 CALL SLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S )
927* tri U at A(0,k+1)
928 CALL SLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )
929* tri L at A(0,k)
930 ELSE
931* A**T is lower
932 DO J = 1, K - 1
933 CALL SLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )
934* U at A(0,1)
935 END DO
936 DO J = K + 1, N
937 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
938* k by k rect. at A(0,k+1)
939 END DO
940 DO J = 0, K - 2
941 CALL SLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S )
942* L at A(0,0)
943 END DO
944 S = S + S
945* double s for the off diagonal elements
946 CALL SLASSQ( K, A( LDA ), LDA+1, SCALE, S )
947* tri L at A(0,1)
948 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
949* tri U at A(0,0)
950 END IF
951 END IF
952 END IF
953 VALUE = SCALE*SQRT( S )
954 END IF
955*
956 SLANSF = VALUE
957 RETURN
958*
959* End of SLANSF
960*
961 END
norm
norm(diag(diag(diag(inv(mat))) -id.SOL), 2) % destroy mumps instance id.JOB
slassq
subroutine slassq(n, x, incx, scl, sumsq)
SLASSQ updates a sum of squares represented in scaled form.
Definition slassq.f90:137
sisnan
logical function sisnan(sin)
SISNAN tests input for NaN.
Definition sisnan.f:59
lsame
logical function lsame(ca, cb)
LSAME
Definition lsame.f:53
slansf
real function slansf(norm, transr, uplo, n, a, work)
SLANSF
Definition slansf.f:209