 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASQ2( N, Z, INFO ) 00002 * 00003 * -- LAPACK routine (version 3.2) -- 00004 * 00005 * -- Contributed by Osni Marques of the Lawrence Berkeley National -- 00006 * -- Laboratory and Beresford Parlett of the Univ. of California at -- 00007 * -- Berkeley -- 00008 * -- November 2008 -- 00009 * 00010 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00011 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00012 * 00013 * .. Scalar Arguments .. 00014 INTEGER INFO, N 00015 * .. 00016 * .. Array Arguments .. 00017 DOUBLE PRECISION Z( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * DLASQ2 computes all the eigenvalues of the symmetric positive 00024 * definite tridiagonal matrix associated with the qd array Z to high 00025 * relative accuracy are computed to high relative accuracy, in the 00026 * absence of denormalization, underflow and overflow. 00027 * 00028 * To see the relation of Z to the tridiagonal matrix, let L be a 00029 * unit lower bidiagonal matrix with subdiagonals Z(2,4,6,,..) and 00030 * let U be an upper bidiagonal matrix with 1's above and diagonal 00031 * Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the 00032 * symmetric tridiagonal to which it is similar. 00033 * 00034 * Note : DLASQ2 defines a logical variable, IEEE, which is true 00035 * on machines which follow ieee-754 floating-point standard in their 00036 * handling of infinities and NaNs, and false otherwise. This variable 00037 * is passed to DLASQ3. 00038 * 00039 * Arguments 00040 * ========= 00041 * 00042 * N (input) INTEGER 00043 * The number of rows and columns in the matrix. N >= 0. 00044 * 00045 * Z (input/output) DOUBLE PRECISION array, dimension ( 4*N ) 00046 * On entry Z holds the qd array. On exit, entries 1 to N hold 00047 * the eigenvalues in decreasing order, Z( 2*N+1 ) holds the 00048 * trace, and Z( 2*N+2 ) holds the sum of the eigenvalues. If 00049 * N > 2, then Z( 2*N+3 ) holds the iteration count, Z( 2*N+4 ) 00050 * holds NDIVS/NIN^2, and Z( 2*N+5 ) holds the percentage of 00051 * shifts that failed. 00052 * 00053 * INFO (output) INTEGER 00054 * = 0: successful exit 00055 * < 0: if the i-th argument is a scalar and had an illegal 00056 * value, then INFO = -i, if the i-th argument is an 00057 * array and the j-entry had an illegal value, then 00058 * INFO = -(i*100+j) 00059 * > 0: the algorithm failed 00060 * = 1, a split was marked by a positive value in E 00061 * = 2, current block of Z not diagonalized after 30*N 00062 * iterations (in inner while loop) 00063 * = 3, termination criterion of outer while loop not met 00064 * (program created more than N unreduced blocks) 00065 * 00066 * Further Details 00067 * =============== 00068 * Local Variables: I0:N0 defines a current unreduced segment of Z. 00069 * The shifts are accumulated in SIGMA. Iteration count is in ITER. 00070 * Ping-pong is controlled by PP (alternates between 0 and 1). 00071 * 00072 * ===================================================================== 00073 * 00074 * .. Parameters .. 00075 DOUBLE PRECISION CBIAS 00076 PARAMETER ( CBIAS = 1.50D0 ) 00077 DOUBLE PRECISION ZERO, HALF, ONE, TWO, FOUR, HUNDRD 00078 PARAMETER ( ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, 00079 $ TWO = 2.0D0, FOUR = 4.0D0, HUNDRD = 100.0D0 ) 00080 * .. 00081 * .. Local Scalars .. 00082 LOGICAL IEEE 00083 INTEGER I0, I4, IINFO, IPN4, ITER, IWHILA, IWHILB, K, 00084 $ KMIN, N0, NBIG, NDIV, NFAIL, PP, SPLT, TTYPE 00085 DOUBLE PRECISION D, DEE, DEEMIN, DESIG, DMIN, DMIN1, DMIN2, DN, 00086 $ DN1, DN2, E, EMAX, EMIN, EPS, G, OLDEMN, QMAX, 00087 $ QMIN, S, SAFMIN, SIGMA, T, TAU, TEMP, TOL, 00088 $ TOL2, TRACE, ZMAX 00089 * .. 00090 * .. External Subroutines .. 00091 EXTERNAL DLASQ3, DLASRT, XERBLA 00092 * .. 00093 * .. External Functions .. 00094 INTEGER ILAENV 00095 DOUBLE PRECISION DLAMCH 00096 EXTERNAL DLAMCH, ILAENV 00097 * .. 00098 * .. Intrinsic Functions .. 00099 INTRINSIC ABS, DBLE, MAX, MIN, SQRT 00100 * .. 00101 * .. Executable Statements .. 00102 * 00103 * Test the input arguments. 00104 * (in case DLASQ2 is not called by DLASQ1) 00105 * 00106 INFO = 0 00107 EPS = DLAMCH( 'Precision' ) 00108 SAFMIN = DLAMCH( 'Safe minimum' ) 00109 TOL = EPS*HUNDRD 00110 TOL2 = TOL**2 00111 * 00112 IF( N.LT.0 ) THEN 00113 INFO = -1 00114 CALL XERBLA( 'DLASQ2', 1 ) 00115 RETURN 00116 ELSE IF( N.EQ.0 ) THEN 00117 RETURN 00118 ELSE IF( N.EQ.1 ) THEN 00119 * 00120 * 1-by-1 case. 00121 * 00122 IF( Z( 1 ).LT.ZERO ) THEN 00123 INFO = -201 00124 CALL XERBLA( 'DLASQ2', 2 ) 00125 END IF 00126 RETURN 00127 ELSE IF( N.EQ.2 ) THEN 00128 * 00129 * 2-by-2 case. 00130 * 00131 IF( Z( 2 ).LT.ZERO .OR. Z( 3 ).LT.ZERO ) THEN 00132 INFO = -2 00133 CALL XERBLA( 'DLASQ2', 2 ) 00134 RETURN 00135 ELSE IF( Z( 3 ).GT.Z( 1 ) ) THEN 00136 D = Z( 3 ) 00137 Z( 3 ) = Z( 1 ) 00138 Z( 1 ) = D 00139 END IF 00140 Z( 5 ) = Z( 1 ) + Z( 2 ) + Z( 3 ) 00141 IF( Z( 2 ).GT.Z( 3 )*TOL2 ) THEN 00142 T = HALF*( ( Z( 1 )-Z( 3 ) )+Z( 2 ) ) 00143 S = Z( 3 )*( Z( 2 ) / T ) 00144 IF( S.LE.T ) THEN 00145 S = Z( 3 )*( Z( 2 ) / ( T*( ONE+SQRT( ONE+S / T ) ) ) ) 00146 ELSE 00147 S = Z( 3 )*( Z( 2 ) / ( T+SQRT( T )*SQRT( T+S ) ) ) 00148 END IF 00149 T = Z( 1 ) + ( S+Z( 2 ) ) 00150 Z( 3 ) = Z( 3 )*( Z( 1 ) / T ) 00151 Z( 1 ) = T 00152 END IF 00153 Z( 2 ) = Z( 3 ) 00154 Z( 6 ) = Z( 2 ) + Z( 1 ) 00155 RETURN 00156 END IF 00157 * 00158 * Check for negative data and compute sums of q's and e's. 00159 * 00160 Z( 2*N ) = ZERO 00161 EMIN = Z( 2 ) 00162 QMAX = ZERO 00163 ZMAX = ZERO 00164 D = ZERO 00165 E = ZERO 00166 * 00167 DO 10 K = 1, 2*( N-1 ), 2 00168 IF( Z( K ).LT.ZERO ) THEN 00169 INFO = -( 200+K ) 00170 CALL XERBLA( 'DLASQ2', 2 ) 00171 RETURN 00172 ELSE IF( Z( K+1 ).LT.ZERO ) THEN 00173 INFO = -( 200+K+1 ) 00174 CALL XERBLA( 'DLASQ2', 2 ) 00175 RETURN 00176 END IF 00177 D = D + Z( K ) 00178 E = E + Z( K+1 ) 00179 QMAX = MAX( QMAX, Z( K ) ) 00180 EMIN = MIN( EMIN, Z( K+1 ) ) 00181 ZMAX = MAX( QMAX, ZMAX, Z( K+1 ) ) 00182 10 CONTINUE 00183 IF( Z( 2*N-1 ).LT.ZERO ) THEN 00184 INFO = -( 200+2*N-1 ) 00185 CALL XERBLA( 'DLASQ2', 2 ) 00186 RETURN 00187 END IF 00188 D = D + Z( 2*N-1 ) 00189 QMAX = MAX( QMAX, Z( 2*N-1 ) ) 00190 ZMAX = MAX( QMAX, ZMAX ) 00191 * 00192 * Check for diagonality. 00193 * 00194 IF( E.EQ.ZERO ) THEN 00195 DO 20 K = 2, N 00196 Z( K ) = Z( 2*K-1 ) 00197 20 CONTINUE 00198 CALL DLASRT( 'D', N, Z, IINFO ) 00199 Z( 2*N-1 ) = D 00200 RETURN 00201 END IF 00202 * 00203 TRACE = D + E 00204 * 00205 * Check for zero data. 00206 * 00207 IF( TRACE.EQ.ZERO ) THEN 00208 Z( 2*N-1 ) = ZERO 00209 RETURN 00210 END IF 00211 * 00212 * Check whether the machine is IEEE conformable. 00213 * 00214 IEEE = ILAENV( 10, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. 00215 $ ILAENV( 11, 'DLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 00216 * 00217 * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). 00218 * 00219 DO 30 K = 2*N, 2, -2 00220 Z( 2*K ) = ZERO 00221 Z( 2*K-1 ) = Z( K ) 00222 Z( 2*K-2 ) = ZERO 00223 Z( 2*K-3 ) = Z( K-1 ) 00224 30 CONTINUE 00225 * 00226 I0 = 1 00227 N0 = N 00228 * 00229 * Reverse the qd-array, if warranted. 00230 * 00231 IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN 00232 IPN4 = 4*( I0+N0 ) 00233 DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4 00234 TEMP = Z( I4-3 ) 00235 Z( I4-3 ) = Z( IPN4-I4-3 ) 00236 Z( IPN4-I4-3 ) = TEMP 00237 TEMP = Z( I4-1 ) 00238 Z( I4-1 ) = Z( IPN4-I4-5 ) 00239 Z( IPN4-I4-5 ) = TEMP 00240 40 CONTINUE 00241 END IF 00242 * 00243 * Initial split checking via dqd and Li's test. 00244 * 00245 PP = 0 00246 * 00247 DO 80 K = 1, 2 00248 * 00249 D = Z( 4*N0+PP-3 ) 00250 DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 00251 IF( Z( I4-1 ).LE.TOL2*D ) THEN 00252 Z( I4-1 ) = -ZERO 00253 D = Z( I4-3 ) 00254 ELSE 00255 D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) 00256 END IF 00257 50 CONTINUE 00258 * 00259 * dqd maps Z to ZZ plus Li's test. 00260 * 00261 EMIN = Z( 4*I0+PP+1 ) 00262 D = Z( 4*I0+PP-3 ) 00263 DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 00264 Z( I4-2*PP-2 ) = D + Z( I4-1 ) 00265 IF( Z( I4-1 ).LE.TOL2*D ) THEN 00266 Z( I4-1 ) = -ZERO 00267 Z( I4-2*PP-2 ) = D 00268 Z( I4-2*PP ) = ZERO 00269 D = Z( I4+1 ) 00270 ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. 00271 $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN 00272 TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) 00273 Z( I4-2*PP ) = Z( I4-1 )*TEMP 00274 D = D*TEMP 00275 ELSE 00276 Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) 00277 D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) 00278 END IF 00279 EMIN = MIN( EMIN, Z( I4-2*PP ) ) 00280 60 CONTINUE 00281 Z( 4*N0-PP-2 ) = D 00282 * 00283 * Now find qmax. 00284 * 00285 QMAX = Z( 4*I0-PP-2 ) 00286 DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4 00287 QMAX = MAX( QMAX, Z( I4 ) ) 00288 70 CONTINUE 00289 * 00290 * Prepare for the next iteration on K. 00291 * 00292 PP = 1 - PP 00293 80 CONTINUE 00294 * 00295 * Initialise variables to pass to DLASQ3. 00296 * 00297 TTYPE = 0 00298 DMIN1 = ZERO 00299 DMIN2 = ZERO 00300 DN = ZERO 00301 DN1 = ZERO 00302 DN2 = ZERO 00303 G = ZERO 00304 TAU = ZERO 00305 * 00306 ITER = 2 00307 NFAIL = 0 00308 NDIV = 2*( N0-I0 ) 00309 * 00310 DO 160 IWHILA = 1, N + 1 00311 IF( N0.LT.1 ) 00312 $ GO TO 170 00313 * 00314 * While array unfinished do 00315 * 00316 * E(N0) holds the value of SIGMA when submatrix in I0:N0 00317 * splits from the rest of the array, but is negated. 00318 * 00319 DESIG = ZERO 00320 IF( N0.EQ.N ) THEN 00321 SIGMA = ZERO 00322 ELSE 00323 SIGMA = -Z( 4*N0-1 ) 00324 END IF 00325 IF( SIGMA.LT.ZERO ) THEN 00326 INFO = 1 00327 RETURN 00328 END IF 00329 * 00330 * Find last unreduced submatrix's top index I0, find QMAX and 00331 * EMIN. Find Gershgorin-type bound if Q's much greater than E's. 00332 * 00333 EMAX = ZERO 00334 IF( N0.GT.I0 ) THEN 00335 EMIN = ABS( Z( 4*N0-5 ) ) 00336 ELSE 00337 EMIN = ZERO 00338 END IF 00339 QMIN = Z( 4*N0-3 ) 00340 QMAX = QMIN 00341 DO 90 I4 = 4*N0, 8, -4 00342 IF( Z( I4-5 ).LE.ZERO ) 00343 $ GO TO 100 00344 IF( QMIN.GE.FOUR*EMAX ) THEN 00345 QMIN = MIN( QMIN, Z( I4-3 ) ) 00346 EMAX = MAX( EMAX, Z( I4-5 ) ) 00347 END IF 00348 QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) 00349 EMIN = MIN( EMIN, Z( I4-5 ) ) 00350 90 CONTINUE 00351 I4 = 4 00352 * 00353 100 CONTINUE 00354 I0 = I4 / 4 00355 PP = 0 00356 * 00357 IF( N0-I0.GT.1 ) THEN 00358 DEE = Z( 4*I0-3 ) 00359 DEEMIN = DEE 00360 KMIN = I0 00361 DO 110 I4 = 4*I0+1, 4*N0-3, 4 00362 DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) ) 00363 IF( DEE.LE.DEEMIN ) THEN 00364 DEEMIN = DEE 00365 KMIN = ( I4+3 )/4 00366 END IF 00367 110 CONTINUE 00368 IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 00369 $ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN 00370 IPN4 = 4*( I0+N0 ) 00371 PP = 2 00372 DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4 00373 TEMP = Z( I4-3 ) 00374 Z( I4-3 ) = Z( IPN4-I4-3 ) 00375 Z( IPN4-I4-3 ) = TEMP 00376 TEMP = Z( I4-2 ) 00377 Z( I4-2 ) = Z( IPN4-I4-2 ) 00378 Z( IPN4-I4-2 ) = TEMP 00379 TEMP = Z( I4-1 ) 00380 Z( I4-1 ) = Z( IPN4-I4-5 ) 00381 Z( IPN4-I4-5 ) = TEMP 00382 TEMP = Z( I4 ) 00383 Z( I4 ) = Z( IPN4-I4-4 ) 00384 Z( IPN4-I4-4 ) = TEMP 00385 120 CONTINUE 00386 END IF 00387 END IF 00388 * 00389 * Put -(initial shift) into DMIN. 00390 * 00391 DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) 00392 * 00393 * Now I0:N0 is unreduced. 00394 * PP = 0 for ping, PP = 1 for pong. 00395 * PP = 2 indicates that flipping was applied to the Z array and 00396 * and that the tests for deflation upon entry in DLASQ3 00397 * should not be performed. 00398 * 00399 NBIG = 30*( N0-I0+1 ) 00400 DO 140 IWHILB = 1, NBIG 00401 IF( I0.GT.N0 ) 00402 $ GO TO 150 00403 * 00404 * While submatrix unfinished take a good dqds step. 00405 * 00406 CALL DLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 00407 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 00408 $ DN2, G, TAU ) 00409 * 00410 PP = 1 - PP 00411 * 00412 * When EMIN is very small check for splits. 00413 * 00414 IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN 00415 IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. 00416 $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN 00417 SPLT = I0 - 1 00418 QMAX = Z( 4*I0-3 ) 00419 EMIN = Z( 4*I0-1 ) 00420 OLDEMN = Z( 4*I0 ) 00421 DO 130 I4 = 4*I0, 4*( N0-3 ), 4 00422 IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. 00423 $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN 00424 Z( I4-1 ) = -SIGMA 00425 SPLT = I4 / 4 00426 QMAX = ZERO 00427 EMIN = Z( I4+3 ) 00428 OLDEMN = Z( I4+4 ) 00429 ELSE 00430 QMAX = MAX( QMAX, Z( I4+1 ) ) 00431 EMIN = MIN( EMIN, Z( I4-1 ) ) 00432 OLDEMN = MIN( OLDEMN, Z( I4 ) ) 00433 END IF 00434 130 CONTINUE 00435 Z( 4*N0-1 ) = EMIN 00436 Z( 4*N0 ) = OLDEMN 00437 I0 = SPLT + 1 00438 END IF 00439 END IF 00440 * 00441 140 CONTINUE 00442 * 00443 INFO = 2 00444 RETURN 00445 * 00446 * end IWHILB 00447 * 00448 150 CONTINUE 00449 * 00450 160 CONTINUE 00451 * 00452 INFO = 3 00453 RETURN 00454 * 00455 * end IWHILA 00456 * 00457 170 CONTINUE 00458 * 00459 * Move q's to the front. 00460 * 00461 DO 180 K = 2, N 00462 Z( K ) = Z( 4*K-3 ) 00463 180 CONTINUE 00464 * 00465 * Sort and compute sum of eigenvalues. 00466 * 00467 CALL DLASRT( 'D', N, Z, IINFO ) 00468 * 00469 E = ZERO 00470 DO 190 K = N, 1, -1 00471 E = E + Z( K ) 00472 190 CONTINUE 00473 * 00474 * Store trace, sum(eigenvalues) and information on performance. 00475 * 00476 Z( 2*N+1 ) = TRACE 00477 Z( 2*N+2 ) = E 00478 Z( 2*N+3 ) = DBLE( ITER ) 00479 Z( 2*N+4 ) = DBLE( NDIV ) / DBLE( N**2 ) 00480 Z( 2*N+5 ) = HUNDRD*NFAIL / DBLE( ITER ) 00481 RETURN 00482 * 00483 * End of DLASQ2 00484 * 00485 END
1.7.3 
