3 * =========== DOCUMENTATION ===========
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
9 *> Download SLANSF + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slansf.f">
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slansf.f">
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slansf.f">
21 * REAL FUNCTION SLANSF( NORM, TRANSR, UPLO, N, A, WORK )
23 * .. Scalar Arguments ..
24 * CHARACTER NORM, TRANSR, UPLO
27 * .. Array Arguments ..
28 * REAL A( 0: * ), WORK( 0: * )
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.
45 *> SLANSF = ( max(abs(A(i,j))), NORM = 'M' or 'm'
47 *> ( norm1(A), NORM = '1', 'O' or 'o'
49 *> ( normI(A), NORM = 'I' or 'i'
51 *> ( normF(A), NORM = 'F', 'f', 'E' or 'e'
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.
64 *> NORM is CHARACTER*1
65 *> Specifies the value to be returned in SLANSF as described
71 *> TRANSR is CHARACTER*1
72 *> Specifies whether the RFP format of A is normal or
74 *> = 'N': RFP format is Normal;
75 *> = 'T': RFP format is Transpose.
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.
90 *> The order of the matrix A. N >= 0. When N = 0, SLANSF is
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.
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.
113 *> \author Univ. of Tennessee
114 *> \author Univ. of California Berkeley
115 *> \author Univ. of Colorado Denver
118 *> \date September 2012
120 *> \ingroup realOTHERcomputational
122 *> \par Further Details:
123 * =====================
127 *> We first consider Rectangular Full Packed (RFP) Format when N is
128 *> even. We give an example where N = 6.
130 *> AP is Upper AP is Lower
132 *> 00 01 02 03 04 05 00
133 *> 11 12 13 14 15 10 11
134 *> 22 23 24 25 20 21 22
135 *> 33 34 35 30 31 32 33
136 *> 44 45 40 41 42 43 44
137 *> 55 50 51 52 53 54 55
140 *> Let TRANSR = 'N'. RFP holds AP as follows:
141 *> For UPLO = 'U' the upper trapezoid A(0:5,0:2) consists of the last
142 *> three columns of AP upper. The lower triangle A(4:6,0:2) consists of
143 *> the transpose of the first three columns of AP upper.
144 *> For UPLO = 'L' the lower trapezoid A(1:6,0:2) consists of the first
145 *> three columns of AP lower. The upper triangle A(0:2,0:2) consists of
146 *> the transpose of the last three columns of AP lower.
147 *> This covers the case N even and TRANSR = 'N'.
159 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
160 *> transpose of RFP A above. One therefore gets:
165 *> 03 13 23 33 00 01 02 33 00 10 20 30 40 50
166 *> 04 14 24 34 44 11 12 43 44 11 21 31 41 51
167 *> 05 15 25 35 45 55 22 53 54 55 22 32 42 52
170 *> We then consider Rectangular Full Packed (RFP) Format when N is
171 *> odd. We give an example where N = 5.
173 *> AP is Upper AP is Lower
182 *> Let TRANSR = 'N'. RFP holds AP as follows:
183 *> For UPLO = 'U' the upper trapezoid A(0:4,0:2) consists of the last
184 *> three columns of AP upper. The lower triangle A(3:4,0:1) consists of
185 *> the transpose of the first two columns of AP upper.
186 *> For UPLO = 'L' the lower trapezoid A(0:4,0:2) consists of the first
187 *> three columns of AP lower. The upper triangle A(0:1,1:2) consists of
188 *> the transpose of the last two columns of AP lower.
189 *> This covers the case N odd and TRANSR = 'N'.
199 *> Now let TRANSR = 'T'. RFP A in both UPLO cases is just the
200 *> transpose of RFP A above. One therefore gets:
204 *> 02 12 22 00 01 00 10 20 30 40 50
205 *> 03 13 23 33 11 33 11 21 31 41 51
206 *> 04 14 24 34 44 43 44 22 32 42 52
209 * =====================================================================
210 REAL FUNCTION SLANSF( NORM, TRANSR, UPLO, N, A, WORK )
212 * -- LAPACK computational routine (version 3.4.2) --
213 * -- LAPACK is a software package provided by Univ. of Tennessee, --
214 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
217 * .. Scalar Arguments ..
218 CHARACTER NORM, TRANSR, UPLO
221 * .. Array Arguments ..
222 REAL A( 0: * ), WORK( 0: * )
225 * =====================================================================
230 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 )
232 * .. Local Scalars ..
233 INTEGER I, J, IFM, ILU, NOE, N1, K, L, LDA
234 REAL SCALE, S, VALUE, AA, TEMP
236 * .. External Functions ..
237 LOGICAL LSAME, SISNAN
238 EXTERNAL LSAME, SISNAN
240 * .. External Subroutines ..
243 * .. Intrinsic Functions ..
246 * .. Executable Statements ..
251 ELSE IF( N.EQ.1 ) THEN
256 * set noe = 1 if n is odd. if n is even set noe=0
259 IF( MOD( N, 2 ).EQ.0 )
262 * set ifm = 0 when form='T or 't' and 1 otherwise
265 IF( LSAME( TRANSR, 'T' ) )
268 * set ilu = 0 when uplo='U or 'u' and 1 otherwise
271 IF( LSAME( UPLO, 'U' ) )
274 * set lda = (n+1)/2 when ifm = 0
275 * set lda = n when ifm = 1 and noe = 1
276 * set lda = n+1 when ifm = 1 and noe = 0
290 IF( LSAME( NORM, 'M' ) ) THEN
292 * Find max(abs(A(i,j))).
302 TEMP = ABS( A( I+J*LDA ) )
303 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
308 * xpose case; A is k by n
311 TEMP = ABS( A( I+J*LDA ) )
312 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
323 TEMP = ABS( A( I+J*LDA ) )
324 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
329 * xpose case; A is k by n+1
332 TEMP = ABS( A( I+J*LDA ) )
333 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
339 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
340 $ ( NORM.EQ.'1' ) ) THEN
342 * Find normI(A) ( = norm1(A), since A is symmetric).
355 AA = ABS( A( I+J*LDA ) )
358 WORK( I ) = WORK( I ) + AA
360 AA = ABS( A( I+J*LDA ) )
366 AA = ABS( A( I+J*LDA ) )
368 WORK( J ) = WORK( J ) + AA
372 AA = ABS( A( I+J*LDA ) )
375 WORK( L ) = WORK( L ) + AA
377 WORK( J ) = WORK( J ) + S
383 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
389 * k=(n+1)/2 for n odd and ilu=1
396 AA = ABS( A( I+J*LDA ) )
399 WORK( I+K ) = WORK( I+K ) + AA
402 AA = ABS( A( I+J*LDA ) )
405 WORK( I+K ) = WORK( I+K ) + S
409 AA = ABS( A( I+J*LDA ) )
415 AA = ABS( A( I+J*LDA ) )
418 WORK( L ) = WORK( L ) + AA
420 WORK( J ) = WORK( J ) + S
425 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
438 AA = ABS( A( I+J*LDA ) )
441 WORK( I ) = WORK( I ) + AA
443 AA = ABS( A( I+J*LDA ) )
447 AA = ABS( A( I+J*LDA ) )
449 WORK( J ) = WORK( J ) + AA
453 AA = ABS( A( I+J*LDA ) )
456 WORK( L ) = WORK( L ) + AA
458 WORK( J ) = WORK( J ) + S
463 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
474 AA = ABS( A( I+J*LDA ) )
477 WORK( I+K ) = WORK( I+K ) + AA
479 AA = ABS( A( I+J*LDA ) )
482 WORK( I+K ) = WORK( I+K ) + S
485 AA = ABS( A( I+J*LDA ) )
491 AA = ABS( A( I+J*LDA ) )
494 WORK( L ) = WORK( L ) + AA
496 WORK( J ) = WORK( J ) + S
501 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
515 * k is the row size and lda
522 AA = ABS( A( I+J*LDA ) )
524 WORK( I+N1 ) = WORK( I+N1 ) + AA
529 * j=n1=k-1 is special
530 S = ABS( A( 0+J*LDA ) )
533 AA = ABS( A( I+J*LDA ) )
535 WORK( I+N1 ) = WORK( I+N1 ) + AA
538 WORK( J ) = WORK( J ) + S
542 AA = ABS( A( I+J*LDA ) )
544 WORK( I ) = WORK( I ) + AA
548 AA = ABS( A( I+J*LDA ) )
551 WORK( J-K ) = WORK( J-K ) + S
553 S = ABS( A( I+J*LDA ) )
557 AA = ABS( A( I+J*LDA ) )
559 WORK( L ) = WORK( L ) + AA
562 WORK( J ) = WORK( J ) + S
567 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
573 * k=(n+1)/2 for n odd and ilu=1
581 AA = ABS( A( I+J*LDA ) )
583 WORK( I ) = WORK( I ) + AA
586 AA = ABS( A( I+J*LDA ) )
587 * i=j so process of A(j,j)
590 * is initialised here
592 * i=j process A(j+k,j+k)
593 AA = ABS( A( I+J*LDA ) )
595 DO L = K + J + 1, N - 1
597 AA = ABS( A( I+J*LDA ) )
600 WORK( L ) = WORK( L ) + AA
602 WORK( K+J ) = WORK( K+J ) + S
604 * j=k-1 is special :process col A(k-1,0:k-1)
607 AA = ABS( A( I+J*LDA ) )
609 WORK( I ) = WORK( I ) + AA
613 AA = ABS( A( I+J*LDA ) )
617 * done with col j=k+1
619 * process col j of A = A(j,0:k-1)
622 AA = ABS( A( I+J*LDA ) )
624 WORK( I ) = WORK( I ) + AA
627 WORK( J ) = WORK( J ) + S
632 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
645 AA = ABS( A( I+J*LDA ) )
647 WORK( I+K ) = WORK( I+K ) + AA
653 AA = ABS( A( 0+J*LDA ) )
657 AA = ABS( A( I+J*LDA ) )
659 WORK( I+K ) = WORK( I+K ) + AA
662 WORK( J ) = WORK( J ) + S
666 AA = ABS( A( I+J*LDA ) )
668 WORK( I ) = WORK( I ) + AA
672 AA = ABS( A( I+J*LDA ) )
675 WORK( J-K-1 ) = WORK( J-K-1 ) + S
677 AA = ABS( A( I+J*LDA ) )
682 AA = ABS( A( I+J*LDA ) )
684 WORK( L ) = WORK( L ) + AA
687 WORK( J ) = WORK( J ) + S
692 AA = ABS( A( I+J*LDA ) )
694 WORK( I ) = WORK( I ) + AA
698 AA = ABS( A( I+J*LDA ) )
701 WORK( I ) = WORK( I ) + S
705 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
713 * j=0 is special :process col A(k:n-1,k)
719 WORK( I+K ) = WORK( I+K ) + AA
722 WORK( K ) = WORK( K ) + S
727 AA = ABS( A( I+J*LDA ) )
729 WORK( I ) = WORK( I ) + AA
732 AA = ABS( A( I+J*LDA ) )
733 * i=j-1 so process of A(j-1,j-1)
736 * is initialised here
738 * i=j process A(j+k,j+k)
739 AA = ABS( A( I+J*LDA ) )
741 DO L = K + J + 1, N - 1
743 AA = ABS( A( I+J*LDA ) )
746 WORK( L ) = WORK( L ) + AA
748 WORK( K+J ) = WORK( K+J ) + S
750 * j=k is special :process col A(k,0:k-1)
753 AA = ABS( A( I+J*LDA ) )
755 WORK( I ) = WORK( I ) + AA
759 AA = ABS( A( I+J*LDA ) )
763 * done with col j=k+1
765 * process col j-1 of A = A(j-1,0:k-1)
768 AA = ABS( A( I+J*LDA ) )
770 WORK( I ) = WORK( I ) + AA
773 WORK( J-1 ) = WORK( J-1 ) + S
778 IF( VALUE .LT. TEMP .OR. SISNAN( TEMP ) )
784 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
798 CALL SLASSQ( K-J-2, A( K+J+1+J*LDA ), 1, SCALE, S )
802 CALL SLASSQ( K+J-1, A( 0+J*LDA ), 1, SCALE, S )
806 * double s for the off diagonal elements
807 CALL SLASSQ( K-1, A( K ), LDA+1, SCALE, S )
809 CALL SLASSQ( K, A( K-1 ), LDA+1, SCALE, S )
814 CALL SLASSQ( N-J-1, A( J+1+J*LDA ), 1, SCALE, S )
818 CALL SLASSQ( J, A( 0+( 1+J )*LDA ), 1, SCALE, S )
822 * double s for the off diagonal elements
823 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
825 CALL SLASSQ( K-1, A( 0+LDA ), LDA+1, SCALE, S )
833 CALL SLASSQ( J, A( 0+( K+J )*LDA ), 1, SCALE, S )
837 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
838 * k by k-1 rect. at A(0,0)
841 CALL SLASSQ( K-J-1, A( J+1+( J+K-1 )*LDA ), 1,
846 * double s for the off diagonal elements
847 CALL SLASSQ( K-1, A( 0+K*LDA ), LDA+1, SCALE, S )
849 CALL SLASSQ( K, A( 0+( K-1 )*LDA ), LDA+1, SCALE, S )
854 CALL SLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
858 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
859 * k by k-1 rect. at A(0,k)
862 CALL SLASSQ( K-J-2, A( J+2+J*LDA ), 1, SCALE, S )
866 * double s for the off diagonal elements
867 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
869 CALL SLASSQ( K-1, A( 1 ), LDA+1, SCALE, S )
880 CALL SLASSQ( K-J-1, A( K+J+2+J*LDA ), 1, SCALE, S )
884 CALL SLASSQ( K+J, A( 0+J*LDA ), 1, SCALE, S )
888 * double s for the off diagonal elements
889 CALL SLASSQ( K, A( K+1 ), LDA+1, SCALE, S )
891 CALL SLASSQ( K, A( K ), LDA+1, SCALE, S )
896 CALL SLASSQ( N-J-1, A( J+2+J*LDA ), 1, SCALE, S )
900 CALL SLASSQ( J, A( 0+J*LDA ), 1, SCALE, S )
904 * double s for the off diagonal elements
905 CALL SLASSQ( K, A( 1 ), LDA+1, SCALE, S )
907 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
915 CALL SLASSQ( J, A( 0+( K+1+J )*LDA ), 1, SCALE, S )
919 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
920 * k by k rect. at A(0,0)
923 CALL SLASSQ( K-J-1, A( J+1+( J+K )*LDA ), 1, SCALE,
928 * double s for the off diagonal elements
929 CALL SLASSQ( K, A( 0+( K+1 )*LDA ), LDA+1, SCALE, S )
931 CALL SLASSQ( K, A( 0+K*LDA ), LDA+1, SCALE, S )
936 CALL SLASSQ( J, A( 0+( J+1 )*LDA ), 1, SCALE, S )
940 CALL SLASSQ( K, A( 0+J*LDA ), 1, SCALE, S )
941 * k by k rect. at A(0,k+1)
944 CALL SLASSQ( K-J-1, A( J+1+J*LDA ), 1, SCALE, S )
948 * double s for the off diagonal elements
949 CALL SLASSQ( K, A( LDA ), LDA+1, SCALE, S )
951 CALL SLASSQ( K, A( 0 ), LDA+1, SCALE, S )
956 VALUE = SCALE*SQRT( S )