1
DOUBLE PRECISION FUNCTION DLANSP( NORM, UPLO, N, AP, WORK )
3
* -- LAPACK auxiliary routine (version 2.0) --
4
* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
5
* Courant Institute, Argonne National Lab, and Rice University
8
* .. Scalar Arguments ..
12
* .. Array Arguments ..
13
DOUBLE PRECISION AP( * ), WORK( * )
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.
26
* DLANSP returns the value
28
* DLANSP = ( max(abs(A(i,j))), NORM = 'M' or 'm'
30
* ( norm1(A), NORM = '1', 'O' or 'o'
32
* ( normI(A), NORM = 'I' or 'i'
34
* ( normF(A), NORM = 'F', 'f', 'E' or 'e'
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 matrix norm.
44
* NORM (input) CHARACTER*1
45
* Specifies the value to be returned in DLANSP as described
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
55
* The order of the matrix A. N >= 0. When N = 0, DLANSP is
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.
65
* WORK (workspace) DOUBLE PRECISION array, dimension (LWORK),
66
* where LWORK >= N when NORM = 'I' or '1' or 'O'; otherwise,
67
* WORK is not referenced.
69
* =====================================================================
72
DOUBLE PRECISION ONE, ZERO
73
PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
77
DOUBLE PRECISION ABSA, SCALE, SUM, VALUE
79
* .. External Subroutines ..
82
* .. External Functions ..
86
* .. Intrinsic Functions ..
87
INTRINSIC ABS, MAX, SQRT
89
* .. Executable Statements ..
93
ELSE IF( LSAME( NORM, 'M' ) ) THEN
95
* Find max(abs(A(i,j))).
98
IF( LSAME( UPLO, 'U' ) ) THEN
101
DO 10 I = K, K + J - 1
102
VALUE = MAX( VALUE, ABS( AP( I ) ) )
109
DO 30 I = K, K + N - J
110
VALUE = MAX( VALUE, ABS( AP( I ) ) )
115
ELSE IF( ( LSAME( NORM, 'I' ) ) .OR. ( LSAME( NORM, 'O' ) ) .OR.
116
$ ( NORM.EQ.'1' ) ) THEN
118
* Find normI(A) ( = norm1(A), since A is symmetric).
122
IF( LSAME( UPLO, 'U' ) ) THEN
126
ABSA = ABS( AP( K ) )
128
WORK( I ) = WORK( I ) + ABSA
131
WORK( J ) = SUM + ABS( AP( K ) )
135
VALUE = MAX( VALUE, WORK( I ) )
142
SUM = WORK( J ) + ABS( AP( K ) )
145
ABSA = ABS( AP( K ) )
147
WORK( I ) = WORK( I ) + ABSA
150
VALUE = MAX( VALUE, SUM )
153
ELSE IF( ( LSAME( NORM, 'F' ) ) .OR. ( LSAME( NORM, 'E' ) ) ) THEN
160
IF( LSAME( UPLO, 'U' ) ) THEN
162
CALL DLASSQ( J-1, AP( K ), 1, SCALE, SUM )
167
CALL DLASSQ( N-J, AP( K ), 1, SCALE, SUM )
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
180
SUM = SUM + ( ABSA / SCALE )**2
183
IF( LSAME( UPLO, 'U' ) ) THEN
189
VALUE = SCALE*SQRT( SUM )