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