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