1 SUBROUTINE DPOEQUB( N, A, LDA, S, SCOND, AMAX, INFO )
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
5 * -- Jason Riedy of Univ. of California Berkeley. --
6 * -- November 2008 --
7 *
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley and NAG Ltd. --
10 *
11 IMPLICIT NONE
12 * ..
13 * .. Scalar Arguments ..
14 INTEGER INFO, LDA, N
15 DOUBLE PRECISION AMAX, SCOND
16 * ..
17 * .. Array Arguments ..
18 DOUBLE PRECISION A( LDA, * ), S( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DPOEQU computes row and column scalings intended to equilibrate a
25 * symmetric positive definite matrix A and reduce its condition number
26 * (with respect to the two-norm). S contains the scale factors,
27 * S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
28 * elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This
29 * choice of S puts the condition number of B within a factor N of the
30 * smallest possible condition number over all possible diagonal
31 * scalings.
32 *
33 * Arguments
34 * =========
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
39 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
40 * The N-by-N symmetric positive definite matrix whose scaling
41 * factors are to be computed. Only the diagonal elements of A
42 * are referenced.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,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 ZERO, ONE
69 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
70 * ..
71 * .. Local Scalars ..
72 INTEGER I
73 DOUBLE PRECISION SMIN, BASE, TMP
74 * ..
75 * .. External Functions ..
76 DOUBLE PRECISION DLAMCH
77 EXTERNAL DLAMCH
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC MAX, MIN, SQRT, LOG, INT
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 * Positive definite only performs 1 pass of equilibration.
90 *
91 INFO = 0
92 IF( N.LT.0 ) THEN
93 INFO = -1
94 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
95 INFO = -3
96 END IF
97 IF( INFO.NE.0 ) THEN
98 CALL XERBLA( 'DPOEQUB', -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 BASE = DLAMCH( 'B' )
111 TMP = -0.5D+0 / LOG ( BASE )
112 *
113 * Find the minimum and maximum diagonal elements.
114 *
115 S( 1 ) = A( 1, 1 )
116 SMIN = S( 1 )
117 AMAX = S( 1 )
118 DO 10 I = 2, N
119 S( I ) = A( I, I )
120 SMIN = MIN( SMIN, S( I ) )
121 AMAX = MAX( AMAX, S( I ) )
122 10 CONTINUE
123 *
124 IF( SMIN.LE.ZERO ) THEN
125 *
126 * Find the first non-positive diagonal element and return.
127 *
128 DO 20 I = 1, N
129 IF( S( I ).LE.ZERO ) THEN
130 INFO = I
131 RETURN
132 END IF
133 20 CONTINUE
134 ELSE
135 *
136 * Set the scale factors to the reciprocals
137 * of the diagonal elements.
138 *
139 DO 30 I = 1, N
140 S( I ) = BASE ** INT( TMP * LOG( S( I ) ) )
141 30 CONTINUE
142 *
143 * Compute SCOND = min(S(I)) / max(S(I)).
144 *
145 SCOND = SQRT( SMIN ) / SQRT( AMAX )
146 END IF
147 *
148 RETURN
149 *
150 * End of DPOEQUB
151 *
152 END
2 *
3 * -- LAPACK routine (version 3.2) --
4 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and --
5 * -- Jason Riedy of Univ. of California Berkeley. --
6 * -- November 2008 --
7 *
8 * -- LAPACK is a software package provided by Univ. of Tennessee, --
9 * -- Univ. of California Berkeley and NAG Ltd. --
10 *
11 IMPLICIT NONE
12 * ..
13 * .. Scalar Arguments ..
14 INTEGER INFO, LDA, N
15 DOUBLE PRECISION AMAX, SCOND
16 * ..
17 * .. Array Arguments ..
18 DOUBLE PRECISION A( LDA, * ), S( * )
19 * ..
20 *
21 * Purpose
22 * =======
23 *
24 * DPOEQU computes row and column scalings intended to equilibrate a
25 * symmetric positive definite matrix A and reduce its condition number
26 * (with respect to the two-norm). S contains the scale factors,
27 * S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with
28 * elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This
29 * choice of S puts the condition number of B within a factor N of the
30 * smallest possible condition number over all possible diagonal
31 * scalings.
32 *
33 * Arguments
34 * =========
35 *
36 * N (input) INTEGER
37 * The order of the matrix A. N >= 0.
38 *
39 * A (input) DOUBLE PRECISION array, dimension (LDA,N)
40 * The N-by-N symmetric positive definite matrix whose scaling
41 * factors are to be computed. Only the diagonal elements of A
42 * are referenced.
43 *
44 * LDA (input) INTEGER
45 * The leading dimension of the array A. LDA >= max(1,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 ZERO, ONE
69 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
70 * ..
71 * .. Local Scalars ..
72 INTEGER I
73 DOUBLE PRECISION SMIN, BASE, TMP
74 * ..
75 * .. External Functions ..
76 DOUBLE PRECISION DLAMCH
77 EXTERNAL DLAMCH
78 * ..
79 * .. External Subroutines ..
80 EXTERNAL XERBLA
81 * ..
82 * .. Intrinsic Functions ..
83 INTRINSIC MAX, MIN, SQRT, LOG, INT
84 * ..
85 * .. Executable Statements ..
86 *
87 * Test the input parameters.
88 *
89 * Positive definite only performs 1 pass of equilibration.
90 *
91 INFO = 0
92 IF( N.LT.0 ) THEN
93 INFO = -1
94 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
95 INFO = -3
96 END IF
97 IF( INFO.NE.0 ) THEN
98 CALL XERBLA( 'DPOEQUB', -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 BASE = DLAMCH( 'B' )
111 TMP = -0.5D+0 / LOG ( BASE )
112 *
113 * Find the minimum and maximum diagonal elements.
114 *
115 S( 1 ) = A( 1, 1 )
116 SMIN = S( 1 )
117 AMAX = S( 1 )
118 DO 10 I = 2, N
119 S( I ) = A( I, I )
120 SMIN = MIN( SMIN, S( I ) )
121 AMAX = MAX( AMAX, S( I ) )
122 10 CONTINUE
123 *
124 IF( SMIN.LE.ZERO ) THEN
125 *
126 * Find the first non-positive diagonal element and return.
127 *
128 DO 20 I = 1, N
129 IF( S( I ).LE.ZERO ) THEN
130 INFO = I
131 RETURN
132 END IF
133 20 CONTINUE
134 ELSE
135 *
136 * Set the scale factors to the reciprocals
137 * of the diagonal elements.
138 *
139 DO 30 I = 1, N
140 S( I ) = BASE ** INT( TMP * LOG( S( I ) ) )
141 30 CONTINUE
142 *
143 * Compute SCOND = min(S(I)) / max(S(I)).
144 *
145 SCOND = SQRT( SMIN ) / SQRT( AMAX )
146 END IF
147 *
148 RETURN
149 *
150 * End of DPOEQUB
151 *
152 END