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