00001 SUBROUTINE DPOEQUB( N, A, LDA, S, SCOND, AMAX, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00005 * -- Jason Riedy of Univ. of California Berkeley. -- 00006 * -- November 2008 -- 00007 * 00008 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00009 * -- Univ. of California Berkeley and NAG Ltd. -- 00010 * 00011 IMPLICIT NONE 00012 * .. 00013 * .. Scalar Arguments .. 00014 INTEGER INFO, LDA, N 00015 DOUBLE PRECISION AMAX, SCOND 00016 * .. 00017 * .. Array Arguments .. 00018 DOUBLE PRECISION A( LDA, * ), S( * ) 00019 * .. 00020 * 00021 * Purpose 00022 * ======= 00023 * 00024 * DPOEQU computes row and column scalings intended to equilibrate a 00025 * symmetric positive definite matrix A and reduce its condition number 00026 * (with respect to the two-norm). S contains the scale factors, 00027 * S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with 00028 * elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This 00029 * choice of S puts the condition number of B within a factor N of the 00030 * smallest possible condition number over all possible diagonal 00031 * scalings. 00032 * 00033 * Arguments 00034 * ========= 00035 * 00036 * N (input) INTEGER 00037 * The order of the matrix A. N >= 0. 00038 * 00039 * A (input) DOUBLE PRECISION array, dimension (LDA,N) 00040 * The N-by-N symmetric positive definite matrix whose scaling 00041 * factors are to be computed. Only the diagonal elements of A 00042 * are referenced. 00043 * 00044 * LDA (input) INTEGER 00045 * The leading dimension of the array A. LDA >= max(1,N). 00046 * 00047 * S (output) DOUBLE PRECISION array, dimension (N) 00048 * If INFO = 0, S contains the scale factors for A. 00049 * 00050 * SCOND (output) DOUBLE PRECISION 00051 * If INFO = 0, S contains the ratio of the smallest S(i) to 00052 * the largest S(i). If SCOND >= 0.1 and AMAX is neither too 00053 * large nor too small, it is not worth scaling by S. 00054 * 00055 * AMAX (output) DOUBLE PRECISION 00056 * Absolute value of largest matrix element. If AMAX is very 00057 * close to overflow or very close to underflow, the matrix 00058 * should be scaled. 00059 * 00060 * INFO (output) INTEGER 00061 * = 0: successful exit 00062 * < 0: if INFO = -i, the i-th argument had an illegal value 00063 * > 0: if INFO = i, the i-th diagonal element is nonpositive. 00064 * 00065 * ===================================================================== 00066 * 00067 * .. Parameters .. 00068 DOUBLE PRECISION ZERO, ONE 00069 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00070 * .. 00071 * .. Local Scalars .. 00072 INTEGER I 00073 DOUBLE PRECISION SMIN, BASE, TMP 00074 * .. 00075 * .. External Functions .. 00076 DOUBLE PRECISION DLAMCH 00077 EXTERNAL DLAMCH 00078 * .. 00079 * .. External Subroutines .. 00080 EXTERNAL XERBLA 00081 * .. 00082 * .. Intrinsic Functions .. 00083 INTRINSIC MAX, MIN, SQRT, LOG, INT 00084 * .. 00085 * .. Executable Statements .. 00086 * 00087 * Test the input parameters. 00088 * 00089 * Positive definite only performs 1 pass of equilibration. 00090 * 00091 INFO = 0 00092 IF( N.LT.0 ) THEN 00093 INFO = -1 00094 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00095 INFO = -3 00096 END IF 00097 IF( INFO.NE.0 ) THEN 00098 CALL XERBLA( 'DPOEQUB', -INFO ) 00099 RETURN 00100 END IF 00101 * 00102 * Quick return if possible. 00103 * 00104 IF( N.EQ.0 ) THEN 00105 SCOND = ONE 00106 AMAX = ZERO 00107 RETURN 00108 END IF 00109 00110 BASE = DLAMCH( 'B' ) 00111 TMP = -0.5D+0 / LOG ( BASE ) 00112 * 00113 * Find the minimum and maximum diagonal elements. 00114 * 00115 S( 1 ) = A( 1, 1 ) 00116 SMIN = S( 1 ) 00117 AMAX = S( 1 ) 00118 DO 10 I = 2, N 00119 S( I ) = A( I, I ) 00120 SMIN = MIN( SMIN, S( I ) ) 00121 AMAX = MAX( AMAX, S( I ) ) 00122 10 CONTINUE 00123 * 00124 IF( SMIN.LE.ZERO ) THEN 00125 * 00126 * Find the first non-positive diagonal element and return. 00127 * 00128 DO 20 I = 1, N 00129 IF( S( I ).LE.ZERO ) THEN 00130 INFO = I 00131 RETURN 00132 END IF 00133 20 CONTINUE 00134 ELSE 00135 * 00136 * Set the scale factors to the reciprocals 00137 * of the diagonal elements. 00138 * 00139 DO 30 I = 1, N 00140 S( I ) = BASE ** INT( TMP * LOG( S( I ) ) ) 00141 30 CONTINUE 00142 * 00143 * Compute SCOND = min(S(I)) / max(S(I)). 00144 * 00145 SCOND = SQRT( SMIN ) / SQRT( AMAX ) 00146 END IF 00147 * 00148 RETURN 00149 * 00150 * End of DPOEQUB 00151 * 00152 END
1.7.2