LAPACK 3.12.1
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ dlaed0()

subroutine dlaed0 ( integer icompq,
integer qsiz,
integer n,
double precision, dimension( * ) d,
double precision, dimension( * ) e,
double precision, dimension( ldq, * ) q,
integer ldq,
double precision, dimension( ldqs, * ) qstore,
integer ldqs,
double precision, dimension( * ) work,
integer, dimension( * ) iwork,
integer info )

DLAED0 used by DSTEDC. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.

Download DLAED0 + dependencies [TGZ] [ZIP] [TXT]

Purpose:
!>
!> DLAED0 computes all eigenvalues and corresponding eigenvectors of a
!> symmetric tridiagonal matrix using the divide and conquer method.
!>
Parameters
[in]ICOMPQ
!>          ICOMPQ is INTEGER
!>          = 0:  Compute eigenvalues only.
!>          = 1:  Compute eigenvectors of original dense symmetric matrix
!>                also.  On entry, Q contains the orthogonal matrix used
!>                to reduce the original matrix to tridiagonal form.
!>          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
!>                matrix.
!>
[in]QSIZ
!>          QSIZ is INTEGER
!>         The dimension of the orthogonal matrix used to reduce
!>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
!>
[in]N
!>          N is INTEGER
!>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
!>
[in,out]D
!>          D is DOUBLE PRECISION array, dimension (N)
!>         On entry, the main diagonal of the tridiagonal matrix.
!>         On exit, its eigenvalues.
!>
[in]E
!>          E is DOUBLE PRECISION array, dimension (N-1)
!>         The off-diagonal elements of the tridiagonal matrix.
!>         On exit, E has been destroyed.
!>
[in,out]Q
!>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
!>         On entry, Q must contain an N-by-N orthogonal matrix.
!>         If ICOMPQ = 0    Q is not referenced.
!>         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
!>                          orthogonal matrix used to reduce the full
!>                          matrix to tridiagonal form corresponding to
!>                          the subset of the full matrix which is being
!>                          decomposed at this time.
!>         If ICOMPQ = 2    On entry, Q will be the identity matrix.
!>                          On exit, Q contains the eigenvectors of the
!>                          tridiagonal matrix.
!>
[in]LDQ
!>          LDQ is INTEGER
!>         The leading dimension of the array Q.  If eigenvectors are
!>         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
!>
[out]QSTORE
!>          QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
!>         Referenced only when ICOMPQ = 1.  Used to store parts of
!>         the eigenvector matrix when the updating matrix multiplies
!>         take place.
!>
[in]LDQS
!>          LDQS is INTEGER
!>         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
!>         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
!>
[out]WORK
!>          WORK is DOUBLE PRECISION array,
!>         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
!>                     1 + 3*N + 2*N*lg N + 3*N**2
!>                     ( lg( N ) = smallest integer k
!>                                 such that 2^k >= N )
!>         If ICOMPQ = 2, the dimension of WORK must be at least
!>                     4*N + N**2.
!>
[out]IWORK
!>          IWORK is INTEGER array,
!>         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
!>                        6 + 6*N + 5*N*lg N.
!>                        ( lg( N ) = smallest integer k
!>                                    such that 2^k >= N )
!>         If ICOMPQ = 2, the dimension of IWORK must be at least
!>                        3 + 5*N.
!>
[out]INFO
!>          INFO is INTEGER
!>          = 0:  successful exit.
!>          < 0:  if INFO = -i, the i-th argument had an illegal value.
!>          > 0:  The algorithm failed to compute an eigenvalue while
!>                working on the submatrix lying in rows and columns
!>                INFO/(N+1) through mod(INFO,N+1).
!>
Author
Univ. of Tennessee
Univ. of California Berkeley
Univ. of Colorado Denver
NAG Ltd.
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 168 of file dlaed0.f.

