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