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