1 SUBROUTINE DPPEQU( UPLO, N, AP, S, SCOND, AMAX, INFO )
2 *
3 * -- LAPACK 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 UPLO
10 INTEGER INFO, N
11 DOUBLE PRECISION AMAX, SCOND
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AP( * ), S( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DPPEQU computes row and column scalings intended to equilibrate a
21 * symmetric positive definite matrix A in packed storage and reduce
22 * its condition number (with respect to the two-norm). S contains the
23 * scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix
24 * B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal.
25 * This choice of S puts the condition number of B within a factor N of
26 * the smallest possible condition number over all possible diagonal
27 * scalings.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER*1
33 * = 'U': Upper triangle of A is stored;
34 * = 'L': Lower triangle of A is stored.
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
39 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
40 * The upper or lower triangle of the symmetric matrix A, packed
41 * columnwise in a linear array. The j-th column of A is stored
42 * in the array AP as follows:
43 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
44 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
45 *
46 * S (output) DOUBLE PRECISION array, dimension (N)
47 * If INFO = 0, S contains the scale factors for A.
48 *
49 * SCOND (output) DOUBLE PRECISION
50 * If INFO = 0, S contains the ratio of the smallest S(i) to
51 * the largest S(i). If SCOND >= 0.1 and AMAX is neither too
52 * large nor too small, it is not worth scaling by S.
53 *
54 * AMAX (output) DOUBLE PRECISION
55 * Absolute value of largest matrix element. If AMAX is very
56 * close to overflow or very close to underflow, the matrix
57 * should be scaled.
58 *
59 * INFO (output) INTEGER
60 * = 0: successful exit
61 * < 0: if INFO = -i, the i-th argument had an illegal value
62 * > 0: if INFO = i, the i-th diagonal element is nonpositive.
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ONE, ZERO
68 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
69 * ..
70 * .. Local Scalars ..
71 LOGICAL UPPER
72 INTEGER I, JJ
73 DOUBLE PRECISION SMIN
74 * ..
75 * .. External Functions ..
76 LOGICAL LSAME
77 EXTERNAL LSAME
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC MAX, MIN, SQRT
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 INFO = 0
90 UPPER = LSAME( UPLO, 'U' )
91 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92 INFO = -1
93 ELSE IF( N.LT.0 ) THEN
94 INFO = -2
95 END IF
96 IF( INFO.NE.0 ) THEN
97 CALL XERBLA( 'DPPEQU', -INFO )
98 RETURN
99 END IF
100 *
101 * Quick return if possible
102 *
103 IF( N.EQ.0 ) THEN
104 SCOND = ONE
105 AMAX = ZERO
106 RETURN
107 END IF
108 *
109 * Initialize SMIN and AMAX.
110 *
111 S( 1 ) = AP( 1 )
112 SMIN = S( 1 )
113 AMAX = S( 1 )
114 *
115 IF( UPPER ) THEN
116 *
117 * UPLO = 'U': Upper triangle of A is stored.
118 * Find the minimum and maximum diagonal elements.
119 *
120 JJ = 1
121 DO 10 I = 2, N
122 JJ = JJ + I
123 S( I ) = AP( JJ )
124 SMIN = MIN( SMIN, S( I ) )
125 AMAX = MAX( AMAX, S( I ) )
126 10 CONTINUE
127 *
128 ELSE
129 *
130 * UPLO = 'L': Lower triangle of A is stored.
131 * Find the minimum and maximum diagonal elements.
132 *
133 JJ = 1
134 DO 20 I = 2, N
135 JJ = JJ + N - I + 2
136 S( I ) = AP( JJ )
137 SMIN = MIN( SMIN, S( I ) )
138 AMAX = MAX( AMAX, S( I ) )
139 20 CONTINUE
140 END IF
141 *
142 IF( SMIN.LE.ZERO ) THEN
143 *
144 * Find the first non-positive diagonal element and return.
145 *
146 DO 30 I = 1, N
147 IF( S( I ).LE.ZERO ) THEN
148 INFO = I
149 RETURN
150 END IF
151 30 CONTINUE
152 ELSE
153 *
154 * Set the scale factors to the reciprocals
155 * of the diagonal elements.
156 *
157 DO 40 I = 1, N
158 S( I ) = ONE / SQRT( S( I ) )
159 40 CONTINUE
160 *
161 * Compute SCOND = min(S(I)) / max(S(I))
162 *
163 SCOND = SQRT( SMIN ) / SQRT( AMAX )
164 END IF
165 RETURN
166 *
167 * End of DPPEQU
168 *
169 END
2 *
3 * -- LAPACK 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 UPLO
10 INTEGER INFO, N
11 DOUBLE PRECISION AMAX, SCOND
12 * ..
13 * .. Array Arguments ..
14 DOUBLE PRECISION AP( * ), S( * )
15 * ..
16 *
17 * Purpose
18 * =======
19 *
20 * DPPEQU computes row and column scalings intended to equilibrate a
21 * symmetric positive definite matrix A in packed storage and reduce
22 * its condition number (with respect to the two-norm). S contains the
23 * scale factors, S(i)=1/sqrt(A(i,i)), chosen so that the scaled matrix
24 * B with elements B(i,j)=S(i)*A(i,j)*S(j) has ones on the diagonal.
25 * This choice of S puts the condition number of B within a factor N of
26 * the smallest possible condition number over all possible diagonal
27 * scalings.
28 *
29 * Arguments
30 * =========
31 *
32 * UPLO (input) CHARACTER*1
33 * = 'U': Upper triangle of A is stored;
34 * = 'L': Lower triangle of A is stored.
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
39 * AP (input) DOUBLE PRECISION array, dimension (N*(N+1)/2)
40 * The upper or lower triangle of the symmetric matrix A, packed
41 * columnwise in a linear array. The j-th column of A is stored
42 * in the array AP as follows:
43 * if UPLO = 'U', AP(i + (j-1)*j/2) = A(i,j) for 1<=i<=j;
44 * if UPLO = 'L', AP(i + (j-1)*(2n-j)/2) = A(i,j) for j<=i<=n.
45 *
46 * S (output) DOUBLE PRECISION array, dimension (N)
47 * If INFO = 0, S contains the scale factors for A.
48 *
49 * SCOND (output) DOUBLE PRECISION
50 * If INFO = 0, S contains the ratio of the smallest S(i) to
51 * the largest S(i). If SCOND >= 0.1 and AMAX is neither too
52 * large nor too small, it is not worth scaling by S.
53 *
54 * AMAX (output) DOUBLE PRECISION
55 * Absolute value of largest matrix element. If AMAX is very
56 * close to overflow or very close to underflow, the matrix
57 * should be scaled.
58 *
59 * INFO (output) INTEGER
60 * = 0: successful exit
61 * < 0: if INFO = -i, the i-th argument had an illegal value
62 * > 0: if INFO = i, the i-th diagonal element is nonpositive.
63 *
64 * =====================================================================
65 *
66 * .. Parameters ..
67 DOUBLE PRECISION ONE, ZERO
68 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 )
69 * ..
70 * .. Local Scalars ..
71 LOGICAL UPPER
72 INTEGER I, JJ
73 DOUBLE PRECISION SMIN
74 * ..
75 * .. External Functions ..
76 LOGICAL LSAME
77 EXTERNAL LSAME
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC MAX, MIN, SQRT
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 INFO = 0
90 UPPER = LSAME( UPLO, 'U' )
91 IF( .NOT.UPPER .AND. .NOT.LSAME( UPLO, 'L' ) ) THEN
92 INFO = -1
93 ELSE IF( N.LT.0 ) THEN
94 INFO = -2
95 END IF
96 IF( INFO.NE.0 ) THEN
97 CALL XERBLA( 'DPPEQU', -INFO )
98 RETURN
99 END IF
100 *
101 * Quick return if possible
102 *
103 IF( N.EQ.0 ) THEN
104 SCOND = ONE
105 AMAX = ZERO
106 RETURN
107 END IF
108 *
109 * Initialize SMIN and AMAX.
110 *
111 S( 1 ) = AP( 1 )
112 SMIN = S( 1 )
113 AMAX = S( 1 )
114 *
115 IF( UPPER ) THEN
116 *
117 * UPLO = 'U': Upper triangle of A is stored.
118 * Find the minimum and maximum diagonal elements.
119 *
120 JJ = 1
121 DO 10 I = 2, N
122 JJ = JJ + I
123 S( I ) = AP( JJ )
124 SMIN = MIN( SMIN, S( I ) )
125 AMAX = MAX( AMAX, S( I ) )
126 10 CONTINUE
127 *
128 ELSE
129 *
130 * UPLO = 'L': Lower triangle of A is stored.
131 * Find the minimum and maximum diagonal elements.
132 *
133 JJ = 1
134 DO 20 I = 2, N
135 JJ = JJ + N - I + 2
136 S( I ) = AP( JJ )
137 SMIN = MIN( SMIN, S( I ) )
138 AMAX = MAX( AMAX, S( I ) )
139 20 CONTINUE
140 END IF
141 *
142 IF( SMIN.LE.ZERO ) THEN
143 *
144 * Find the first non-positive diagonal element and return.
145 *
146 DO 30 I = 1, N
147 IF( S( I ).LE.ZERO ) THEN
148 INFO = I
149 RETURN
150 END IF
151 30 CONTINUE
152 ELSE
153 *
154 * Set the scale factors to the reciprocals
155 * of the diagonal elements.
156 *
157 DO 40 I = 1, N
158 S( I ) = ONE / SQRT( S( I ) )
159 40 CONTINUE
160 *
161 * Compute SCOND = min(S(I)) / max(S(I))
162 *
163 SCOND = SQRT( SMIN ) / SQRT( AMAX )
164 END IF
165 RETURN
166 *
167 * End of DPPEQU
168 *
169 END