1*> \brief \b SLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method. 2* 3* =========== DOCUMENTATION =========== 4* 5* Online html documentation available at 6* http://www.netlib.org/lapack/explore-html/ 7* 8*> \htmlonly 9*> Download SLAED0 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed0.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed0.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed0.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, 22* WORK, IWORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ 26* .. 27* .. Array Arguments .. 28* INTEGER IWORK( * ) 29* REAL D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ), 30* $ WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SLAED0 computes all eigenvalues and corresponding eigenvectors of a 40*> symmetric tridiagonal matrix using the divide and conquer method. 41*> \endverbatim 42* 43* Arguments: 44* ========== 45* 46*> \param[in] ICOMPQ 47*> \verbatim 48*> ICOMPQ is INTEGER 49*> = 0: Compute eigenvalues only. 50*> = 1: Compute eigenvectors of original dense symmetric matrix 51*> also. On entry, Q contains the orthogonal matrix used 52*> to reduce the original matrix to tridiagonal form. 53*> = 2: Compute eigenvalues and eigenvectors of tridiagonal 54*> matrix. 55*> \endverbatim 56*> 57*> \param[in] QSIZ 58*> \verbatim 59*> QSIZ is INTEGER 60*> The dimension of the orthogonal matrix used to reduce 61*> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. 62*> \endverbatim 63*> 64*> \param[in] N 65*> \verbatim 66*> N is INTEGER 67*> The dimension of the symmetric tridiagonal matrix. N >= 0. 68*> \endverbatim 69*> 70*> \param[in,out] D 71*> \verbatim 72*> D is REAL array, dimension (N) 73*> On entry, the main diagonal of the tridiagonal matrix. 74*> On exit, its eigenvalues. 75*> \endverbatim 76*> 77*> \param[in] E 78*> \verbatim 79*> E is REAL array, dimension (N-1) 80*> The off-diagonal elements of the tridiagonal matrix. 81*> On exit, E has been destroyed. 82*> \endverbatim 83*> 84*> \param[in,out] Q 85*> \verbatim 86*> Q is REAL array, dimension (LDQ, N) 87*> On entry, Q must contain an N-by-N orthogonal matrix. 88*> If ICOMPQ = 0 Q is not referenced. 89*> If ICOMPQ = 1 On entry, Q is a subset of the columns of the 90*> orthogonal matrix used to reduce the full 91*> matrix to tridiagonal form corresponding to 92*> the subset of the full matrix which is being 93*> decomposed at this time. 94*> If ICOMPQ = 2 On entry, Q will be the identity matrix. 95*> On exit, Q contains the eigenvectors of the 96*> tridiagonal matrix. 97*> \endverbatim 98*> 99*> \param[in] LDQ 100*> \verbatim 101*> LDQ is INTEGER 102*> The leading dimension of the array Q. If eigenvectors are 103*> desired, then LDQ >= max(1,N). In any case, LDQ >= 1. 104*> \endverbatim 105*> 106*> \param[out] QSTORE 107*> \verbatim 108*> QSTORE is REAL array, dimension (LDQS, N) 109*> Referenced only when ICOMPQ = 1. Used to store parts of 110*> the eigenvector matrix when the updating matrix multiplies 111*> take place. 112*> \endverbatim 113*> 114*> \param[in] LDQS 115*> \verbatim 116*> LDQS is INTEGER 117*> The leading dimension of the array QSTORE. If ICOMPQ = 1, 118*> then LDQS >= max(1,N). In any case, LDQS >= 1. 119*> \endverbatim 120*> 121*> \param[out] WORK 122*> \verbatim 123*> WORK is REAL array, 124*> If ICOMPQ = 0 or 1, the dimension of WORK must be at least 125*> 1 + 3*N + 2*N*lg N + 3*N**2 126*> ( lg( N ) = smallest integer k 127*> such that 2^k >= N ) 128*> If ICOMPQ = 2, the dimension of WORK must be at least 129*> 4*N + N**2. 130*> \endverbatim 131*> 132*> \param[out] IWORK 133*> \verbatim 134*> IWORK is INTEGER array, 135*> If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 136*> 6 + 6*N + 5*N*lg N. 137*> ( lg( N ) = smallest integer k 138*> such that 2^k >= N ) 139*> If ICOMPQ = 2, the dimension of IWORK must be at least 140*> 3 + 5*N. 141*> \endverbatim 142*> 143*> \param[out] INFO 144*> \verbatim 145*> INFO is INTEGER 146*> = 0: successful exit. 147*> < 0: if INFO = -i, the i-th argument had an illegal value. 148*> > 0: The algorithm failed to compute an eigenvalue while 149*> working on the submatrix lying in rows and columns 150*> INFO/(N+1) through mod(INFO,N+1). 151*> \endverbatim 152* 153* Authors: 154* ======== 155* 156*> \author Univ. of Tennessee 157*> \author Univ. of California Berkeley 158*> \author Univ. of Colorado Denver 159*> \author NAG Ltd. 160* 161*> \date September 2012 162* 163*> \ingroup auxOTHERcomputational 164* 165*> \par Contributors: 166* ================== 167*> 168*> Jeff Rutter, Computer Science Division, University of California 169*> at Berkeley, USA 170* 171* ===================================================================== 172 SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, 173 $ WORK, IWORK, INFO ) 174* 175* -- LAPACK computational routine (version 3.4.2) -- 176* -- LAPACK is a software package provided by Univ. of Tennessee, -- 177* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 178* September 2012 179* 180* .. Scalar Arguments .. 181 INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ 182* .. 183* .. Array Arguments .. 184 INTEGER IWORK( * ) 185 REAL D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ), 186 $ WORK( * ) 187* .. 188* 189* ===================================================================== 190* 191* .. Parameters .. 192 REAL ZERO, ONE, TWO 193 PARAMETER ( ZERO = 0.E0, ONE = 1.E0, TWO = 2.E0 ) 194* .. 195* .. Local Scalars .. 196 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM, 197 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM, 198 $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1, 199 $ SPM2, SUBMAT, SUBPBS, TLVLS 200 REAL TEMP 201* .. 202* .. External Subroutines .. 203 EXTERNAL SCOPY, SGEMM, SLACPY, SLAED1, SLAED7, SSTEQR, 204 $ XERBLA 205* .. 206* .. External Functions .. 207 INTEGER ILAENV 208 EXTERNAL ILAENV 209* .. 210* .. Intrinsic Functions .. 211 INTRINSIC ABS, INT, LOG, MAX, REAL 212* .. 213* .. Executable Statements .. 214* 215* Test the input parameters. 216* 217 INFO = 0 218* 219 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN 220 INFO = -1 221 ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN 222 INFO = -2 223 ELSE IF( N.LT.0 ) THEN 224 INFO = -3 225 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 226 INFO = -7 227 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN 228 INFO = -9 229 END IF 230 IF( INFO.NE.0 ) THEN 231 CALL XERBLA( 'SLAED0', -INFO ) 232 RETURN 233 END IF 234* 235* Quick return if possible 236* 237 IF( N.EQ.0 ) 238 $ RETURN 239* 240 SMLSIZ = ILAENV( 9, 'SLAED0', ' ', 0, 0, 0, 0 ) 241* 242* Determine the size and placement of the submatrices, and save in 243* the leading elements of IWORK. 244* 245 IWORK( 1 ) = N 246 SUBPBS = 1 247 TLVLS = 0 248 10 CONTINUE 249 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN 250 DO 20 J = SUBPBS, 1, -1 251 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 252 IWORK( 2*J-1 ) = IWORK( J ) / 2 253 20 CONTINUE 254 TLVLS = TLVLS + 1 255 SUBPBS = 2*SUBPBS 256 GO TO 10 257 END IF 258 DO 30 J = 2, SUBPBS 259 IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 260 30 CONTINUE 261* 262* Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 263* using rank-1 modifications (cuts). 264* 265 SPM1 = SUBPBS - 1 266 DO 40 I = 1, SPM1 267 SUBMAT = IWORK( I ) + 1 268 SMM1 = SUBMAT - 1 269 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) 270 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) 271 40 CONTINUE 272* 273 INDXQ = 4*N + 3 274 IF( ICOMPQ.NE.2 ) THEN 275* 276* Set up workspaces for eigenvalues only/accumulate new vectors 277* routine 278* 279 TEMP = LOG( REAL( N ) ) / LOG( TWO ) 280 LGN = INT( TEMP ) 281 IF( 2**LGN.LT.N ) 282 $ LGN = LGN + 1 283 IF( 2**LGN.LT.N ) 284 $ LGN = LGN + 1 285 IPRMPT = INDXQ + N + 1 286 IPERM = IPRMPT + N*LGN 287 IQPTR = IPERM + N*LGN 288 IGIVPT = IQPTR + N + 2 289 IGIVCL = IGIVPT + N*LGN 290* 291 IGIVNM = 1 292 IQ = IGIVNM + 2*N*LGN 293 IWREM = IQ + N**2 + 1 294* 295* Initialize pointers 296* 297 DO 50 I = 0, SUBPBS 298 IWORK( IPRMPT+I ) = 1 299 IWORK( IGIVPT+I ) = 1 300 50 CONTINUE 301 IWORK( IQPTR ) = 1 302 END IF 303* 304* Solve each submatrix eigenproblem at the bottom of the divide and 305* conquer tree. 306* 307 CURR = 0 308 DO 70 I = 0, SPM1 309 IF( I.EQ.0 ) THEN 310 SUBMAT = 1 311 MATSIZ = IWORK( 1 ) 312 ELSE 313 SUBMAT = IWORK( I ) + 1 314 MATSIZ = IWORK( I+1 ) - IWORK( I ) 315 END IF 316 IF( ICOMPQ.EQ.2 ) THEN 317 CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 318 $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO ) 319 IF( INFO.NE.0 ) 320 $ GO TO 130 321 ELSE 322 CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 323 $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK, 324 $ INFO ) 325 IF( INFO.NE.0 ) 326 $ GO TO 130 327 IF( ICOMPQ.EQ.1 ) THEN 328 CALL SGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE, 329 $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+ 330 $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ), 331 $ LDQS ) 332 END IF 333 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2 334 CURR = CURR + 1 335 END IF 336 K = 1 337 DO 60 J = SUBMAT, IWORK( I+1 ) 338 IWORK( INDXQ+J ) = K 339 K = K + 1 340 60 CONTINUE 341 70 CONTINUE 342* 343* Successively merge eigensystems of adjacent submatrices 344* into eigensystem for the corresponding larger matrix. 345* 346* while ( SUBPBS > 1 ) 347* 348 CURLVL = 1 349 80 CONTINUE 350 IF( SUBPBS.GT.1 ) THEN 351 SPM2 = SUBPBS - 2 352 DO 90 I = 0, SPM2, 2 353 IF( I.EQ.0 ) THEN 354 SUBMAT = 1 355 MATSIZ = IWORK( 2 ) 356 MSD2 = IWORK( 1 ) 357 CURPRB = 0 358 ELSE 359 SUBMAT = IWORK( I ) + 1 360 MATSIZ = IWORK( I+2 ) - IWORK( I ) 361 MSD2 = MATSIZ / 2 362 CURPRB = CURPRB + 1 363 END IF 364* 365* Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) 366* into an eigensystem of size MATSIZ. 367* SLAED1 is used only for the full eigensystem of a tridiagonal 368* matrix. 369* SLAED7 handles the cases in which eigenvalues only or eigenvalues 370* and eigenvectors of a full symmetric matrix (which was reduced to 371* tridiagonal form) are desired. 372* 373 IF( ICOMPQ.EQ.2 ) THEN 374 CALL SLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ), 375 $ LDQ, IWORK( INDXQ+SUBMAT ), 376 $ E( SUBMAT+MSD2-1 ), MSD2, WORK, 377 $ IWORK( SUBPBS+1 ), INFO ) 378 ELSE 379 CALL SLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB, 380 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, 381 $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ), 382 $ MSD2, WORK( IQ ), IWORK( IQPTR ), 383 $ IWORK( IPRMPT ), IWORK( IPERM ), 384 $ IWORK( IGIVPT ), IWORK( IGIVCL ), 385 $ WORK( IGIVNM ), WORK( IWREM ), 386 $ IWORK( SUBPBS+1 ), INFO ) 387 END IF 388 IF( INFO.NE.0 ) 389 $ GO TO 130 390 IWORK( I / 2+1 ) = IWORK( I+2 ) 391 90 CONTINUE 392 SUBPBS = SUBPBS / 2 393 CURLVL = CURLVL + 1 394 GO TO 80 395 END IF 396* 397* end while 398* 399* Re-merge the eigenvalues/vectors which were deflated at the final 400* merge step. 401* 402 IF( ICOMPQ.EQ.1 ) THEN 403 DO 100 I = 1, N 404 J = IWORK( INDXQ+I ) 405 WORK( I ) = D( J ) 406 CALL SCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 ) 407 100 CONTINUE 408 CALL SCOPY( N, WORK, 1, D, 1 ) 409 ELSE IF( ICOMPQ.EQ.2 ) THEN 410 DO 110 I = 1, N 411 J = IWORK( INDXQ+I ) 412 WORK( I ) = D( J ) 413 CALL SCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 ) 414 110 CONTINUE 415 CALL SCOPY( N, WORK, 1, D, 1 ) 416 CALL SLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ ) 417 ELSE 418 DO 120 I = 1, N 419 J = IWORK( INDXQ+I ) 420 WORK( I ) = D( J ) 421 120 CONTINUE 422 CALL SCOPY( N, WORK, 1, D, 1 ) 423 END IF 424 GO TO 140 425* 426 130 CONTINUE 427 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 428* 429 140 CONTINUE 430 RETURN 431* 432* End of SLAED0 433* 434 END 435c $Id$ 436