1 SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, 2 $ WORK, IWORK, INFO ) 3* 4* -- LAPACK routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ 11* .. 12* .. Array Arguments .. 13 INTEGER IWORK( * ) 14 DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ), 15 $ WORK( * ) 16* .. 17* 18* Purpose 19* ======= 20* 21* DLAED0 computes all eigenvalues and corresponding eigenvectors of a 22* symmetric tridiagonal matrix using the divide and conquer method. 23* 24* Arguments 25* ========= 26* 27* ICOMPQ (input) INTEGER 28* = 0: Compute eigenvalues only. 29* = 1: Compute eigenvectors of original dense symmetric matrix 30* also. On entry, Q contains the orthogonal matrix used 31* to reduce the original matrix to tridiagonal form. 32* = 2: Compute eigenvalues and eigenvectors of tridiagonal 33* matrix. 34* 35* QSIZ (input) INTEGER 36* The dimension of the orthogonal matrix used to reduce 37* the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. 38* 39* N (input) INTEGER 40* The dimension of the symmetric tridiagonal matrix. N >= 0. 41* 42* D (input/output) DOUBLE PRECISION array, dimension (N) 43* On entry, the main diagonal of the tridiagonal matrix. 44* On exit, its eigenvalues. 45* 46* E (input) DOUBLE PRECISION array, dimension (N-1) 47* The off-diagonal elements of the tridiagonal matrix. 48* On exit, E has been destroyed. 49* 50* Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) 51* On entry, Q must contain an N-by-N orthogonal matrix. 52* If ICOMPQ = 0 Q is not referenced. 53* If ICOMPQ = 1 On entry, Q is a subset of the columns of the 54* orthogonal matrix used to reduce the full 55* matrix to tridiagonal form corresponding to 56* the subset of the full matrix which is being 57* decomposed at this time. 58* If ICOMPQ = 2 On entry, Q will be the identity matrix. 59* On exit, Q contains the eigenvectors of the 60* tridiagonal matrix. 61* 62* LDQ (input) INTEGER 63* The leading dimension of the array Q. If eigenvectors are 64* desired, then LDQ >= max(1,N). In any case, LDQ >= 1. 65* 66* QSTORE (workspace) DOUBLE PRECISION array, dimension (LDQS, N) 67* Referenced only when ICOMPQ = 1. Used to store parts of 68* the eigenvector matrix when the updating matrix multiplies 69* take place. 70* 71* LDQS (input) INTEGER 72* The leading dimension of the array QSTORE. If ICOMPQ = 1, 73* then LDQS >= max(1,N). In any case, LDQS >= 1. 74* 75* WORK (workspace) DOUBLE PRECISION array, 76* If ICOMPQ = 0 or 1, the dimension of WORK must be at least 77* 1 + 3*N + 2*N*lg N + 2*N**2 78* ( lg( N ) = smallest integer k 79* such that 2^k >= N ) 80* If ICOMPQ = 2, the dimension of WORK must be at least 81* 4*N + N**2. 82* 83* IWORK (workspace) INTEGER array, 84* If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 85* 6 + 6*N + 5*N*lg N. 86* ( lg( N ) = smallest integer k 87* such that 2^k >= N ) 88* If ICOMPQ = 2, the dimension of IWORK must be at least 89* 3 + 5*N. 90* 91* INFO (output) INTEGER 92* = 0: successful exit. 93* < 0: if INFO = -i, the i-th argument had an illegal value. 94* > 0: The algorithm failed to compute an eigenvalue while 95* working on the submatrix lying in rows and columns 96* INFO/(N+1) through mod(INFO,N+1). 97* 98* Further Details 99* =============== 100* 101* Based on contributions by 102* Jeff Rutter, Computer Science Division, University of California 103* at Berkeley, USA 104* 105* ===================================================================== 106* 107* .. Parameters .. 108 DOUBLE PRECISION ZERO, ONE, TWO 109 PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 ) 110* .. 111* .. Local Scalars .. 112 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM, 113 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM, 114 $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1, 115 $ SPM2, SUBMAT, SUBPBS, TLVLS 116 DOUBLE PRECISION TEMP 117* .. 118* .. External Subroutines .. 119 EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR, 120 $ XERBLA 121* .. 122* .. External Functions .. 123 INTEGER ILAENV 124 EXTERNAL ILAENV 125* .. 126* .. Intrinsic Functions .. 127 INTRINSIC ABS, DBLE, INT, LOG, MAX 128* .. 129* .. Executable Statements .. 130* 131* Test the input parameters. 132* 133 INFO = 0 134* 135 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN 136 INFO = -1 137 ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN 138 INFO = -2 139 ELSE IF( N.LT.0 ) THEN 140 INFO = -3 141 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 142 INFO = -7 143 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN 144 INFO = -9 145 END IF 146 IF( INFO.NE.0 ) THEN 147 CALL XERBLA( 'DLAED0', -INFO ) 148 RETURN 149 END IF 150* 151* Quick return if possible 152* 153 IF( N.EQ.0 ) 154 $ RETURN 155* 156 SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 ) 157* 158* Determine the size and placement of the submatrices, and save in 159* the leading elements of IWORK. 160* 161 IWORK( 1 ) = N 162 SUBPBS = 1 163 TLVLS = 0 164 10 CONTINUE 165 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN 166 DO 20 J = SUBPBS, 1, -1 167 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 168 IWORK( 2*J-1 ) = IWORK( J ) / 2 169 20 CONTINUE 170 TLVLS = TLVLS + 1 171 SUBPBS = 2*SUBPBS 172 GO TO 10 173 END IF 174 DO 30 J = 2, SUBPBS 175 IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 176 30 CONTINUE 177* 178* Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 179* using rank-1 modifications (cuts). 180* 181 SPM1 = SUBPBS - 1 182 DO 40 I = 1, SPM1 183 SUBMAT = IWORK( I ) + 1 184 SMM1 = SUBMAT - 1 185 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) 186 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) 187 40 CONTINUE 188* 189 INDXQ = 4*N + 3 190 IF( ICOMPQ.NE.2 ) THEN 191* 192* Set up workspaces for eigenvalues only/accumulate new vectors 193* routine 194* 195 TEMP = LOG( DBLE( N ) ) / LOG( TWO ) 196 LGN = INT( TEMP ) 197 IF( 2**LGN.LT.N ) 198 $ LGN = LGN + 1 199 IF( 2**LGN.LT.N ) 200 $ LGN = LGN + 1 201 IPRMPT = INDXQ + N + 1 202 IPERM = IPRMPT + N*LGN 203 IQPTR = IPERM + N*LGN 204 IGIVPT = IQPTR + N + 2 205 IGIVCL = IGIVPT + N*LGN 206* 207 IGIVNM = 1 208 IQ = IGIVNM + 2*N*LGN 209 IWREM = IQ + N**2 + 1 210* 211* Initialize pointers 212* 213 DO 50 I = 0, SUBPBS 214 IWORK( IPRMPT+I ) = 1 215 IWORK( IGIVPT+I ) = 1 216 50 CONTINUE 217 IWORK( IQPTR ) = 1 218 END IF 219* 220* Solve each submatrix eigenproblem at the bottom of the divide and 221* conquer tree. 222* 223 CURR = 0 224 DO 70 I = 0, SPM1 225 IF( I.EQ.0 ) THEN 226 SUBMAT = 1 227 MATSIZ = IWORK( 1 ) 228 ELSE 229 SUBMAT = IWORK( I ) + 1 230 MATSIZ = IWORK( I+1 ) - IWORK( I ) 231 END IF 232 IF( ICOMPQ.EQ.2 ) THEN 233 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 234 $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO ) 235 IF( INFO.NE.0 ) 236 $ GO TO 130 237 ELSE 238 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 239 $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK, 240 $ INFO ) 241 IF( INFO.NE.0 ) 242 $ GO TO 130 243 IF( ICOMPQ.EQ.1 ) THEN 244 CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE, 245 $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+ 246 $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ), 247 $ LDQS ) 248 END IF 249 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2 250 CURR = CURR + 1 251 END IF 252 K = 1 253 DO 60 J = SUBMAT, IWORK( I+1 ) 254 IWORK( INDXQ+J ) = K 255 K = K + 1 256 60 CONTINUE 257 70 CONTINUE 258* 259* Successively merge eigensystems of adjacent submatrices 260* into eigensystem for the corresponding larger matrix. 261* 262* while ( SUBPBS > 1 ) 263* 264 CURLVL = 1 265 80 CONTINUE 266 IF( SUBPBS.GT.1 ) THEN 267 SPM2 = SUBPBS - 2 268 DO 90 I = 0, SPM2, 2 269 IF( I.EQ.0 ) THEN 270 SUBMAT = 1 271 MATSIZ = IWORK( 2 ) 272 MSD2 = IWORK( 1 ) 273 CURPRB = 0 274 ELSE 275 SUBMAT = IWORK( I ) + 1 276 MATSIZ = IWORK( I+2 ) - IWORK( I ) 277 MSD2 = MATSIZ / 2 278 CURPRB = CURPRB + 1 279 END IF 280* 281* Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) 282* into an eigensystem of size MATSIZ. 283* DLAED1 is used only for the full eigensystem of a tridiagonal 284* matrix. 285* DLAED7 handles the cases in which eigenvalues only or eigenvalues 286* and eigenvectors of a full symmetric matrix (which was reduced to 287* tridiagonal form) are desired. 288* 289 IF( ICOMPQ.EQ.2 ) THEN 290 CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ), 291 $ LDQ, IWORK( INDXQ+SUBMAT ), 292 $ E( SUBMAT+MSD2-1 ), MSD2, WORK, 293 $ IWORK( SUBPBS+1 ), INFO ) 294 ELSE 295 CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB, 296 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, 297 $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ), 298 $ MSD2, WORK( IQ ), IWORK( IQPTR ), 299 $ IWORK( IPRMPT ), IWORK( IPERM ), 300 $ IWORK( IGIVPT ), IWORK( IGIVCL ), 301 $ WORK( IGIVNM ), WORK( IWREM ), 302 $ IWORK( SUBPBS+1 ), INFO ) 303 END IF 304 IF( INFO.NE.0 ) 305 $ GO TO 130 306 IWORK( I / 2+1 ) = IWORK( I+2 ) 307 90 CONTINUE 308 SUBPBS = SUBPBS / 2 309 CURLVL = CURLVL + 1 310 GO TO 80 311 END IF 312* 313* end while 314* 315* Re-merge the eigenvalues/vectors which were deflated at the final 316* merge step. 317* 318 IF( ICOMPQ.EQ.1 ) THEN 319 DO 100 I = 1, N 320 J = IWORK( INDXQ+I ) 321 WORK( I ) = D( J ) 322 CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 ) 323 100 CONTINUE 324 CALL DCOPY( N, WORK, 1, D, 1 ) 325 ELSE IF( ICOMPQ.EQ.2 ) THEN 326 DO 110 I = 1, N 327 J = IWORK( INDXQ+I ) 328 WORK( I ) = D( J ) 329 CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 ) 330 110 CONTINUE 331 CALL DCOPY( N, WORK, 1, D, 1 ) 332 CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ ) 333 ELSE 334 DO 120 I = 1, N 335 J = IWORK( INDXQ+I ) 336 WORK( I ) = D( J ) 337 120 CONTINUE 338 CALL DCOPY( N, WORK, 1, D, 1 ) 339 END IF 340 GO TO 140 341* 342 130 CONTINUE 343 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 344* 345 140 CONTINUE 346 RETURN 347* 348* End of DLAED0 349* 350 END 351