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