 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE CLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK, 00002 $ 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 INFO, LDQ, LDQS, N, QSIZ 00011 * .. 00012 * .. Array Arguments .. 00013 INTEGER IWORK( * ) 00014 REAL D( * ), E( * ), RWORK( * ) 00015 COMPLEX Q( LDQ, * ), QSTORE( LDQS, * ) 00016 * .. 00017 * 00018 * Purpose 00019 * ======= 00020 * 00021 * Using the divide and conquer method, CLAED0 computes all eigenvalues 00022 * of a symmetric tridiagonal matrix which is one diagonal block of 00023 * those from reducing a dense or band Hermitian matrix and 00024 * corresponding eigenvectors of the dense or band matrix. 00025 * 00026 * Arguments 00027 * ========= 00028 * 00029 * QSIZ (input) INTEGER 00030 * The dimension of the unitary matrix used to reduce 00031 * the full matrix to tridiagonal form. QSIZ >= N if ICOMPQ = 1. 00032 * 00033 * N (input) INTEGER 00034 * The dimension of the symmetric tridiagonal matrix. N >= 0. 00035 * 00036 * D (input/output) REAL array, dimension (N) 00037 * On entry, the diagonal elements of the tridiagonal matrix. 00038 * On exit, the eigenvalues in ascending order. 00039 * 00040 * E (input/output) REAL array, dimension (N-1) 00041 * On entry, the off-diagonal elements of the tridiagonal matrix. 00042 * On exit, E has been destroyed. 00043 * 00044 * Q (input/output) COMPLEX array, dimension (LDQ,N) 00045 * On entry, Q must contain an QSIZ x N matrix whose columns 00046 * unitarily orthonormal. It is a part of the unitary matrix 00047 * that reduces the full dense Hermitian matrix to a 00048 * (reducible) symmetric tridiagonal matrix. 00049 * 00050 * LDQ (input) INTEGER 00051 * The leading dimension of the array Q. LDQ >= max(1,N). 00052 * 00053 * IWORK (workspace) INTEGER array, 00054 * the dimension of IWORK must be at least 00055 * 6 + 6*N + 5*N*lg N 00056 * ( lg( N ) = smallest integer k 00057 * such that 2^k >= N ) 00058 * 00059 * RWORK (workspace) REAL array, 00060 * dimension (1 + 3*N + 2*N*lg N + 3*N**2) 00061 * ( lg( N ) = smallest integer k 00062 * such that 2^k >= N ) 00063 * 00064 * QSTORE (workspace) COMPLEX array, dimension (LDQS, N) 00065 * Used to store parts of 00066 * the eigenvector matrix when the updating matrix multiplies 00067 * take place. 00068 * 00069 * LDQS (input) INTEGER 00070 * The leading dimension of the array QSTORE. 00071 * LDQS >= max(1,N). 00072 * 00073 * INFO (output) INTEGER 00074 * = 0: successful exit. 00075 * < 0: if INFO = -i, the i-th argument had an illegal value. 00076 * > 0: The algorithm failed to compute an eigenvalue while 00077 * working on the submatrix lying in rows and columns 00078 * INFO/(N+1) through mod(INFO,N+1). 00079 * 00080 * ===================================================================== 00081 * 00082 * Warning: N could be as big as QSIZ! 00083 * 00084 * .. Parameters .. 00085 REAL TWO 00086 PARAMETER ( TWO = 2.E+0 ) 00087 * .. 00088 * .. Local Scalars .. 00089 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM, 00090 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM, 00091 $ J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1, 00092 $ SPM1, SPM2, SUBMAT, SUBPBS, TLVLS 00093 REAL TEMP 00094 * .. 00095 * .. External Subroutines .. 00096 EXTERNAL CCOPY, CLACRM, CLAED7, SCOPY, SSTEQR, XERBLA 00097 * .. 00098 * .. External Functions .. 00099 INTEGER ILAENV 00100 EXTERNAL ILAENV 00101 * .. 00102 * .. Intrinsic Functions .. 00103 INTRINSIC ABS, INT, LOG, MAX, REAL 00104 * .. 00105 * .. Executable Statements .. 00106 * 00107 * Test the input parameters. 00108 * 00109 INFO = 0 00110 * 00111 * IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN 00112 * INFO = -1 00113 * ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) ) 00114 * $ THEN 00115 IF( QSIZ.LT.MAX( 0, N ) ) THEN 00116 INFO = -1 00117 ELSE IF( N.LT.0 ) THEN 00118 INFO = -2 00119 ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN 00120 INFO = -6 00121 ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN 00122 INFO = -8 00123 END IF 00124 IF( INFO.NE.0 ) THEN 00125 CALL XERBLA( 'CLAED0', -INFO ) 00126 RETURN 00127 END IF 00128 * 00129 * Quick return if possible 00130 * 00131 IF( N.EQ.0 ) 00132 $ RETURN 00133 * 00134 SMLSIZ = ILAENV( 9, 'CLAED0', ' ', 0, 0, 0, 0 ) 00135 * 00136 * Determine the size and placement of the submatrices, and save in 00137 * the leading elements of IWORK. 00138 * 00139 IWORK( 1 ) = N 00140 SUBPBS = 1 00141 TLVLS = 0 00142 10 CONTINUE 00143 IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN 00144 DO 20 J = SUBPBS, 1, -1 00145 IWORK( 2*J ) = ( IWORK( J )+1 ) / 2 00146 IWORK( 2*J-1 ) = IWORK( J ) / 2 00147 20 CONTINUE 00148 TLVLS = TLVLS + 1 00149 SUBPBS = 2*SUBPBS 00150 GO TO 10 00151 END IF 00152 DO 30 J = 2, SUBPBS 00153 IWORK( J ) = IWORK( J ) + IWORK( J-1 ) 00154 30 CONTINUE 00155 * 00156 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1 00157 * using rank-1 modifications (cuts). 00158 * 00159 SPM1 = SUBPBS - 1 00160 DO 40 I = 1, SPM1 00161 SUBMAT = IWORK( I ) + 1 00162 SMM1 = SUBMAT - 1 00163 D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) ) 00164 D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) ) 00165 40 CONTINUE 00166 * 00167 INDXQ = 4*N + 3 00168 * 00169 * Set up workspaces for eigenvalues only/accumulate new vectors 00170 * routine 00171 * 00172 TEMP = LOG( REAL( N ) ) / LOG( TWO ) 00173 LGN = INT( TEMP ) 00174 IF( 2**LGN.LT.N ) 00175 $ LGN = LGN + 1 00176 IF( 2**LGN.LT.N ) 00177 $ LGN = LGN + 1 00178 IPRMPT = INDXQ + N + 1 00179 IPERM = IPRMPT + N*LGN 00180 IQPTR = IPERM + N*LGN 00181 IGIVPT = IQPTR + N + 2 00182 IGIVCL = IGIVPT + N*LGN 00183 * 00184 IGIVNM = 1 00185 IQ = IGIVNM + 2*N*LGN 00186 IWREM = IQ + N**2 + 1 00187 * Initialize pointers 00188 DO 50 I = 0, SUBPBS 00189 IWORK( IPRMPT+I ) = 1 00190 IWORK( IGIVPT+I ) = 1 00191 50 CONTINUE 00192 IWORK( IQPTR ) = 1 00193 * 00194 * Solve each submatrix eigenproblem at the bottom of the divide and 00195 * conquer tree. 00196 * 00197 CURR = 0 00198 DO 70 I = 0, SPM1 00199 IF( I.EQ.0 ) THEN 00200 SUBMAT = 1 00201 MATSIZ = IWORK( 1 ) 00202 ELSE 00203 SUBMAT = IWORK( I ) + 1 00204 MATSIZ = IWORK( I+1 ) - IWORK( I ) 00205 END IF 00206 LL = IQ - 1 + IWORK( IQPTR+CURR ) 00207 CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ), 00208 $ RWORK( LL ), MATSIZ, RWORK, INFO ) 00209 CALL CLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ), 00210 $ MATSIZ, QSTORE( 1, SUBMAT ), LDQS, 00211 $ RWORK( IWREM ) ) 00212 IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2 00213 CURR = CURR + 1 00214 IF( INFO.GT.0 ) THEN 00215 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 00216 RETURN 00217 END IF 00218 K = 1 00219 DO 60 J = SUBMAT, IWORK( I+1 ) 00220 IWORK( INDXQ+J ) = K 00221 K = K + 1 00222 60 CONTINUE 00223 70 CONTINUE 00224 * 00225 * Successively merge eigensystems of adjacent submatrices 00226 * into eigensystem for the corresponding larger matrix. 00227 * 00228 * while ( SUBPBS > 1 ) 00229 * 00230 CURLVL = 1 00231 80 CONTINUE 00232 IF( SUBPBS.GT.1 ) THEN 00233 SPM2 = SUBPBS - 2 00234 DO 90 I = 0, SPM2, 2 00235 IF( I.EQ.0 ) THEN 00236 SUBMAT = 1 00237 MATSIZ = IWORK( 2 ) 00238 MSD2 = IWORK( 1 ) 00239 CURPRB = 0 00240 ELSE 00241 SUBMAT = IWORK( I ) + 1 00242 MATSIZ = IWORK( I+2 ) - IWORK( I ) 00243 MSD2 = MATSIZ / 2 00244 CURPRB = CURPRB + 1 00245 END IF 00246 * 00247 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2) 00248 * into an eigensystem of size MATSIZ. CLAED7 handles the case 00249 * when the eigenvectors of a full or band Hermitian matrix (which 00250 * was reduced to tridiagonal form) are desired. 00251 * 00252 * I am free to use Q as a valuable working space until Loop 150. 00253 * 00254 CALL CLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB, 00255 $ D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS, 00256 $ E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ), 00257 $ RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ), 00258 $ IWORK( IPERM ), IWORK( IGIVPT ), 00259 $ IWORK( IGIVCL ), RWORK( IGIVNM ), 00260 $ Q( 1, SUBMAT ), RWORK( IWREM ), 00261 $ IWORK( SUBPBS+1 ), INFO ) 00262 IF( INFO.GT.0 ) THEN 00263 INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1 00264 RETURN 00265 END IF 00266 IWORK( I / 2+1 ) = IWORK( I+2 ) 00267 90 CONTINUE 00268 SUBPBS = SUBPBS / 2 00269 CURLVL = CURLVL + 1 00270 GO TO 80 00271 END IF 00272 * 00273 * end while 00274 * 00275 * Re-merge the eigenvalues/vectors which were deflated at the final 00276 * merge step. 00277 * 00278 DO 100 I = 1, N 00279 J = IWORK( INDXQ+I ) 00280 RWORK( I ) = D( J ) 00281 CALL CCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 ) 00282 100 CONTINUE 00283 CALL SCOPY( N, RWORK, 1, D, 1 ) 00284 * 00285 RETURN 00286 * 00287 * End of CLAED0 00288 * 00289 END
1.7.3 
