LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
subroutine slaed0 ( integer  ICOMPQ,
integer  QSIZ,
integer  N,
real, dimension( * )  D,
real, dimension( * )  E,
real, dimension( ldq, * )  Q,
integer  LDQ,
real, dimension( ldqs, * )  QSTORE,
integer  LDQS,
real, dimension( * )  WORK,
integer, dimension( * )  IWORK,
integer  INFO 
)

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

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

Purpose: 
 SLAED0 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 REAL array, dimension (N)
         On entry, the main diagonal of the tridiagonal matrix.
         On exit, its eigenvalues.
[in]E
          E is REAL array, dimension (N-1)
         The off-diagonal elements of the tridiagonal matrix.
         On exit, E has been destroyed.
[in,out]Q
          Q is REAL 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 REAL 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 REAL 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.
Date
September 2012
Contributors:
Jeff Rutter, Computer Science Division, University of California at Berkeley, USA

Definition at line 174 of file slaed0.f.

174 *
175 * -- LAPACK computational routine (version 3.4.2) --
176 * -- LAPACK is a software package provided by Univ. of Tennessee, --
177 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
178 * September 2012
179 *
180 * .. Scalar Arguments ..
181  INTEGER icompq, info, ldq, ldqs, n, qsiz
182 * ..
183 * .. Array Arguments ..
184  INTEGER iwork( * )
185  REAL d( * ), e( * ), q( ldq, * ), qstore( ldqs, * ),
186  $ work( * )
187 * ..
188 *
189 * =====================================================================
190 *
191 * .. Parameters ..
192  REAL zero, one, two
193  parameter ( zero = 0.e0, one = 1.e0, two = 2.e0 )
194 * ..
195 * .. Local Scalars ..
196  INTEGER curlvl, curprb, curr, i, igivcl, igivnm,
197  $ igivpt, indxq, iperm, iprmpt, iq, iqptr, iwrem,
198  $ j, k, lgn, matsiz, msd2, smlsiz, smm1, spm1,
199  $ spm2, submat, subpbs, tlvls
200  REAL temp
201 * ..
202 * .. External Subroutines ..
203  EXTERNAL scopy, sgemm, slacpy, slaed1, slaed7, ssteqr,
204  $ xerbla
205 * ..
206 * .. External Functions ..
207  INTEGER ilaenv
208  EXTERNAL ilaenv
209 * ..
210 * .. Intrinsic Functions ..
211  INTRINSIC abs, int, log, max, real
212 * ..
213 * .. Executable Statements ..
214 *
215 * Test the input parameters.
216 *
217  info = 0
218 *
219  IF( icompq.LT.0 .OR. icompq.GT.2 ) THEN
220  info = -1
221  ELSE IF( ( icompq.EQ.1 ) .AND. ( qsiz.LT.max( 0, n ) ) ) THEN
222  info = -2
223  ELSE IF( n.LT.0 ) THEN
224  info = -3
225  ELSE IF( ldq.LT.max( 1, n ) ) THEN
226  info = -7
227  ELSE IF( ldqs.LT.max( 1, n ) ) THEN
228  info = -9
229  END IF
230  IF( info.NE.0 ) THEN
231  CALL xerbla( 'SLAED0', -info )
232  RETURN
233  END IF
234 *
235 * Quick return if possible
236 *
237  IF( n.EQ.0 )
238  $ RETURN
239 *
240  smlsiz = ilaenv( 9, 'SLAED0', ' ', 0, 0, 0, 0 )
241 *
242 * Determine the size and placement of the submatrices, and save in
243 * the leading elements of IWORK.
244 *
245  iwork( 1 ) = n
246  subpbs = 1
247  tlvls = 0
248  10 CONTINUE
249  IF( iwork( subpbs ).GT.smlsiz ) THEN
250  DO 20 j = subpbs, 1, -1
251  iwork( 2*j ) = ( iwork( j )+1 ) / 2
252  iwork( 2*j-1 ) = iwork( j ) / 2
253  20 CONTINUE
254  tlvls = tlvls + 1
255  subpbs = 2*subpbs
256  GO TO 10
257  END IF
258  DO 30 j = 2, subpbs
259  iwork( j ) = iwork( j ) + iwork( j-1 )
260  30 CONTINUE
261 *
262 * Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
263 * using rank-1 modifications (cuts).
264 *
265  spm1 = subpbs - 1
266  DO 40 i = 1, spm1
267  submat = iwork( i ) + 1
268  smm1 = submat - 1
269  d( smm1 ) = d( smm1 ) - abs( e( smm1 ) )
270  d( submat ) = d( submat ) - abs( e( smm1 ) )
271  40 CONTINUE
272 *
273  indxq = 4*n + 3
274  IF( icompq.NE.2 ) THEN
275 *
276 * Set up workspaces for eigenvalues only/accumulate new vectors
277 * routine
278 *
279  temp = log( REAL( N ) ) / log( two )
280  lgn = int( temp )
281  IF( 2**lgn.LT.n )
282  $ lgn = lgn + 1
283  IF( 2**lgn.LT.n )
284  $ lgn = lgn + 1
285  iprmpt = indxq + n + 1
286  iperm = iprmpt + n*lgn
287  iqptr = iperm + n*lgn
288  igivpt = iqptr + n + 2
289  igivcl = igivpt + n*lgn
290 *
291  igivnm = 1
292  iq = igivnm + 2*n*lgn
293  iwrem = iq + n**2 + 1
294 *
295 * Initialize pointers
296 *
297  DO 50 i = 0, subpbs
298  iwork( iprmpt+i ) = 1
299  iwork( igivpt+i ) = 1
300  50 CONTINUE
301  iwork( iqptr ) = 1
302  END IF
303 *
304 * Solve each submatrix eigenproblem at the bottom of the divide and
305 * conquer tree.
306 *
307  curr = 0
308  DO 70 i = 0, spm1
309  IF( i.EQ.0 ) THEN
310  submat = 1
311  matsiz = iwork( 1 )
312  ELSE
313  submat = iwork( i ) + 1
314  matsiz = iwork( i+1 ) - iwork( i )
315  END IF
316  IF( icompq.EQ.2 ) THEN
317  CALL ssteqr( 'I', matsiz, d( submat ), e( submat ),
318  $ q( submat, submat ), ldq, work, info )
319  IF( info.NE.0 )
320  $ GO TO 130
321  ELSE
322  CALL ssteqr( 'I', matsiz, d( submat ), e( submat ),
323  $ work( iq-1+iwork( iqptr+curr ) ), matsiz, work,
324  $ info )
325  IF( info.NE.0 )
326  $ GO TO 130
327  IF( icompq.EQ.1 ) THEN
328  CALL sgemm( 'N', 'N', qsiz, matsiz, matsiz, one,
329  $ q( 1, submat ), ldq, work( iq-1+iwork( iqptr+
330  $ curr ) ), matsiz, zero, qstore( 1, submat ),
331  $ ldqs )
332  END IF
333  iwork( iqptr+curr+1 ) = iwork( iqptr+curr ) + matsiz**2
334  curr = curr + 1
335  END IF
336  k = 1
337  DO 60 j = submat, iwork( i+1 )
338  iwork( indxq+j ) = k
339  k = k + 1
340  60 CONTINUE
341  70 CONTINUE
342 *
343 * Successively merge eigensystems of adjacent submatrices
344 * into eigensystem for the corresponding larger matrix.
345 *
346 * while ( SUBPBS > 1 )
347 *
348  curlvl = 1
349  80 CONTINUE
350  IF( subpbs.GT.1 ) THEN
351  spm2 = subpbs - 2
352  DO 90 i = 0, spm2, 2
353  IF( i.EQ.0 ) THEN
354  submat = 1
355  matsiz = iwork( 2 )
356  msd2 = iwork( 1 )
357  curprb = 0
358  ELSE
359  submat = iwork( i ) + 1
360  matsiz = iwork( i+2 ) - iwork( i )
361  msd2 = matsiz / 2
362  curprb = curprb + 1
363  END IF
364 *
365 * Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
366 * into an eigensystem of size MATSIZ.
367 * SLAED1 is used only for the full eigensystem of a tridiagonal
368 * matrix.
369 * SLAED7 handles the cases in which eigenvalues only or eigenvalues
370 * and eigenvectors of a full symmetric matrix (which was reduced to
371 * tridiagonal form) are desired.
372 *
373  IF( icompq.EQ.2 ) THEN
374  CALL slaed1( matsiz, d( submat ), q( submat, submat ),
375  $ ldq, iwork( indxq+submat ),
376  $ e( submat+msd2-1 ), msd2, work,
377  $ iwork( subpbs+1 ), info )
378  ELSE
379  CALL slaed7( icompq, matsiz, qsiz, tlvls, curlvl, curprb,
380  $ d( submat ), qstore( 1, submat ), ldqs,
381  $ iwork( indxq+submat ), e( submat+msd2-1 ),
382  $ msd2, work( iq ), iwork( iqptr ),
383  $ iwork( iprmpt ), iwork( iperm ),
384  $ iwork( igivpt ), iwork( igivcl ),
385  $ work( igivnm ), work( iwrem ),
386  $ iwork( subpbs+1 ), info )
387  END IF
388  IF( info.NE.0 )
389  $ GO TO 130
390  iwork( i / 2+1 ) = iwork( i+2 )
391  90 CONTINUE
392  subpbs = subpbs / 2
393  curlvl = curlvl + 1
394  GO TO 80
395  END IF
396 *
397 * end while
398 *
399 * Re-merge the eigenvalues/vectors which were deflated at the final
400 * merge step.
401 *
402  IF( icompq.EQ.1 ) THEN
403  DO 100 i = 1, n
404  j = iwork( indxq+i )
405  work( i ) = d( j )
406  CALL scopy( qsiz, qstore( 1, j ), 1, q( 1, i ), 1 )
407  100 CONTINUE
408  CALL scopy( n, work, 1, d, 1 )
409  ELSE IF( icompq.EQ.2 ) THEN
410  DO 110 i = 1, n
411  j = iwork( indxq+i )
412  work( i ) = d( j )
413  CALL scopy( n, q( 1, j ), 1, work( n*i+1 ), 1 )
414  110 CONTINUE
415  CALL scopy( n, work, 1, d, 1 )
416  CALL slacpy( 'A', n, n, work( n+1 ), n, q, ldq )
417  ELSE
418  DO 120 i = 1, n
419  j = iwork( indxq+i )
420  work( i ) = d( j )
421  120 CONTINUE
422  CALL scopy( n, work, 1, d, 1 )
423  END IF
424  GO TO 140
425 *
426  130 CONTINUE
427  info = submat*( n+1 ) + submat + matsiz - 1
428 *
429  140 CONTINUE
430  RETURN
431 *
432 * End of SLAED0
433 *
sgemm
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
Definition: sgemm.f:189
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
slacpy
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
slaed1
subroutine slaed1(N, D, Q, LDQ, INDXQ, RHO, CUTPNT, WORK, IWORK, INFO)
SLAED1 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: slaed1.f:165
ssteqr
subroutine ssteqr(COMPZ, N, D, E, Z, LDZ, WORK, INFO)
SSTEQR
Definition: ssteqr.f:133
ilaenv
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
Definition: tstiee.f:83
slaed7
subroutine slaed7(ICOMPQ, N, QSIZ, TLVLS, CURLVL, CURPBM, D, Q, LDQ, INDXQ, RHO, CUTPNT, QSTORE, QPTR, PRMPTR, PERM, GIVPTR, GIVCOL, GIVNUM, WORK, IWORK, INFO)
SLAED7 used by sstedc. Computes the updated eigensystem of a diagonal matrix after modification by a ...
Definition: slaed7.f:262
scopy
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
Definition: scopy.f:53

Here is the call graph for this function:

Here is the caller graph for this function: