 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLASQ2( 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 REAL Z( * ) 00018 * .. 00019 * 00020 * Purpose 00021 * ======= 00022 * 00023 * SLASQ2 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 : SLASQ2 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 SLASQ3. 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) REAL 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 REAL CBIAS 00076 PARAMETER ( CBIAS = 1.50E0 ) 00077 REAL ZERO, HALF, ONE, TWO, FOUR, HUNDRD 00078 PARAMETER ( ZERO = 0.0E0, HALF = 0.5E0, ONE = 1.0E0, 00079 $ TWO = 2.0E0, FOUR = 4.0E0, HUNDRD = 100.0E0 ) 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 REAL 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 SLASQ3, SLASRT, XERBLA 00092 * .. 00093 * .. External Functions .. 00094 INTEGER ILAENV 00095 REAL SLAMCH 00096 EXTERNAL ILAENV, SLAMCH 00097 * .. 00098 * .. Intrinsic Functions .. 00099 INTRINSIC ABS, MAX, MIN, REAL, SQRT 00100 * .. 00101 * .. Executable Statements .. 00102 * 00103 * Test the input arguments. 00104 * (in case SLASQ2 is not called by SLASQ1) 00105 * 00106 INFO = 0 00107 EPS = SLAMCH( 'Precision' ) 00108 SAFMIN = SLAMCH( 'Safe minimum' ) 00109 TOL = EPS*HUNDRD 00110 TOL2 = TOL**2 00111 * 00112 IF( N.LT.0 ) THEN 00113 INFO = -1 00114 CALL XERBLA( 'SLASQ2', 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( 'SLASQ2', 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( 'SLASQ2', 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( 'SLASQ2', 2 ) 00171 RETURN 00172 ELSE IF( Z( K+1 ).LT.ZERO ) THEN 00173 INFO = -( 200+K+1 ) 00174 CALL XERBLA( 'SLASQ2', 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( 'SLASQ2', 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 SLASRT( '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, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 .AND. 00215 * $ ILAENV( 11, 'SLASQ2', 'N', 1, 2, 3, 4 ).EQ.1 00216 * 00217 * [11/15/2008] The case IEEE=.TRUE. has a problem in single precision with 00218 * some the test matrices of type 16. The double precision code is fine. 00219 * 00220 IEEE = .FALSE. 00221 * 00222 * Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...). 00223 * 00224 DO 30 K = 2*N, 2, -2 00225 Z( 2*K ) = ZERO 00226 Z( 2*K-1 ) = Z( K ) 00227 Z( 2*K-2 ) = ZERO 00228 Z( 2*K-3 ) = Z( K-1 ) 00229 30 CONTINUE 00230 * 00231 I0 = 1 00232 N0 = N 00233 * 00234 * Reverse the qd-array, if warranted. 00235 * 00236 IF( CBIAS*Z( 4*I0-3 ).LT.Z( 4*N0-3 ) ) THEN 00237 IPN4 = 4*( I0+N0 ) 00238 DO 40 I4 = 4*I0, 2*( I0+N0-1 ), 4 00239 TEMP = Z( I4-3 ) 00240 Z( I4-3 ) = Z( IPN4-I4-3 ) 00241 Z( IPN4-I4-3 ) = TEMP 00242 TEMP = Z( I4-1 ) 00243 Z( I4-1 ) = Z( IPN4-I4-5 ) 00244 Z( IPN4-I4-5 ) = TEMP 00245 40 CONTINUE 00246 END IF 00247 * 00248 * Initial split checking via dqd and Li's test. 00249 * 00250 PP = 0 00251 * 00252 DO 80 K = 1, 2 00253 * 00254 D = Z( 4*N0+PP-3 ) 00255 DO 50 I4 = 4*( N0-1 ) + PP, 4*I0 + PP, -4 00256 IF( Z( I4-1 ).LE.TOL2*D ) THEN 00257 Z( I4-1 ) = -ZERO 00258 D = Z( I4-3 ) 00259 ELSE 00260 D = Z( I4-3 )*( D / ( D+Z( I4-1 ) ) ) 00261 END IF 00262 50 CONTINUE 00263 * 00264 * dqd maps Z to ZZ plus Li's test. 00265 * 00266 EMIN = Z( 4*I0+PP+1 ) 00267 D = Z( 4*I0+PP-3 ) 00268 DO 60 I4 = 4*I0 + PP, 4*( N0-1 ) + PP, 4 00269 Z( I4-2*PP-2 ) = D + Z( I4-1 ) 00270 IF( Z( I4-1 ).LE.TOL2*D ) THEN 00271 Z( I4-1 ) = -ZERO 00272 Z( I4-2*PP-2 ) = D 00273 Z( I4-2*PP ) = ZERO 00274 D = Z( I4+1 ) 00275 ELSE IF( SAFMIN*Z( I4+1 ).LT.Z( I4-2*PP-2 ) .AND. 00276 $ SAFMIN*Z( I4-2*PP-2 ).LT.Z( I4+1 ) ) THEN 00277 TEMP = Z( I4+1 ) / Z( I4-2*PP-2 ) 00278 Z( I4-2*PP ) = Z( I4-1 )*TEMP 00279 D = D*TEMP 00280 ELSE 00281 Z( I4-2*PP ) = Z( I4+1 )*( Z( I4-1 ) / Z( I4-2*PP-2 ) ) 00282 D = Z( I4+1 )*( D / Z( I4-2*PP-2 ) ) 00283 END IF 00284 EMIN = MIN( EMIN, Z( I4-2*PP ) ) 00285 60 CONTINUE 00286 Z( 4*N0-PP-2 ) = D 00287 * 00288 * Now find qmax. 00289 * 00290 QMAX = Z( 4*I0-PP-2 ) 00291 DO 70 I4 = 4*I0 - PP + 2, 4*N0 - PP - 2, 4 00292 QMAX = MAX( QMAX, Z( I4 ) ) 00293 70 CONTINUE 00294 * 00295 * Prepare for the next iteration on K. 00296 * 00297 PP = 1 - PP 00298 80 CONTINUE 00299 * 00300 * Initialise variables to pass to SLASQ3. 00301 * 00302 TTYPE = 0 00303 DMIN1 = ZERO 00304 DMIN2 = ZERO 00305 DN = ZERO 00306 DN1 = ZERO 00307 DN2 = ZERO 00308 G = ZERO 00309 TAU = ZERO 00310 * 00311 ITER = 2 00312 NFAIL = 0 00313 NDIV = 2*( N0-I0 ) 00314 * 00315 DO 160 IWHILA = 1, N + 1 00316 IF( N0.LT.1 ) 00317 $ GO TO 170 00318 * 00319 * While array unfinished do 00320 * 00321 * E(N0) holds the value of SIGMA when submatrix in I0:N0 00322 * splits from the rest of the array, but is negated. 00323 * 00324 DESIG = ZERO 00325 IF( N0.EQ.N ) THEN 00326 SIGMA = ZERO 00327 ELSE 00328 SIGMA = -Z( 4*N0-1 ) 00329 END IF 00330 IF( SIGMA.LT.ZERO ) THEN 00331 INFO = 1 00332 RETURN 00333 END IF 00334 * 00335 * Find last unreduced submatrix's top index I0, find QMAX and 00336 * EMIN. Find Gershgorin-type bound if Q's much greater than E's. 00337 * 00338 EMAX = ZERO 00339 IF( N0.GT.I0 ) THEN 00340 EMIN = ABS( Z( 4*N0-5 ) ) 00341 ELSE 00342 EMIN = ZERO 00343 END IF 00344 QMIN = Z( 4*N0-3 ) 00345 QMAX = QMIN 00346 DO 90 I4 = 4*N0, 8, -4 00347 IF( Z( I4-5 ).LE.ZERO ) 00348 $ GO TO 100 00349 IF( QMIN.GE.FOUR*EMAX ) THEN 00350 QMIN = MIN( QMIN, Z( I4-3 ) ) 00351 EMAX = MAX( EMAX, Z( I4-5 ) ) 00352 END IF 00353 QMAX = MAX( QMAX, Z( I4-7 )+Z( I4-5 ) ) 00354 EMIN = MIN( EMIN, Z( I4-5 ) ) 00355 90 CONTINUE 00356 I4 = 4 00357 * 00358 100 CONTINUE 00359 I0 = I4 / 4 00360 PP = 0 00361 * 00362 IF( N0-I0.GT.1 ) THEN 00363 DEE = Z( 4*I0-3 ) 00364 DEEMIN = DEE 00365 KMIN = I0 00366 DO 110 I4 = 4*I0+1, 4*N0-3, 4 00367 DEE = Z( I4 )*( DEE /( DEE+Z( I4-2 ) ) ) 00368 IF( DEE.LE.DEEMIN ) THEN 00369 DEEMIN = DEE 00370 KMIN = ( I4+3 )/4 00371 END IF 00372 110 CONTINUE 00373 IF( (KMIN-I0)*2.LT.N0-KMIN .AND. 00374 $ DEEMIN.LE.HALF*Z(4*N0-3) ) THEN 00375 IPN4 = 4*( I0+N0 ) 00376 PP = 2 00377 DO 120 I4 = 4*I0, 2*( I0+N0-1 ), 4 00378 TEMP = Z( I4-3 ) 00379 Z( I4-3 ) = Z( IPN4-I4-3 ) 00380 Z( IPN4-I4-3 ) = TEMP 00381 TEMP = Z( I4-2 ) 00382 Z( I4-2 ) = Z( IPN4-I4-2 ) 00383 Z( IPN4-I4-2 ) = TEMP 00384 TEMP = Z( I4-1 ) 00385 Z( I4-1 ) = Z( IPN4-I4-5 ) 00386 Z( IPN4-I4-5 ) = TEMP 00387 TEMP = Z( I4 ) 00388 Z( I4 ) = Z( IPN4-I4-4 ) 00389 Z( IPN4-I4-4 ) = TEMP 00390 120 CONTINUE 00391 END IF 00392 END IF 00393 * 00394 * Put -(initial shift) into DMIN. 00395 * 00396 DMIN = -MAX( ZERO, QMIN-TWO*SQRT( QMIN )*SQRT( EMAX ) ) 00397 * 00398 * Now I0:N0 is unreduced. 00399 * PP = 0 for ping, PP = 1 for pong. 00400 * PP = 2 indicates that flipping was applied to the Z array and 00401 * and that the tests for deflation upon entry in SLASQ3 00402 * should not be performed. 00403 * 00404 NBIG = 30*( N0-I0+1 ) 00405 DO 140 IWHILB = 1, NBIG 00406 IF( I0.GT.N0 ) 00407 $ GO TO 150 00408 * 00409 * While submatrix unfinished take a good dqds step. 00410 * 00411 CALL SLASQ3( I0, N0, Z, PP, DMIN, SIGMA, DESIG, QMAX, NFAIL, 00412 $ ITER, NDIV, IEEE, TTYPE, DMIN1, DMIN2, DN, DN1, 00413 $ DN2, G, TAU ) 00414 * 00415 PP = 1 - PP 00416 * 00417 * When EMIN is very small check for splits. 00418 * 00419 IF( PP.EQ.0 .AND. N0-I0.GE.3 ) THEN 00420 IF( Z( 4*N0 ).LE.TOL2*QMAX .OR. 00421 $ Z( 4*N0-1 ).LE.TOL2*SIGMA ) THEN 00422 SPLT = I0 - 1 00423 QMAX = Z( 4*I0-3 ) 00424 EMIN = Z( 4*I0-1 ) 00425 OLDEMN = Z( 4*I0 ) 00426 DO 130 I4 = 4*I0, 4*( N0-3 ), 4 00427 IF( Z( I4 ).LE.TOL2*Z( I4-3 ) .OR. 00428 $ Z( I4-1 ).LE.TOL2*SIGMA ) THEN 00429 Z( I4-1 ) = -SIGMA 00430 SPLT = I4 / 4 00431 QMAX = ZERO 00432 EMIN = Z( I4+3 ) 00433 OLDEMN = Z( I4+4 ) 00434 ELSE 00435 QMAX = MAX( QMAX, Z( I4+1 ) ) 00436 EMIN = MIN( EMIN, Z( I4-1 ) ) 00437 OLDEMN = MIN( OLDEMN, Z( I4 ) ) 00438 END IF 00439 130 CONTINUE 00440 Z( 4*N0-1 ) = EMIN 00441 Z( 4*N0 ) = OLDEMN 00442 I0 = SPLT + 1 00443 END IF 00444 END IF 00445 * 00446 140 CONTINUE 00447 * 00448 INFO = 2 00449 RETURN 00450 * 00451 * end IWHILB 00452 * 00453 150 CONTINUE 00454 * 00455 160 CONTINUE 00456 * 00457 INFO = 3 00458 RETURN 00459 * 00460 * end IWHILA 00461 * 00462 170 CONTINUE 00463 * 00464 * Move q's to the front. 00465 * 00466 DO 180 K = 2, N 00467 Z( K ) = Z( 4*K-3 ) 00468 180 CONTINUE 00469 * 00470 * Sort and compute sum of eigenvalues. 00471 * 00472 CALL SLASRT( 'D', N, Z, IINFO ) 00473 * 00474 E = ZERO 00475 DO 190 K = N, 1, -1 00476 E = E + Z( K ) 00477 190 CONTINUE 00478 * 00479 * Store trace, sum(eigenvalues) and information on performance. 00480 * 00481 Z( 2*N+1 ) = TRACE 00482 Z( 2*N+2 ) = E 00483 Z( 2*N+3 ) = REAL( ITER ) 00484 Z( 2*N+4 ) = REAL( NDIV ) / REAL( N**2 ) 00485 Z( 2*N+5 ) = HUNDRD*NFAIL / REAL( ITER ) 00486 RETURN 00487 * 00488 * End of SLASQ2 00489 * 00490 END
1.7.3 