170*
171* -- LAPACK computational routine --
172* -- LAPACK is a software package provided by Univ. of Tennessee, --
173* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
174*
175* .. Scalar Arguments ..
176 INTEGER ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
177* ..
178* .. Array Arguments ..
179 INTEGER IWORK( * )
180 DOUBLE PRECISION D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
181 $ WORK( * )
182* ..
183*
184* =====================================================================
185*
186* .. Parameters ..
187 DOUBLE PRECISION ZERO, ONE, TWO
188 parameter( zero = 0.d0, one = 1.d0, two = 2.d0 )
189* ..
190* .. Local Scalars ..
191 INTEGER CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
192 $ IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
193 $ J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
194 $ SPM2, SUBMAT, SUBPBS, TLVLS
195 DOUBLE PRECISION TEMP
196* ..
197* .. External Subroutines ..
198 EXTERNAL dcopy, dgemm, dlacpy, dlaed1, dlaed7,
199 $ dsteqr,
200 $ xerbla
201* ..
202* .. External Functions ..
203 INTEGER ILAENV
204 EXTERNAL ilaenv
205* ..
206* .. Intrinsic Functions ..
207 INTRINSIC abs, dble, int, log, max
208* ..
209* .. Executable Statements ..
210*
211* Test the input parameters.
212*
213 info = 0
214*
215 IF( icompq.LT.0 .OR. icompq.GT.2 ) THEN
216 info = -1
217 ELSE IF( ( icompq.EQ.1 ) .AND. ( qsiz.LT.max( 0, n ) ) ) THEN
218 info = -2
219 ELSE IF( n.LT.0 ) THEN
220 info = -3
221 ELSE IF( ldq.LT.max( 1, n ) ) THEN
222 info = -7
223 ELSE IF( ldqs.LT.max( 1, n ) ) THEN
224 info = -9
225 END IF
226 IF( info.NE.0 ) THEN
227 CALL xerbla( 'DLAED0', -info )
228 RETURN
229 END IF
230*
231* Quick return if possible
232*
233 IF( n.EQ.0 )
234 $ RETURN
235*
236 smlsiz = ilaenv( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
237*
238* Determine the size and placement of the submatrices, and save in
239* the leading elements of IWORK.
240*
241 iwork( 1 ) = n
242 subpbs = 1
243 tlvls = 0
244 10 CONTINUE
245 IF( iwork( subpbs ).GT.smlsiz ) THEN
246 DO 20 j = subpbs, 1, -1
247 iwork( 2*j ) = ( iwork( j )+1 ) / 2
248 iwork( 2*j-1 ) = iwork( j ) / 2
249 20 CONTINUE
250 tlvls = tlvls + 1
251 subpbs = 2*subpbs
252 GO TO 10
253 END IF
254 DO 30 j = 2, subpbs
255 iwork( j ) = iwork( j ) + iwork( j-1 )
256 30 CONTINUE
257*
258* Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
259* using rank-1 modifications (cuts).
260*
261 spm1 = subpbs - 1
262 DO 40 i = 1, spm1
263 submat = iwork( i ) + 1
264 smm1 = submat - 1
265 d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
266 d( submat ) = d( submat ) - abs( e( smm1 ) )
267 40 CONTINUE
268*
269 indxq = 4*n + 3
270 IF( icompq.NE.2 ) THEN
271*
272* Set up workspaces for eigenvalues only/accumulate new vectors
273* routine
274*
275 temp = log( dble( n ) ) / log( two )
276 lgn = int( temp )
277 IF( 2**lgn.LT.n )
278 $ lgn = lgn + 1
279 IF( 2**lgn.LT.n )
280 $ lgn = lgn + 1
281 iprmpt = indxq + n + 1
282 iperm = iprmpt + n*lgn
283 iqptr = iperm + n*lgn
284 igivpt = iqptr + n + 2
285 igivcl = igivpt + n*lgn
286*
287 igivnm = 1
288 iq = igivnm + 2*n*lgn
289 iwrem = iq + n**2 + 1
290*
291* Initialize pointers
292*
293 DO 50 i = 0, subpbs
294 iwork( iprmpt+i ) = 1
295 iwork( igivpt+i ) = 1
296 50 CONTINUE
297 iwork( iqptr ) = 1
298 END IF
299*
300* Solve each submatrix eigenproblem at the bottom of the divide and
301* conquer tree.
302*
303 curr = 0
304 DO 70 i = 0, spm1
305 IF( i.EQ.0 ) THEN
306 submat = 1
307 matsiz = iwork( 1 )
308 ELSE
309 submat = iwork( i ) + 1
310 matsiz = iwork( i+1 ) - iwork( i )
311 END IF
312 IF( icompq.EQ.2 ) THEN
313 CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
314 $ q( submat, submat ), ldq, work, info )
315 IF( info.NE.0 )
316 $ GO TO 130
317 ELSE
318 CALL dsteqr( 'I', matsiz, d( submat ), e( submat ),
319 $ work( iq-1+iwork( iqptr+curr ) ), matsiz, work,
320 $ info )
321 IF( info.NE.0 )
322 $ GO TO 130
323 IF( icompq.EQ.1 ) THEN
324 CALL dgemm( 'N', 'N', qsiz, matsiz, matsiz, one,
325 $ q( 1, submat ), ldq, work( iq-1+iwork( iqptr+
326 $ curr ) ), matsiz, zero, qstore( 1, submat ),
327 $ ldqs )
328 END IF
329 iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
330 curr = curr + 1
331 END IF
332 k = 1
333 DO 60 j = submat, iwork( i+1 )
334 iwork( indxq+j ) = k
335 k = k + 1
336 60 CONTINUE
337 70 CONTINUE
338*
339* Successively merge eigensystems of adjacent submatrices
340* into eigensystem for the corresponding larger matrix.
341*
342* while ( SUBPBS > 1 )
343*
344 curlvl = 1
345 80 CONTINUE
346 IF( subpbs.GT.1 ) THEN
347 spm2 = subpbs - 2
348 DO 90 i = 0, spm2, 2
349 IF( i.EQ.0 ) THEN
350 submat = 1
351 matsiz = iwork( 2 )
352 msd2 = iwork( 1 )
353 curprb = 0
354 ELSE
355 submat = iwork( i ) + 1
356 matsiz = iwork( i+2 ) - iwork( i )
357 msd2 = matsiz / 2
358 curprb = curprb + 1
359 END IF
360*
361* Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
362* into an eigensystem of size MATSIZ.
363* DLAED1 is used only for the full eigensystem of a tridiagonal
364* matrix.
365* DLAED7 handles the cases in which eigenvalues only or eigenvalues
366* and eigenvectors of a full symmetric matrix (which was reduced to
367* tridiagonal form) are desired.
368*
369 IF( icompq.EQ.2 ) THEN
370 CALL dlaed1( matsiz, d( submat ), q( submat, submat ),
371 $ ldq, iwork( indxq+submat ),
372 $ e( submat+msd2-1 ), msd2, work,
373 $ iwork( subpbs+1 ), info )
374 ELSE
375 CALL dlaed7( icompq, matsiz, qsiz, tlvls, curlvl,
376 $ curprb,
377 $ d( submat ), qstore( 1, submat ), ldqs,
378 $ iwork( indxq+submat ), e( submat+msd2-1 ),
379 $ msd2, work( iq ), iwork( iqptr ),
380 $ iwork( iprmpt ), iwork( iperm ),
381 $ iwork( igivpt ), iwork( igivcl ),
382 $ work( igivnm ), work( iwrem ),
383 $ iwork( subpbs+1 ), info )
384 END IF
385 IF( info.NE.0 )
386 $ GO TO 130
387 iwork( i / 2+1 ) = iwork( i+2 )
388 90 CONTINUE
389 subpbs = subpbs / 2
390 curlvl = curlvl + 1
391 GO TO 80
392 END IF
393*
394* end while
395*
396* Re-merge the eigenvalues/vectors which were deflated at the final
397* merge step.
398*
399 IF( icompq.EQ.1 ) THEN
400 DO 100 i = 1, n
401 j = iwork( indxq+i )
402 work( i ) = d( j )
403 CALL dcopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
404 100 CONTINUE
405 CALL dcopy( n, work, 1, d, 1 )
406 ELSE IF( icompq.EQ.2 ) THEN
407 DO 110 i = 1, n
408 j = iwork( indxq+i )
409 work( i ) = d( j )
410 CALL dcopy( n, q( 1, j ), 1, work( n*i+1 ), 1 )
411 110 CONTINUE
412 CALL dcopy( n, work, 1, d, 1 )
413 CALL dlacpy( 'A', n, n, work( n+1 ), n, q, ldq )
414 ELSE
415 DO 120 i = 1, n
416 j = iwork( indxq+i )
417 work( i ) = d( j )
418 120 CONTINUE
419 CALL dcopy( n, work, 1, d, 1 )
420 END IF
421 GO TO 140
422*
423 130 CONTINUE
424 info = submat*( n+1 ) + submat + matsiz - 1
425*
426 140 CONTINUE
427 RETURN
428*
429* End of DLAED0
430*
xerbla
subroutine xerbla(srname, info)
Definition cblat2.f:3285
dcopy
subroutine dcopy(n, dx, incx, dy, incy)
DCOPY
Definition dcopy.f:82
dgemm
subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
DGEMM
Definition dgemm.f:188
ilaenv
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
Definition ilaenv.f:160
dlacpy
subroutine dlacpy(uplo, m, n, a, lda, b, ldb)
DLACPY copies all or part of one two-dimensional array to another.
Definition dlacpy.f:101
dlaed1
subroutine dlaed1(n, d, q, ldq, indxq, rho, cutpnt, work, iwork, info)
DLAED1 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition dlaed1.f:162
dlaed7
subroutine dlaed7(icompq, n, qsiz, tlvls, curlvl, curpbm, d, q, ldq, indxq, rho, cutpnt, qstore, qptr, prmptr, perm, givptr, givcol, givnum, work, iwork, info)
DLAED7 used by DSTEDC. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition dlaed7.f:259
dsteqr
subroutine dsteqr(compz, n, d, e, z, ldz, work, info)
DSTEQR
Definition dsteqr.f:129
Here is the call graph for this function:
Here is the caller graph for this function: