![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b ZHEEQUB 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download ZHEEQUB + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zheequb.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zheequb.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zheequb.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER INFO, LDA, N 00025 * DOUBLE PRECISION AMAX, SCOND 00026 * CHARACTER UPLO 00027 * .. 00028 * .. Array Arguments .. 00029 * COMPLEX*16 A( LDA, * ), WORK( * ) 00030 * DOUBLE PRECISION S( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> ZHEEQUB computes row and column scalings intended to equilibrate a 00040 *> Hermitian matrix A and reduce its condition number 00041 *> (with respect to the two-norm). S contains the scale factors, 00042 *> S(i) = 1/sqrt(A(i,i)), chosen so that the scaled matrix B with 00043 *> elements B(i,j) = S(i)*A(i,j)*S(j) has ones on the diagonal. This 00044 *> choice of S puts the condition number of B within a factor N of the 00045 *> smallest possible condition number over all possible diagonal 00046 *> scalings. 00047 *> \endverbatim 00048 * 00049 * Arguments: 00050 * ========== 00051 * 00052 *> \param[in] UPLO 00053 *> \verbatim 00054 *> UPLO is CHARACTER*1 00055 *> = 'U': Upper triangles of A and B are stored; 00056 *> = 'L': Lower triangles of A and B are stored. 00057 *> \endverbatim 00058 *> 00059 *> \param[in] N 00060 *> \verbatim 00061 *> N is INTEGER 00062 *> The order of the matrix A. N >= 0. 00063 *> \endverbatim 00064 *> 00065 *> \param[in] A 00066 *> \verbatim 00067 *> A is COMPLEX*16 array, dimension (LDA,N) 00068 *> The N-by-N Hermitian matrix whose scaling 00069 *> factors are to be computed. Only the diagonal elements of A 00070 *> are referenced. 00071 *> \endverbatim 00072 *> 00073 *> \param[in] LDA 00074 *> \verbatim 00075 *> LDA is INTEGER 00076 *> The leading dimension of the array A. LDA >= max(1,N). 00077 *> \endverbatim 00078 *> 00079 *> \param[out] S 00080 *> \verbatim 00081 *> S is DOUBLE PRECISION array, dimension (N) 00082 *> If INFO = 0, S contains the scale factors for A. 00083 *> \endverbatim 00084 *> 00085 *> \param[out] SCOND 00086 *> \verbatim 00087 *> SCOND is DOUBLE PRECISION 00088 *> If INFO = 0, S contains the ratio of the smallest S(i) to 00089 *> the largest S(i). If SCOND >= 0.1 and AMAX is neither too 00090 *> large nor too small, it is not worth scaling by S. 00091 *> \endverbatim 00092 *> 00093 *> \param[out] AMAX 00094 *> \verbatim 00095 *> AMAX is DOUBLE PRECISION 00096 *> Absolute value of largest matrix element. If AMAX is very 00097 *> close to overflow or very close to underflow, the matrix 00098 *> should be scaled. 00099 *> \endverbatim 00100 *> 00101 *> \param[out] WORK 00102 *> \verbatim 00103 *> WORK is DOUBLE PRECISION array, dimension (3*N) 00104 *> \endverbatim 00105 *> 00106 *> \param[out] INFO 00107 *> \verbatim 00108 *> INFO is INTEGER 00109 *> = 0: successful exit 00110 *> < 0: if INFO = -i, the i-th argument had an illegal value 00111 *> > 0: if INFO = i, the i-th diagonal element is nonpositive. 00112 *> \endverbatim 00113 * 00114 * Authors: 00115 * ======== 00116 * 00117 *> \author Univ. of Tennessee 00118 *> \author Univ. of California Berkeley 00119 *> \author Univ. of Colorado Denver 00120 *> \author NAG Ltd. 00121 * 00122 *> \date April 2012 00123 * 00124 *> \ingroup complex16HEcomputational 00125 * 00126 * ===================================================================== 00127 SUBROUTINE ZHEEQUB( UPLO, N, A, LDA, S, SCOND, AMAX, WORK, INFO ) 00128 * 00129 * -- LAPACK computational routine (version 3.4.1) -- 00130 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00131 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00132 * April 2012 00133 * 00134 * .. Scalar Arguments .. 00135 INTEGER INFO, LDA, N 00136 DOUBLE PRECISION AMAX, SCOND 00137 CHARACTER UPLO 00138 * .. 00139 * .. Array Arguments .. 00140 COMPLEX*16 A( LDA, * ), WORK( * ) 00141 DOUBLE PRECISION S( * ) 00142 * .. 00143 * 00144 * ===================================================================== 00145 * 00146 * .. Parameters .. 00147 DOUBLE PRECISION ONE, ZERO 00148 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00149 INTEGER MAX_ITER 00150 PARAMETER ( MAX_ITER = 100 ) 00151 * .. 00152 * .. Local Scalars .. 00153 INTEGER I, J, ITER 00154 DOUBLE PRECISION AVG, STD, TOL, C0, C1, C2, T, U, SI, D, 00155 $ BASE, SMIN, SMAX, SMLNUM, BIGNUM, SCALE, SUMSQ 00156 LOGICAL UP 00157 COMPLEX*16 ZDUM 00158 * .. 00159 * .. External Functions .. 00160 DOUBLE PRECISION DLAMCH 00161 LOGICAL LSAME 00162 EXTERNAL DLAMCH, LSAME 00163 * .. 00164 * .. External Subroutines .. 00165 EXTERNAL ZLASSQ 00166 * .. 00167 * .. Intrinsic Functions .. 00168 INTRINSIC ABS, DBLE, DIMAG, INT, LOG, MAX, MIN, SQRT 00169 * .. 00170 * .. Statement Functions .. 00171 DOUBLE PRECISION CABS1 00172 * .. 00173 * .. Statement Function Definitions .. 00174 CABS1( ZDUM ) = ABS( DBLE( ZDUM ) ) + ABS( DIMAG( ZDUM ) ) 00175 * 00176 * Test input parameters. 00177 * 00178 INFO = 0 00179 IF (.NOT. ( LSAME( UPLO, 'U' ) .OR. LSAME( UPLO, 'L' ) ) ) THEN 00180 INFO = -1 00181 ELSE IF ( N .LT. 0 ) THEN 00182 INFO = -2 00183 ELSE IF ( LDA .LT. MAX( 1, N ) ) THEN 00184 INFO = -4 00185 END IF 00186 IF ( INFO .NE. 0 ) THEN 00187 CALL XERBLA( 'ZHEEQUB', -INFO ) 00188 RETURN 00189 END IF 00190 00191 UP = LSAME( UPLO, 'U' ) 00192 AMAX = ZERO 00193 * 00194 * Quick return if possible. 00195 * 00196 IF ( N .EQ. 0 ) THEN 00197 SCOND = ONE 00198 RETURN 00199 END IF 00200 00201 DO I = 1, N 00202 S( I ) = ZERO 00203 END DO 00204 00205 AMAX = ZERO 00206 IF ( UP ) THEN 00207 DO J = 1, N 00208 DO I = 1, J-1 00209 S( I ) = MAX( S( I ), CABS1( A( I, J ) ) ) 00210 S( J ) = MAX( S( J ), CABS1( A( I, J ) ) ) 00211 AMAX = MAX( AMAX, CABS1( A( I, J ) ) ) 00212 END DO 00213 S( J ) = MAX( S( J ), CABS1( A( J, J ) ) ) 00214 AMAX = MAX( AMAX, CABS1( A( J, J ) ) ) 00215 END DO 00216 ELSE 00217 DO J = 1, N 00218 S( J ) = MAX( S( J ), CABS1( A( J, J ) ) ) 00219 AMAX = MAX( AMAX, CABS1( A( J, J ) ) ) 00220 DO I = J+1, N 00221 S( I ) = MAX( S( I ), CABS1( A( I, J ) ) ) 00222 S( J ) = MAX( S( J ), CABS1( A( I, J ) ) ) 00223 AMAX = MAX( AMAX, CABS1( A(I, J ) ) ) 00224 END DO 00225 END DO 00226 END IF 00227 DO J = 1, N 00228 S( J ) = 1.0D+0 / S( J ) 00229 END DO 00230 00231 TOL = ONE / SQRT( 2.0D0 * N ) 00232 00233 DO ITER = 1, MAX_ITER 00234 SCALE = 0.0D+0 00235 SUMSQ = 0.0D+0 00236 * beta = |A|s 00237 DO I = 1, N 00238 WORK( I ) = ZERO 00239 END DO 00240 IF ( UP ) THEN 00241 DO J = 1, N 00242 DO I = 1, J-1 00243 T = CABS1( A( I, J ) ) 00244 WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) 00245 WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) 00246 END DO 00247 WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) 00248 END DO 00249 ELSE 00250 DO J = 1, N 00251 WORK( J ) = WORK( J ) + CABS1( A( J, J ) ) * S( J ) 00252 DO I = J+1, N 00253 T = CABS1( A( I, J ) ) 00254 WORK( I ) = WORK( I ) + CABS1( A( I, J ) ) * S( J ) 00255 WORK( J ) = WORK( J ) + CABS1( A( I, J ) ) * S( I ) 00256 END DO 00257 END DO 00258 END IF 00259 00260 * avg = s^T beta / n 00261 AVG = 0.0D+0 00262 DO I = 1, N 00263 AVG = AVG + S( I )*WORK( I ) 00264 END DO 00265 AVG = AVG / N 00266 00267 STD = 0.0D+0 00268 DO I = 2*N+1, 3*N 00269 WORK( I ) = S( I-2*N ) * WORK( I-2*N ) - AVG 00270 END DO 00271 CALL ZLASSQ( N, WORK( 2*N+1 ), 1, SCALE, SUMSQ ) 00272 STD = SCALE * SQRT( SUMSQ / N ) 00273 00274 IF ( STD .LT. TOL * AVG ) GOTO 999 00275 00276 DO I = 1, N 00277 T = CABS1( A( I, I ) ) 00278 SI = S( I ) 00279 C2 = ( N-1 ) * T 00280 C1 = ( N-2 ) * ( WORK( I ) - T*SI ) 00281 C0 = -(T*SI)*SI + 2*WORK( I )*SI - N*AVG 00282 00283 D = C1*C1 - 4*C0*C2 00284 IF ( D .LE. 0 ) THEN 00285 INFO = -1 00286 RETURN 00287 END IF 00288 SI = -2*C0 / ( C1 + SQRT( D ) ) 00289 00290 D = SI - S(I) 00291 U = ZERO 00292 IF ( UP ) THEN 00293 DO J = 1, I 00294 T = CABS1( A( J, I ) ) 00295 U = U + S( J )*T 00296 WORK( J ) = WORK( J ) + D*T 00297 END DO 00298 DO J = I+1,N 00299 T = CABS1( A( I, J ) ) 00300 U = U + S( J )*T 00301 WORK( J ) = WORK( J ) + D*T 00302 END DO 00303 ELSE 00304 DO J = 1, I 00305 T = CABS1( A( I, J ) ) 00306 U = U + S( J )*T 00307 WORK( J ) = WORK( J ) + D*T 00308 END DO 00309 DO J = I+1,N 00310 T = CABS1( A( J, I ) ) 00311 U = U + S( J )*T 00312 WORK( J ) = WORK( J ) + D*T 00313 END DO 00314 END IF 00315 AVG = AVG + ( U + WORK( I ) ) * D / N 00316 S( I ) = SI 00317 END DO 00318 00319 END DO 00320 00321 999 CONTINUE 00322 00323 SMLNUM = DLAMCH( 'SAFEMIN' ) 00324 BIGNUM = ONE / SMLNUM 00325 SMIN = BIGNUM 00326 SMAX = ZERO 00327 T = ONE / SQRT( AVG ) 00328 BASE = DLAMCH( 'B' ) 00329 U = ONE / LOG( BASE ) 00330 DO I = 1, N 00331 S( I ) = BASE ** INT( U * LOG( S( I ) * T ) ) 00332 SMIN = MIN( SMIN, S( I ) ) 00333 SMAX = MAX( SMAX, S( I ) ) 00334 END DO 00335 SCOND = MAX( SMIN, SMLNUM ) / MIN( SMAX, BIGNUM ) 00336 00337 END