 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2.2) -- 00004 * -- Contributed by James Demmel, Deaglan Halligan, Yozo Hida and -- 00005 * -- Jason Riedy of Univ. of California Berkeley. -- 00006 * -- June 2010 -- 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 CHARACTER UPLO 00017 * .. 00018 * .. Array Arguments .. 00019 COMPLEX*16 A( LDA, * ), WORK( * ) 00020 DOUBLE PRECISION S( * ) 00021 * .. 00022 * 00023 * Purpose 00024 * ======= 00025 * 00026 * ZSYEQUB computes row and column scalings intended to equilibrate a 00027 * symmetric matrix A and reduce its condition number 00028 * (with respect to the two-norm). S contains the scale factors, 00029 * S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with 00030 * elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This 00031 * choice of S puts the condition number of B within a factor N of the 00032 * smallest possible condition number over all possible diagonal 00033 * scalings. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * N (input) INTEGER 00039 * The order of the matrix A. N >= 0. 00040 * 00041 * A (input) COMPLEX*16 array, dimension (LDA,N) 00042 * The N-by-N symmetric matrix whose scaling 00043 * factors are to be computed. Only the diagonal elements of A 00044 * are referenced. 00045 * 00046 * LDA (input) INTEGER 00047 * The leading dimension of the array A. LDA >= max(1,N). 00048 * 00049 * S (output) DOUBLE PRECISION array, dimension (N) 00050 * If INFO = 0, S contains the scale factors for A. 00051 * 00052 * SCOND (output) DOUBLE PRECISION 00053 * If INFO = 0, S contains the ratio of the smallest S(i) to 00054 * the largest S(i). If SCOND >= 0.1 and AMAX is neither too 00055 * large nor too small, it is not worth scaling by S. 00056 * 00057 * AMAX (output) DOUBLE PRECISION 00058 * Absolute value of largest matrix element. If AMAX is very 00059 * close to overflow or very close to underflow, the matrix 00060 * should be scaled. 00061 * INFO (output) INTEGER 00062 * = 0: successful exit 00063 * < 0: if INFO = -i, the i-th argument had an illegal value 00064 * > 0: if INFO = i, the i-th diagonal element is nonpositive. 00065 * 00066 * ===================================================================== 00067 * 00068 * .. Parameters .. 00069 DOUBLE PRECISION ONE, ZERO 00070 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00071 INTEGER MAX_ITER 00072 PARAMETER ( MAX_ITER = 100 ) 00073 * .. 00074 * .. Local Scalars .. 00075 INTEGER I, J, ITER 00076 DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, 00077 $ BASE, SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ 00078 LOGICAL UP 00079 COMPLEX*16 ZDUM 00080 * .. 00081 * .. External Functions .. 00082 DOUBLE PRECISION DLAMCH 00083 LOGICAL LSAME 00084 EXTERNAL DLAMCH, LSAME 00085 * .. 00086 * .. External Subroutines .. 00087 EXTERNAL ZLASSQ 00088 * .. 00089 * .. Intrinsic Functions .. 00090 INTRINSIC ABS, DBLE, DIMAG, INT, LOG, MAX, MIN, SQRT 00091 * .. 00092 * .. Statement Functions .. 00093 DOUBLE PRECISION CABS1 00094 * .. 00095 * .. Statement Function Definitions .. 00096 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00097 * 00098 * Test input parameters. 00099 * 00100 INFO = 0 00101 IF (.NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN 00102 INFO = -1 00103 ELSE IF ( N .LT. 0 ) THEN 00104 INFO = -2 00105 ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN 00106 INFO = -4 00107 END IF 00108 IF ( INFO .NE. 0 ) THEN 00109 CALL XERBLA( 'ZHEEQUB', -INFO ) 00110 RETURN 00111 END IF 00112 00113 UP = LSAME( UPLO, 'U' ) 00114 AMAX = ZERO 00115 * 00116 * Quick return if possible. 00117 * 00118 IF ( N .EQ. 0 ) THEN 00119 SCOND = ONE 00120 RETURN 00121 END IF 00122 00123 DO I = 1, N 00124 S( I ) = ZERO 00125 END DO 00126 00127 AMAX = ZERO 00128 IF ( UP ) THEN 00129 DO J = 1, N 00130 DO I = 1, J-1 00131 S( I ) = MAX( S( I ), CABS1( A( I, J ) ) ) 00132 S( J ) = MAX( S( J ), CABS1( A( I, J ) ) ) 00133 AMAX = MAX( AMAX, CABS1( A( I, J ) ) ) 00134 END DO 00135 S( J ) = MAX( S( J ), CABS1( A( J, J ) ) ) 00136 AMAX = MAX( AMAX, CABS1( A( J, J ) ) ) 00137 END DO 00138 ELSE 00139 DO J = 1, N 00140 S( J ) = MAX( S( J ), CABS1( A( J, J ) ) ) 00141 AMAX = MAX( AMAX, CABS1( A( J, J ) ) ) 00142 DO I = J+1, N 00143 S( I ) = MAX( S( I ), CABS1( A( I, J ) ) ) 00144 S( J ) = MAX( S( J ), CABS1( A( I, J ) ) ) 00145 AMAX = MAX( AMAX, CABS1( A(I, J ) ) ) 00146 END DO 00147 END DO 00148 END IF 00149 DO J = 1, N 00150 S( J ) = 1.0D+0 / S( J ) 00151 END DO 00152 00153 TOL = ONE / SQRT( 2.0D0 * N ) 00154 00155 DO ITER = 1, MAX_ITER 00156 SCALE = 0.0D+0 00157 SUMSQ = 0.0D+0 00158 * beta = |A|s 00159 DO I = 1, N 00160 WORK( I ) = ZERO 00161 END DO 00162 IF ( UP ) THEN 00163 DO J = 1, N 00164 DO I = 1, J-1 00165 T = CABS1( A( I, J ) ) 00166 WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) 00167 WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) 00168 END DO 00169 WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) 00170 END DO 00171 ELSE 00172 DO J = 1, N 00173 WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) 00174 DO I = J+1, N 00175 T = CABS1( A( I, J ) ) 00176 WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) 00177 WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) 00178 END DO 00179 END DO 00180 END IF 00181 00182 * avg = s^T beta / n 00183 AVG = 0.0D+0 00184 DO I = 1, N 00185 AVG = AVG + S( I )*WORK( I ) 00186 END DO 00187 AVG = AVG / N 00188 00189 STD = 0.0D+0 00190 DO I = 2*N+1, 3*N 00191 WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG 00192 END DO 00193 CALL ZLASSQ( N, WORK( 2*N+1 ), 1, SCALE, SUMSQ ) 00194 STD = SCALE * SQRT( SUMSQ / N ) 00195 00196 IF ( STD .LT. TOL * AVG ) GOTO 999 00197 00198 DO I = 1, N 00199 T = CABS1( A( I, I ) ) 00200 SI = S( I ) 00201 C2 = ( N-1 ) * T 00202 C1 = ( N-2 ) * ( WORK( I ) - T*SI ) 00203 C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG 00204 00205 D = C1*C1 - 4*C0*C2 00206 IF ( D .LE. 0 ) THEN 00207 INFO = -1 00208 RETURN 00209 END IF 00210 SI = -2*C0 / ( C1 + SQRT( D ) ) 00211 00212 D = SI - S(I) 00213 U = ZERO 00214 IF ( UP ) THEN 00215 DO J = 1, I 00216 T = CABS1( A( J, I ) ) 00217 U = U + S( J )*T 00218 WORK( J ) = WORK( J ) + D*T 00219 END DO 00220 DO J = I+1,N 00221 T = CABS1( A( I, J ) ) 00222 U = U + S( J )*T 00223 WORK( J ) = WORK( J ) + D*T 00224 END DO 00225 ELSE 00226 DO J = 1, I 00227 T = CABS1( A( I, J ) ) 00228 U = U + S( J )*T 00229 WORK( J ) = WORK( J ) + D*T 00230 END DO 00231 DO J = I+1,N 00232 T = CABS1( A( J, I ) ) 00233 U = U + S( J )*T 00234 WORK( J ) = WORK( J ) + D*T 00235 END DO 00236 END IF 00237 AVG = AVG + ( U + WORK( I ) ) * D / N 00238 S( I ) = SI 00239 END DO 00240 00241 END DO 00242 00243 999 CONTINUE 00244 00245 SMLNUM = DLAMCH( 'SAFEMIN' ) 00246 BIGNUM = ONE / SMLNUM 00247 SMIN = BIGNUM 00248 SMAX = ZERO 00249 T = ONE / SQRT( AVG ) 00250 BASE = DLAMCH( 'B' ) 00251 U = ONE / LOG( BASE ) 00252 DO I = 1, N 00253 S( I ) = BASE ** INT( U * LOG( S( I ) * T ) ) 00254 SMIN = MIN( SMIN, S( I ) ) 00255 SMAX = MAX( SMAX, S( I ) ) 00256 END DO 00257 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 00258 00259 END
1.7.3 
