1 DOUBLE PRECISION FUNCTION DLANSP( NORM, UPLO, N, AP, WORK )
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER NORM, UPLO
10 INTEGER N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLANSP returns the value of the one norm, or the Frobenius norm, or
20 * the infinity norm, or the element of largest absolute value of a
21 * real symmetric matrix A, supplied in packed form.
22 *
23 * Description
24 * ===========
25 *
26 * DLANSP returns the value
27 *
28 * DLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
29 * (
30 * ( norm1(A), NORM = '1', 'O' or 'o'
31 * (
32 * ( normI(A), NORM = 'I' or 'i'
33 * (
34 * ( normF(A), NORM = 'F', 'f', 'E' or 'e'
35 *
36 * where norm1 denotes the one norm of a matrix (maximum column sum),
37 * normI denotes the infinity norm of a matrix (maximum row sum) and
38 * normF denotes the Frobenius norm of a matrix (square root of sum of
39 * squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.
40 *
41 * Arguments
42 * =========
43 *
44 * NORM (input) CHARACTER*1
45 * Specifies the value to be returned in DLANSP as described
46 * above.
47 *
48 * UPLO (input) CHARACTER*1
49 * Specifies whether the upper or lower triangular part of the
50 * symmetric matrix A is supplied.
51 * = 'U': Upper triangular part of A is supplied
52 * = 'L': Lower triangular part of A is supplied
53 *
54 * N (input) INTEGER
55 * The order of the matrix A. N >= 0. When N = 0, DLANSP is
56 * set to zero.
57 *
58 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
59 * The upper or lower triangle of the symmetric matrix A, packed
60 * columnwise in a linear array. The j-th column of A is stored
61 * in the array AP as follows:
62 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
63 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
64 *
65 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
66 * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
67 * WORK is not referenced.
68 *
69 * =====================================================================
70 *
71 * .. Parameters ..
72 DOUBLE PRECISION ONE, ZERO
73 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
74 * ..
75 * .. Local Scalars ..
76 INTEGER I, J, K
77 DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL DLASSQ
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 EXTERNAL LSAME
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC ABS, MAX, SQRT
88 * ..
89 * .. Executable Statements ..
90 *
91 IF( N.EQ.0 ) THEN
92 VALUE = ZERO
93 ELSE IF( LSAME( NORM, 'M' ) ) THEN
94 *
95 * Find max(abs(A(i,j))).
96 *
97 VALUE = ZERO
98 IF( LSAME( UPLO, 'U' ) ) THEN
99 K = 1
100 DO 20 J = 1, N
101 DO 10 I = K, K + J - 1
102 VALUE = MAX( VALUE, ABS( AP( I ) ) )
103 10 CONTINUE
104 K = K + J
105 20 CONTINUE
106 ELSE
107 K = 1
108 DO 40 J = 1, N
109 DO 30 I = K, K + N - J
110 VALUE = MAX( VALUE, ABS( AP( I ) ) )
111 30 CONTINUE
112 K = K + N - J + 1
113 40 CONTINUE
114 END IF
115 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
116 $ ( NORM.EQ.'1' ) ) THEN
117 *
118 * Find normI(A) ( = norm1(A), since A is symmetric).
119 *
120 VALUE = ZERO
121 K = 1
122 IF( LSAME( UPLO, 'U' ) ) THEN
123 DO 60 J = 1, N
124 SUM = ZERO
125 DO 50 I = 1, J - 1
126 ABSA = ABS( AP( K ) )
127 SUM = SUM + ABSA
128 WORK( I ) = WORK( I ) + ABSA
129 K = K + 1
130 50 CONTINUE
131 WORK( J ) = SUM + ABS( AP( K ) )
132 K = K + 1
133 60 CONTINUE
134 DO 70 I = 1, N
135 VALUE = MAX( VALUE, WORK( I ) )
136 70 CONTINUE
137 ELSE
138 DO 80 I = 1, N
139 WORK( I ) = ZERO
140 80 CONTINUE
141 DO 100 J = 1, N
142 SUM = WORK( J ) + ABS( AP( K ) )
143 K = K + 1
144 DO 90 I = J + 1, N
145 ABSA = ABS( AP( K ) )
146 SUM = SUM + ABSA
147 WORK( I ) = WORK( I ) + ABSA
148 K = K + 1
149 90 CONTINUE
150 VALUE = MAX( VALUE, SUM )
151 100 CONTINUE
152 END IF
153 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
154 *
155 * Find normF(A).
156 *
157 SCALE = ZERO
158 SUM = ONE
159 K = 2
160 IF( LSAME( UPLO, 'U' ) ) THEN
161 DO 110 J = 2, N
162 CALL DLASSQ( J-1, AP( K ), 1, SCALE, SUM )
163 K = K + J
164 110 CONTINUE
165 ELSE
166 DO 120 J = 1, N - 1
167 CALL DLASSQ( N-J, AP( K ), 1, SCALE, SUM )
168 K = K + N - J + 1
169 120 CONTINUE
170 END IF
171 SUM = 2*SUM
172 K = 1
173 DO 130 I = 1, N
174 IF( AP( K ).NE.ZERO ) THEN
175 ABSA = ABS( AP( K ) )
176 IF( SCALE.LT.ABSA ) THEN
177 SUM = ONE + SUM*( SCALE / ABSA )**2
178 SCALE = ABSA
179 ELSE
180 SUM = SUM + ( ABSA / SCALE )**2
181 END IF
182 END IF
183 IF( LSAME( UPLO, 'U' ) ) THEN
184 K = K + I + 1
185 ELSE
186 K = K + N - I + 1
187 END IF
188 130 CONTINUE
189 VALUE = SCALE*SQRT( SUM )
190 END IF
191 *
192 DLANSP = VALUE
193 RETURN
194 *
195 * End of DLANSP
196 *
197 END
2 *
3 * -- LAPACK auxiliary routine (version 3.2) --
4 * -- LAPACK is a software package provided by Univ. of Tennessee, --
5 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
6 * November 2006
7 *
8 * .. Scalar Arguments ..
9 CHARACTER NORM, UPLO
10 INTEGER N
11 * ..
12 * .. Array Arguments ..
13 DOUBLE PRECISION AP( * ), WORK( * )
14 * ..
15 *
16 * Purpose
17 * =======
18 *
19 * DLANSP returns the value of the one norm, or the Frobenius norm, or
20 * the infinity norm, or the element of largest absolute value of a
21 * real symmetric matrix A, supplied in packed form.
22 *
23 * Description
24 * ===========
25 *
26 * DLANSP returns the value
27 *
28 * DLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
29 * (
30 * ( norm1(A), NORM = '1', 'O' or 'o'
31 * (
32 * ( normI(A), NORM = 'I' or 'i'
33 * (
34 * ( normF(A), NORM = 'F', 'f', 'E' or 'e'
35 *
36 * where norm1 denotes the one norm of a matrix (maximum column sum),
37 * normI denotes the infinity norm of a matrix (maximum row sum) and
38 * normF denotes the Frobenius norm of a matrix (square root of sum of
39 * squares). Note that max(abs(A(i,j))) is not a consistent matrix norm.
40 *
41 * Arguments
42 * =========
43 *
44 * NORM (input) CHARACTER*1
45 * Specifies the value to be returned in DLANSP as described
46 * above.
47 *
48 * UPLO (input) CHARACTER*1
49 * Specifies whether the upper or lower triangular part of the
50 * symmetric matrix A is supplied.
51 * = 'U': Upper triangular part of A is supplied
52 * = 'L': Lower triangular part of A is supplied
53 *
54 * N (input) INTEGER
55 * The order of the matrix A. N >= 0. When N = 0, DLANSP is
56 * set to zero.
57 *
58 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
59 * The upper or lower triangle of the symmetric matrix A, packed
60 * columnwise in a linear array. The j-th column of A is stored
61 * in the array AP as follows:
62 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
63 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
64 *
65 * WORK (workspace) DOUBLE PRECISION array, dimension (MAX(1,LWORK)),
66 * where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
67 * WORK is not referenced.
68 *
69 * =====================================================================
70 *
71 * .. Parameters ..
72 DOUBLE PRECISION ONE, ZERO
73 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
74 * ..
75 * .. Local Scalars ..
76 INTEGER I, J, K
77 DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL DLASSQ
81 * ..
82 * .. External Functions ..
83 LOGICAL LSAME
84 EXTERNAL LSAME
85 * ..
86 * .. Intrinsic Functions ..
87 INTRINSIC ABS, MAX, SQRT
88 * ..
89 * .. Executable Statements ..
90 *
91 IF( N.EQ.0 ) THEN
92 VALUE = ZERO
93 ELSE IF( LSAME( NORM, 'M' ) ) THEN
94 *
95 * Find max(abs(A(i,j))).
96 *
97 VALUE = ZERO
98 IF( LSAME( UPLO, 'U' ) ) THEN
99 K = 1
100 DO 20 J = 1, N
101 DO 10 I = K, K + J - 1
102 VALUE = MAX( VALUE, ABS( AP( I ) ) )
103 10 CONTINUE
104 K = K + J
105 20 CONTINUE
106 ELSE
107 K = 1
108 DO 40 J = 1, N
109 DO 30 I = K, K + N - J
110 VALUE = MAX( VALUE, ABS( AP( I ) ) )
111 30 CONTINUE
112 K = K + N - J + 1
113 40 CONTINUE
114 END IF
115 ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
116 $ ( NORM.EQ.'1' ) ) THEN
117 *
118 * Find normI(A) ( = norm1(A), since A is symmetric).
119 *
120 VALUE = ZERO
121 K = 1
122 IF( LSAME( UPLO, 'U' ) ) THEN
123 DO 60 J = 1, N
124 SUM = ZERO
125 DO 50 I = 1, J - 1
126 ABSA = ABS( AP( K ) )
127 SUM = SUM + ABSA
128 WORK( I ) = WORK( I ) + ABSA
129 K = K + 1
130 50 CONTINUE
131 WORK( J ) = SUM + ABS( AP( K ) )
132 K = K + 1
133 60 CONTINUE
134 DO 70 I = 1, N
135 VALUE = MAX( VALUE, WORK( I ) )
136 70 CONTINUE
137 ELSE
138 DO 80 I = 1, N
139 WORK( I ) = ZERO
140 80 CONTINUE
141 DO 100 J = 1, N
142 SUM = WORK( J ) + ABS( AP( K ) )
143 K = K + 1
144 DO 90 I = J + 1, N
145 ABSA = ABS( AP( K ) )
146 SUM = SUM + ABSA
147 WORK( I ) = WORK( I ) + ABSA
148 K = K + 1
149 90 CONTINUE
150 VALUE = MAX( VALUE, SUM )
151 100 CONTINUE
152 END IF
153 ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
154 *
155 * Find normF(A).
156 *
157 SCALE = ZERO
158 SUM = ONE
159 K = 2
160 IF( LSAME( UPLO, 'U' ) ) THEN
161 DO 110 J = 2, N
162 CALL DLASSQ( J-1, AP( K ), 1, SCALE, SUM )
163 K = K + J
164 110 CONTINUE
165 ELSE
166 DO 120 J = 1, N - 1
167 CALL DLASSQ( N-J, AP( K ), 1, SCALE, SUM )
168 K = K + N - J + 1
169 120 CONTINUE
170 END IF
171 SUM = 2*SUM
172 K = 1
173 DO 130 I = 1, N
174 IF( AP( K ).NE.ZERO ) THEN
175 ABSA = ABS( AP( K ) )
176 IF( SCALE.LT.ABSA ) THEN
177 SUM = ONE + SUM*( SCALE / ABSA )**2
178 SCALE = ABSA
179 ELSE
180 SUM = SUM + ( ABSA / SCALE )**2
181 END IF
182 END IF
183 IF( LSAME( UPLO, 'U' ) ) THEN
184 K = K + I + 1
185 ELSE
186 K = K + N - I + 1
187 END IF
188 130 CONTINUE
189 VALUE = SCALE*SQRT( SUM )
190 END IF
191 *
192 DLANSP = VALUE
193 RETURN
194 *
195 * End of DLANSP
196 *
197 END