![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLAED0 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLAED0 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, 00022 * WORK, IWORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ 00026 * .. 00027 * .. Array Arguments .. 00028 * INTEGER IWORK( * ) 00029 * DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ), 00030 * $ WORK( * ) 00031 * .. 00032 * 00033 * 00034 *> \par Purpose: 00035 * ============= 00036 *> 00037 *> \verbatim 00038 *> 00039 *> DLAED0 computes all eigenvalues and corresponding eigenvectors of a 00040 *> symmetric tridiagonal matrix using the divide and conquer method. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in] ICOMPQ 00047 *> \verbatim 00048 *> ICOMPQ is INTEGER 00049 *> = 0: Compute eigenvalues only. 00050 *> = 1: Compute eigenvectors of original dense symmetric matrix 00051 *> also. On entry, Q contains the orthogonal matrix used 00052 *> to reduce the original matrix to tridiagonal form. 00053 *> = 2: Compute eigenvalues and eigenvectors of tridiagonal 00054 *> matrix. 00055 *> \endverbatim 00056 *> 00057 *> \param[in] QSIZ 00058 *> \verbatim 00059 *> QSIZ is INTEGER 00060 *> The dimension of the orthogonal matrix used to reduce 00061 *> the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. 00062 *> \endverbatim 00063 *> 00064 *> \param[in] N 00065 *> \verbatim 00066 *> N is INTEGER 00067 *> The dimension of the symmetric tridiagonal matrix. N >= 0. 00068 *> \endverbatim 00069 *> 00070 *> \param[in,out] D 00071 *> \verbatim 00072 *> D is DOUBLE PRECISION array, dimension (N) 00073 *> On entry, the main diagonal of the tridiagonal matrix. 00074 *> On exit, its eigenvalues. 00075 *> \endverbatim 00076 *> 00077 *> \param[in] E 00078 *> \verbatim 00079 *> E is DOUBLE PRECISION array, dimension (N-1) 00080 *> The off-diagonal elements of the tridiagonal matrix. 00081 *> On exit, E has been destroyed. 00082 *> \endverbatim 00083 *> 00084 *> \param[in,out] Q 00085 *> \verbatim 00086 *> Q is DOUBLE PRECISION array, dimension (LDQ, N) 00087 *> On entry, Q must contain an N-by-N orthogonal matrix. 00088 *> If ICOMPQ = 0 Q is not referenced. 00089 *> If ICOMPQ = 1 On entry, Q is a subset of the columns of the 00090 *> orthogonal matrix used to reduce the full 00091 *> matrix to tridiagonal form corresponding to 00092 *> the subset of the full matrix which is being 00093 *> decomposed at this time. 00094 *> If ICOMPQ = 2 On entry, Q will be the identity matrix. 00095 *> On exit, Q contains the eigenvectors of the 00096 *> tridiagonal matrix. 00097 *> \endverbatim 00098 *> 00099 *> \param[in] LDQ 00100 *> \verbatim 00101 *> LDQ is INTEGER 00102 *> The leading dimension of the array Q. If eigenvectors are 00103 *> desired, then LDQ >= max(1,N). In any case, LDQ >= 1. 00104 *> \endverbatim 00105 *> 00106 *> \param[out] QSTORE 00107 *> \verbatim 00108 *> QSTORE is DOUBLE PRECISION array, dimension (LDQS, N) 00109 *> Referenced only when ICOMPQ = 1. Used to store parts of 00110 *> the eigenvector matrix when the updating matrix multiplies 00111 *> take place. 00112 *> \endverbatim 00113 *> 00114 *> \param[in] LDQS 00115 *> \verbatim 00116 *> LDQS is INTEGER 00117 *> The leading dimension of the array QSTORE. If ICOMPQ = 1, 00118 *> then LDQS >= max(1,N). In any case, LDQS >= 1. 00119 *> \endverbatim 00120 *> 00121 *> \param[out] WORK 00122 *> \verbatim 00123 *> WORK is DOUBLE PRECISION array, 00124 *> If ICOMPQ = 0 or 1, the dimension of WORK must be at least 00125 *> 1 + 3*N + 2*N*lg N + 3*N**2 00126 *> ( lg( N ) = smallest integer k 00127 *> such that 2^k >= N ) 00128 *> If ICOMPQ = 2, the dimension of WORK must be at least 00129 *> 4*N + N**2. 00130 *> \endverbatim 00131 *> 00132 *> \param[out] IWORK 00133 *> \verbatim 00134 *> IWORK is INTEGER array, 00135 *> If ICOMPQ = 0 or 1, the dimension of IWORK must be at least 00136 *> 6 + 6*N + 5*N*lg N. 00137 *> ( lg( N ) = smallest integer k 00138 *> such that 2^k >= N ) 00139 *> If ICOMPQ = 2, the dimension of IWORK must be at least 00140 *> 3 + 5*N. 00141 *> \endverbatim 00142 *> 00143 *> \param[out] INFO 00144 *> \verbatim 00145 *> INFO is INTEGER 00146 *> = 0: successful exit. 00147 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00148 *> > 0: The algorithm failed to compute an eigenvalue while 00149 *> working on the submatrix lying in rows and columns 00150 *> INFO/(N+1) through mod(INFO,N+1). 00151 *> \endverbatim 00152 * 00153 * Authors: 00154 * ======== 00155 * 00156 *> \author Univ. of Tennessee 00157 *> \author Univ. of California Berkeley 00158 *> \author Univ. of Colorado Denver 00159 *> \author NAG Ltd. 00160 * 00161 *> \date November 2011 00162 * 00163 *> \ingroup auxOTHERcomputational 00164 * 00165 *> \par Contributors: 00166 * ================== 00167 *> 00168 *> Jeff Rutter, Computer Science Division, University of California 00169 *> at Berkeley, USA 00170 * 00171 * ===================================================================== 00172 SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, 00173 $ WORK, IWORK, INFO ) 00174 * 00175 * -- LAPACK computational routine (version 3.4.0) -- 00176 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00178 * November 2011 00179 * 00180 * .. Scalar Arguments .. 00181 INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ 00182 * .. 00183 * .. Array Arguments .. 00184 INTEGER IWORK( * ) 00185 DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ), 00186 $ WORK( * ) 00187 * .. 00188 * 00189 * ===================================================================== 00190 * 00191 * .. Parameters .. 00192 DOUBLE PRECISION ZERO, ONE, TWO 00193 PARAMETER ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 ) 00194 * .. 00195 * .. Local Scalars .. 00196 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM, 00197 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM, 00198 $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1, 00199 $ SPM2, SUBMAT, SUBPBS, TLVLS 00200 DOUBLE PRECISION TEMP 00201 * .. 00202 * .. External Subroutines .. 00203 EXTERNAL DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR, 00204 $ XERBLA 00205 * .. 00206 * .. External Functions .. 00207 INTEGER ILAENV 00208 EXTERNAL ILAENV 00209 * .. 00210 * .. Intrinsic Functions .. 00211 INTRINSIC ABS, DBLE, INT, LOG, MAX 00212 * .. 00213 * .. Executable Statements .. 00214 * 00215 * Test the input parameters. 00216 * 00217 INFO = 0 00218 * 00219 IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN 00220 INFO = -1 00221 ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN 00222 INFO = -2 00223 ELSE IF( N.LT.0 ) THEN 00224 INFO = -3 00225 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 00226 INFO = -7 00227 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN 00228 INFO = -9 00229 END IF 00230 IF( INFO.NE.0 ) THEN 00231 CALL XERBLA( 'DLAED0', -INFO ) 00232 RETURN 00233 END IF 00234 * 00235 * Quick return if possible 00236 * 00237 IF( N.EQ.0 ) 00238 $ RETURN 00239 * 00240 SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 ) 00241 * 00242 * Determine the size and placement of the submatrices, and save in 00243 * the leading elements of IWORK. 00244 * 00245 IWORK( 1 ) = N 00246 SUBPBS = 1 00247 TLVLS = 0 00248 10 CONTINUE 00249 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN 00250 DO 20 J = SUBPBS, 1, -1 00251 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 00252 IWORK( 2*J-1 ) = IWORK( J ) / 2 00253 20 CONTINUE 00254 TLVLS = TLVLS + 1 00255 SUBPBS = 2*SUBPBS 00256 GO TO 10 00257 END IF 00258 DO 30 J = 2, SUBPBS 00259 IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 00260 30 CONTINUE 00261 * 00262 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 00263 * using rank-1 modifications (cuts). 00264 * 00265 SPM1 = SUBPBS - 1 00266 DO 40 I = 1, SPM1 00267 SUBMAT = IWORK( I ) + 1 00268 SMM1 = SUBMAT - 1 00269 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) 00270 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) 00271 40 CONTINUE 00272 * 00273 INDXQ = 4*N + 3 00274 IF( ICOMPQ.NE.2 ) THEN 00275 * 00276 * Set up workspaces for eigenvalues only/accumulate new vectors 00277 * routine 00278 * 00279 TEMP = LOG( DBLE( N ) ) / LOG( TWO ) 00280 LGN = INT( TEMP ) 00281 IF( 2**LGN.LT.N ) 00282 $ LGN = LGN + 1 00283 IF( 2**LGN.LT.N ) 00284 $ LGN = LGN + 1 00285 IPRMPT = INDXQ + N + 1 00286 IPERM = IPRMPT + N*LGN 00287 IQPTR = IPERM + N*LGN 00288 IGIVPT = IQPTR + N + 2 00289 IGIVCL = IGIVPT + N*LGN 00290 * 00291 IGIVNM = 1 00292 IQ = IGIVNM + 2*N*LGN 00293 IWREM = IQ + N**2 + 1 00294 * 00295 * Initialize pointers 00296 * 00297 DO 50 I = 0, SUBPBS 00298 IWORK( IPRMPT+I ) = 1 00299 IWORK( IGIVPT+I ) = 1 00300 50 CONTINUE 00301 IWORK( IQPTR ) = 1 00302 END IF 00303 * 00304 * Solve each submatrix eigenproblem at the bottom of the divide and 00305 * conquer tree. 00306 * 00307 CURR = 0 00308 DO 70 I = 0, SPM1 00309 IF( I.EQ.0 ) THEN 00310 SUBMAT = 1 00311 MATSIZ = IWORK( 1 ) 00312 ELSE 00313 SUBMAT = IWORK( I ) + 1 00314 MATSIZ = IWORK( I+1 ) - IWORK( I ) 00315 END IF 00316 IF( ICOMPQ.EQ.2 ) THEN 00317 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 00318 $ Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO ) 00319 IF( INFO.NE.0 ) 00320 $ GO TO 130 00321 ELSE 00322 CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 00323 $ WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK, 00324 $ INFO ) 00325 IF( INFO.NE.0 ) 00326 $ GO TO 130 00327 IF( ICOMPQ.EQ.1 ) THEN 00328 CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE, 00329 $ Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+ 00330 $ CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ), 00331 $ LDQS ) 00332 END IF 00333 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2 00334 CURR = CURR + 1 00335 END IF 00336 K = 1 00337 DO 60 J = SUBMAT, IWORK( I+1 ) 00338 IWORK( INDXQ+J ) = K 00339 K = K + 1 00340 60 CONTINUE 00341 70 CONTINUE 00342 * 00343 * Successively merge eigensystems of adjacent submatrices 00344 * into eigensystem for the corresponding larger matrix. 00345 * 00346 * while ( SUBPBS > 1 ) 00347 * 00348 CURLVL = 1 00349 80 CONTINUE 00350 IF( SUBPBS.GT.1 ) THEN 00351 SPM2 = SUBPBS - 2 00352 DO 90 I = 0, SPM2, 2 00353 IF( I.EQ.0 ) THEN 00354 SUBMAT = 1 00355 MATSIZ = IWORK( 2 ) 00356 MSD2 = IWORK( 1 ) 00357 CURPRB = 0 00358 ELSE 00359 SUBMAT = IWORK( I ) + 1 00360 MATSIZ = IWORK( I+2 ) - IWORK( I ) 00361 MSD2 = MATSIZ / 2 00362 CURPRB = CURPRB + 1 00363 END IF 00364 * 00365 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) 00366 * into an eigensystem of size MATSIZ. 00367 * DLAED1 is used only for the full eigensystem of a tridiagonal 00368 * matrix. 00369 * DLAED7 handles the cases in which eigenvalues only or eigenvalues 00370 * and eigenvectors of a full symmetric matrix (which was reduced to 00371 * tridiagonal form) are desired. 00372 * 00373 IF( ICOMPQ.EQ.2 ) THEN 00374 CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ), 00375 $ LDQ, IWORK( INDXQ+SUBMAT ), 00376 $ E( SUBMAT+MSD2-1 ), MSD2, WORK, 00377 $ IWORK( SUBPBS+1 ), INFO ) 00378 ELSE 00379 CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB, 00380 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, 00381 $ IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ), 00382 $ MSD2, WORK( IQ ), IWORK( IQPTR ), 00383 $ IWORK( IPRMPT ), IWORK( IPERM ), 00384 $ IWORK( IGIVPT ), IWORK( IGIVCL ), 00385 $ WORK( IGIVNM ), WORK( IWREM ), 00386 $ IWORK( SUBPBS+1 ), INFO ) 00387 END IF 00388 IF( INFO.NE.0 ) 00389 $ GO TO 130 00390 IWORK( I / 2+1 ) = IWORK( I+2 ) 00391 90 CONTINUE 00392 SUBPBS = SUBPBS / 2 00393 CURLVL = CURLVL + 1 00394 GO TO 80 00395 END IF 00396 * 00397 * end while 00398 * 00399 * Re-merge the eigenvalues/vectors which were deflated at the final 00400 * merge step. 00401 * 00402 IF( ICOMPQ.EQ.1 ) THEN 00403 DO 100 I = 1, N 00404 J = IWORK( INDXQ+I ) 00405 WORK( I ) = D( J ) 00406 CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 ) 00407 100 CONTINUE 00408 CALL DCOPY( N, WORK, 1, D, 1 ) 00409 ELSE IF( ICOMPQ.EQ.2 ) THEN 00410 DO 110 I = 1, N 00411 J = IWORK( INDXQ+I ) 00412 WORK( I ) = D( J ) 00413 CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 ) 00414 110 CONTINUE 00415 CALL DCOPY( N, WORK, 1, D, 1 ) 00416 CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ ) 00417 ELSE 00418 DO 120 I = 1, N 00419 J = IWORK( INDXQ+I ) 00420 WORK( I ) = D( J ) 00421 120 CONTINUE 00422 CALL DCOPY( N, WORK, 1, D, 1 ) 00423 END IF 00424 GO TO 140 00425 * 00426 130 CONTINUE 00427 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 00428 * 00429 140 CONTINUE 00430 RETURN 00431 * 00432 * End of DLAED0 00433 * 00434 END