SRC/dstedc.f

   1 *> \brief \b DSTEDC
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download DSTEDC + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dstedc.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dstedc.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dstedc.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
  22 *                          LIWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       CHARACTER          COMPZ
  26 *       INTEGER            INFO, LDZ, LIWORK, LWORK, N
  27 *       ..
  28 *       .. Array Arguments ..
  29 *       INTEGER            IWORK( * )
  30 *       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> DSTEDC computes all eigenvalues and, optionally, eigenvectors of a
  40 *> symmetric tridiagonal matrix using the divide and conquer method.
  41 *> The eigenvectors of a full or band real symmetric matrix can also be
  42 *> found if DSYTRD or DSPTRD or DSBTRD has been used to reduce this
  43 *> matrix to tridiagonal form.
  44 *>
  45 *> This code makes very mild assumptions about floating point
  46 *> arithmetic. It will work on machines with a guard digit in
  47 *> add/subtract, or on those binary machines without guard digits
  48 *> which subtract like the Cray X-MP, Cray Y-MP, Cray C-90, or Cray-2.
  49 *> It could conceivably fail on hexadecimal or decimal machines
  50 *> without guard digits, but we know of none.  See DLAED3 for details.
  51 *> \endverbatim
  52 *
  53 *  Arguments:
  54 *  ==========
  55 *
  56 *> \param[in] COMPZ
  57 *> \verbatim
  58 *>          COMPZ is CHARACTER*1
  59 *>          = 'N':  Compute eigenvalues only.
  60 *>          = 'I':  Compute eigenvectors of tridiagonal matrix also.
  61 *>          = 'V':  Compute eigenvectors of original dense symmetric
  62 *>                  matrix also.  On entry, Z contains the orthogonal
  63 *>                  matrix used to reduce the original matrix to
  64 *>                  tridiagonal form.
  65 *> \endverbatim
  66 *>
  67 *> \param[in] N
  68 *> \verbatim
  69 *>          N is INTEGER
  70 *>          The dimension of the symmetric tridiagonal matrix.  N >= 0.
  71 *> \endverbatim
  72 *>
  73 *> \param[in,out] D
  74 *> \verbatim
  75 *>          D is DOUBLE PRECISION array, dimension (N)
  76 *>          On entry, the diagonal elements of the tridiagonal matrix.
  77 *>          On exit, if INFO = 0, the eigenvalues in ascending order.
  78 *> \endverbatim
  79 *>
  80 *> \param[in,out] E
  81 *> \verbatim
  82 *>          E is DOUBLE PRECISION array, dimension (N-1)
  83 *>          On entry, the subdiagonal elements of the tridiagonal matrix.
  84 *>          On exit, E has been destroyed.
  85 *> \endverbatim
  86 *>
  87 *> \param[in,out] Z
  88 *> \verbatim
  89 *>          Z is DOUBLE PRECISION array, dimension (LDZ,N)
  90 *>          On entry, if COMPZ = 'V', then Z contains the orthogonal
  91 *>          matrix used in the reduction to tridiagonal form.
  92 *>          On exit, if INFO = 0, then if COMPZ = 'V', Z contains the
  93 *>          orthonormal eigenvectors of the original symmetric matrix,
  94 *>          and if COMPZ = 'I', Z contains the orthonormal eigenvectors
  95 *>          of the symmetric tridiagonal matrix.
  96 *>          If  COMPZ = 'N', then Z is not referenced.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDZ
 100 *> \verbatim
 101 *>          LDZ is INTEGER
 102 *>          The leading dimension of the array Z.  LDZ >= 1.
 103 *>          If eigenvectors are desired, then LDZ >= max(1,N).
 104 *> \endverbatim
 105 *>
 106 *> \param[out] WORK
 107 *> \verbatim
 108 *>          WORK is DOUBLE PRECISION array,
 109 *>                                         dimension (LWORK)
 110 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
 111 *> \endverbatim
 112 *>
 113 *> \param[in] LWORK
 114 *> \verbatim
 115 *>          LWORK is INTEGER
 116 *>          The dimension of the array WORK.
 117 *>          If COMPZ = 'N' or N <= 1 then LWORK must be at least 1.
 118 *>          If COMPZ = 'V' and N > 1 then LWORK must be at least
 119 *>                         ( 1 + 3*N + 2*N*lg N + 4*N**2 ),
 120 *>                         where lg( N ) = smallest integer k such
 121 *>                         that 2**k >= N.
 122 *>          If COMPZ = 'I' and N > 1 then LWORK must be at least
 123 *>                         ( 1 + 4*N + N**2 ).
 124 *>          Note that for COMPZ = 'I' or 'V', then if N is less than or
 125 *>          equal to the minimum divide size, usually 25, then LWORK need
 126 *>          only be max(1,2*(N-1)).
 127 *>
 128 *>          If LWORK = -1, then a workspace query is assumed; the routine
 129 *>          only calculates the optimal size of the WORK array, returns
 130 *>          this value as the first entry of the WORK array, and no error
 131 *>          message related to LWORK is issued by XERBLA.
 132 *> \endverbatim
 133 *>
 134 *> \param[out] IWORK
 135 *> \verbatim
 136 *>          IWORK is INTEGER array, dimension (MAX(1,LIWORK))
 137 *>          On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
 138 *> \endverbatim
 139 *>
 140 *> \param[in] LIWORK
 141 *> \verbatim
 142 *>          LIWORK is INTEGER
 143 *>          The dimension of the array IWORK.
 144 *>          If COMPZ = 'N' or N <= 1 then LIWORK must be at least 1.
 145 *>          If COMPZ = 'V' and N > 1 then LIWORK must be at least
 146 *>                         ( 6 + 6*N + 5*N*lg N ).
 147 *>          If COMPZ = 'I' and N > 1 then LIWORK must be at least
 148 *>                         ( 3 + 5*N ).
 149 *>          Note that for COMPZ = 'I' or 'V', then if N is less than or
 150 *>          equal to the minimum divide size, usually 25, then LIWORK
 151 *>          need only be 1.
 152 *>
 153 *>          If LIWORK = -1, then a workspace query is assumed; the
 154 *>          routine only calculates the optimal size of the IWORK array,
 155 *>          returns this value as the first entry of the IWORK array, and
 156 *>          no error message related to LIWORK is issued by XERBLA.
 157 *> \endverbatim
 158 *>
 159 *> \param[out] INFO
 160 *> \verbatim
 161 *>          INFO is INTEGER
 162 *>          = 0:  successful exit.
 163 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 164 *>          > 0:  The algorithm failed to compute an eigenvalue while
 165 *>                working on the submatrix lying in rows and columns
 166 *>                INFO/(N+1) through mod(INFO,N+1).
 167 *> \endverbatim
 168 *
 169 *  Authors:
 170 *  ========
 171 *
 172 *> \author Univ. of Tennessee
 173 *> \author Univ. of California Berkeley
 174 *> \author Univ. of Colorado Denver
 175 *> \author NAG Ltd.
 176 *
 177 *> \date November 2015
 178 *
 179 *> \ingroup auxOTHERcomputational
 180 *
 181 *> \par Contributors:
 182 *  ==================
 183 *>
 184 *> Jeff Rutter, Computer Science Division, University of California
 185 *> at Berkeley, USA \n
 186 *>  Modified by Francoise Tisseur, University of Tennessee
 187 *>
 188 *  =====================================================================
 189       SUBROUTINE DSTEDC( COMPZ, N, D, E, Z, LDZ, WORK, LWORK, IWORK,
 190      $                   LIWORK, INFO )
 191 *
 192 *  -- LAPACK computational routine (version 3.6.0) --
 193 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 194 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 195 *     November 2015
 196 *
 197 *     .. Scalar Arguments ..
 198       CHARACTER          COMPZ
 199       INTEGER            INFO, LDZ, LIWORK, LWORK, N
 200 *     ..
 201 *     .. Array Arguments ..
 202       INTEGER            IWORK( * )
 203       DOUBLE PRECISION   D( * ), E( * ), WORK( * ), Z( LDZ, * )
 204 *     ..
 205 *
 206 *  =====================================================================
 207 *
 208 *     .. Parameters ..
 209       DOUBLE PRECISION   ZERO, ONE, TWO
 210       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
 211 *     ..
 212 *     .. Local Scalars ..
 213       LOGICAL            LQUERY
 214       INTEGER            FINISH, I, ICOMPZ, II, J, K, LGN, LIWMIN,
 215      $                   LWMIN, M, SMLSIZ, START, STOREZ, STRTRW
 216       DOUBLE PRECISION   EPS, ORGNRM, P, TINY
 217 *     ..
 218 *     .. External Functions ..
 219       LOGICAL            LSAME
 220       INTEGER            ILAENV
 221       DOUBLE PRECISION   DLAMCH, DLANST
 222       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANST
 223 *     ..
 224 *     .. External Subroutines ..
 225       EXTERNAL           DGEMM, DLACPY, DLAED0, DLASCL, DLASET, DLASRT,
 226      $                   DSTEQR, DSTERF, DSWAP, XERBLA
 227 *     ..
 228 *     .. Intrinsic Functions ..
 229       INTRINSIC          ABS, DBLE, INT, LOG, MAX, MOD, SQRT
 230 *     ..
 231 *     .. Executable Statements ..
 232 *
 233 *     Test the input parameters.
 234 *
 235       INFO = 0
 236       LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
 237 *
 238       IF( LSAME( COMPZ, 'N' ) ) THEN
 239          ICOMPZ = 0
 240       ELSE IF( LSAME( COMPZ, 'V' ) ) THEN
 241          ICOMPZ = 1
 242       ELSE IF( LSAME( COMPZ, 'I' ) ) THEN
 243          ICOMPZ = 2
 244       ELSE
 245          ICOMPZ = -1
 246       END IF
 247       IF( ICOMPZ.LT.0 ) THEN
 248          INFO = -1
 249       ELSE IF( N.LT.0 ) THEN
 250          INFO = -2
 251       ELSE IF( ( LDZ.LT.1 ) .OR.
 252      $         ( ICOMPZ.GT.0 .AND. LDZ.LT.MAX( 1, N ) ) ) THEN
 253          INFO = -6
 254       END IF
 255 *
 256       IF( INFO.EQ.0 ) THEN
 257 *
 258 *        Compute the workspace requirements
 259 *
 260          SMLSIZ = ILAENV( 9, 'DSTEDC', ' ', 0, 0, 0, 0 )
 261          IF( N.LE.1 .OR. ICOMPZ.EQ.0 ) THEN
 262             LIWMIN = 1
 263             LWMIN = 1
 264          ELSE IF( N.LE.SMLSIZ ) THEN
 265             LIWMIN = 1
 266             LWMIN = 2*( N - 1 )
 267          ELSE
 268             LGN = INT( LOG( DBLE( N ) )/LOG( TWO ) )
 269             IF( 2**LGN.LT.N )
 270      $         LGN = LGN + 1
 271             IF( 2**LGN.LT.N )
 272      $         LGN = LGN + 1
 273             IF( ICOMPZ.EQ.1 ) THEN
 274                LWMIN = 1 + 3*N + 2*N*LGN + 4*N**2
 275                LIWMIN = 6 + 6*N + 5*N*LGN
 276             ELSE IF( ICOMPZ.EQ.2 ) THEN
 277                LWMIN = 1 + 4*N + N**2
 278                LIWMIN = 3 + 5*N
 279             END IF
 280          END IF
 281          WORK( 1 ) = LWMIN
 282          IWORK( 1 ) = LIWMIN
 283 *
 284          IF( LWORK.LT.LWMIN .AND. .NOT. LQUERY ) THEN
 285             INFO = -8
 286          ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT. LQUERY ) THEN
 287             INFO = -10
 288          END IF
 289       END IF
 290 *
 291       IF( INFO.NE.0 ) THEN
 292          CALL XERBLA( 'DSTEDC', -INFO )
 293          RETURN
 294       ELSE IF (LQUERY) THEN
 295          RETURN
 296       END IF
 297 *
 298 *     Quick return if possible
 299 *
 300       IF( N.EQ.0 )
 301      $   RETURN
 302       IF( N.EQ.1 ) THEN
 303          IF( ICOMPZ.NE.0 )
 304      $      Z( 1, 1 ) = ONE
 305          RETURN
 306       END IF
 307 *
 308 *     If the following conditional clause is removed, then the routine
 309 *     will use the Divide and Conquer routine to compute only the
 310 *     eigenvalues, which requires (3N + 3N**2) real workspace and
 311 *     (2 + 5N + 2N lg(N)) integer workspace.
 312 *     Since on many architectures DSTERF is much faster than any other
 313 *     algorithm for finding eigenvalues only, it is used here
 314 *     as the default. If the conditional clause is removed, then
 315 *     information on the size of workspace needs to be changed.
 316 *
 317 *     If COMPZ = 'N', use DSTERF to compute the eigenvalues.
 318 *
 319       IF( ICOMPZ.EQ.0 ) THEN
 320          CALL DSTERF( N, D, E, INFO )
 321          GO TO 50
 322       END IF
 323 *
 324 *     If N is smaller than the minimum divide size (SMLSIZ+1), then
 325 *     solve the problem with another solver.
 326 *
 327       IF( N.LE.SMLSIZ ) THEN
 328 *
 329          CALL DSTEQR( COMPZ, N, D, E, Z, LDZ, WORK, INFO )
 330 *
 331       ELSE
 332 *
 333 *        If COMPZ = 'V', the Z matrix must be stored elsewhere for later
 334 *        use.
 335 *
 336          IF( ICOMPZ.EQ.1 ) THEN
 337             STOREZ = 1 + N*N
 338          ELSE
 339             STOREZ = 1
 340          END IF
 341 *
 342          IF( ICOMPZ.EQ.2 ) THEN
 343             CALL DLASET( 'Full', N, N, ZERO, ONE, Z, LDZ )
 344          END IF
 345 *
 346 *        Scale.
 347 *
 348          ORGNRM = DLANST( 'M', N, D, E )
 349          IF( ORGNRM.EQ.ZERO )
 350      $      GO TO 50
 351 *
 352          EPS = DLAMCH( 'Epsilon' )
 353 *
 354          START = 1
 355 *
 356 *        while ( START <= N )
 357 *
 358    10    CONTINUE
 359          IF( START.LE.N ) THEN
 360 *
 361 *           Let FINISH be the position of the next subdiagonal entry
 362 *           such that E( FINISH ) <= TINY or FINISH = N if no such
 363 *           subdiagonal exists.  The matrix identified by the elements
 364 *           between START and FINISH constitutes an independent
 365 *           sub-problem.
 366 *
 367             FINISH = START
 368    20       CONTINUE
 369             IF( FINISH.LT.N ) THEN
 370                TINY = EPS*SQRT( ABS( D( FINISH ) ) )*
 371      $                    SQRT( ABS( D( FINISH+1 ) ) )
 372                IF( ABS( E( FINISH ) ).GT.TINY ) THEN
 373                   FINISH = FINISH + 1
 374                   GO TO 20
 375                END IF
 376             END IF
 377 *
 378 *           (Sub) Problem determined.  Compute its size and solve it.
 379 *
 380             M = FINISH - START + 1
 381             IF( M.EQ.1 ) THEN
 382                START = FINISH + 1
 383                GO TO 10
 384             END IF
 385             IF( M.GT.SMLSIZ ) THEN
 386 *
 387 *              Scale.
 388 *
 389                ORGNRM = DLANST( 'M', M, D( START ), E( START ) )
 390                CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M, 1, D( START ), M,
 391      $                      INFO )
 392                CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, M-1, 1, E( START ),
 393      $                      M-1, INFO )
 394 *
 395                IF( ICOMPZ.EQ.1 ) THEN
 396                   STRTRW = 1
 397                ELSE
 398                   STRTRW = START
 399                END IF
 400                CALL DLAED0( ICOMPZ, N, M, D( START ), E( START ),
 401      $                      Z( STRTRW, START ), LDZ, WORK( 1 ), N,
 402      $                      WORK( STOREZ ), IWORK, INFO )
 403                IF( INFO.NE.0 ) THEN
 404                   INFO = ( INFO / ( M+1 )+START-1 )*( N+1 ) +
 405      $                   MOD( INFO, ( M+1 ) ) + START - 1
 406                   GO TO 50
 407                END IF
 408 *
 409 *              Scale back.
 410 *
 411                CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, M, 1, D( START ), M,
 412      $                      INFO )
 413 *
 414             ELSE
 415                IF( ICOMPZ.EQ.1 ) THEN
 416 *
 417 *                 Since QR won't update a Z matrix which is larger than
 418 *                 the length of D, we must solve the sub-problem in a
 419 *                 workspace and then multiply back into Z.
 420 *
 421                   CALL DSTEQR( 'I', M, D( START ), E( START ), WORK, M,
 422      $                         WORK( M*M+1 ), INFO )
 423                   CALL DLACPY( 'A', N, M, Z( 1, START ), LDZ,
 424      $                         WORK( STOREZ ), N )
 425                   CALL DGEMM( 'N', 'N', N, M, M, ONE,
 426      $                        WORK( STOREZ ), N, WORK, M, ZERO,
 427      $                        Z( 1, START ), LDZ )
 428                ELSE IF( ICOMPZ.EQ.2 ) THEN
 429                   CALL DSTEQR( 'I', M, D( START ), E( START ),
 430      $                         Z( START, START ), LDZ, WORK, INFO )
 431                ELSE
 432                   CALL DSTERF( M, D( START ), E( START ), INFO )
 433                END IF
 434                IF( INFO.NE.0 ) THEN
 435                   INFO = START*( N+1 ) + FINISH
 436                   GO TO 50
 437                END IF
 438             END IF
 439 *
 440             START = FINISH + 1
 441             GO TO 10
 442          END IF
 443 *
 444 *        endwhile
 445 *
 446          IF( ICOMPZ.EQ.0 ) THEN
 447 *
 448 *          Use Quick Sort
 449 *
 450            CALL DLASRT( 'I', N, D, INFO )
 451 *
 452          ELSE
 453 *
 454 *          Use Selection Sort to minimize swaps of eigenvectors
 455 *
 456            DO 40 II = 2, N
 457               I = II - 1
 458               K = I
 459               P = D( I )
 460               DO 30 J = II, N
 461                  IF( D( J ).LT.P ) THEN
 462                     K = J
 463                     P = D( J )
 464                  END IF
 465    30         CONTINUE
 466               IF( K.NE.I ) THEN
 467                  D( K ) = D( I )
 468                  D( I ) = P
 469                  CALL DSWAP( N, Z( 1, I ), 1, Z( 1, K ), 1 )
 470               END IF
 471    40      CONTINUE
 472          END IF
 473       END IF
 474 *
 475    50 CONTINUE
 476       WORK( 1 ) = LWMIN
 477       IWORK( 1 ) = LIWMIN
 478 *
 479       RETURN
 480 *
 481 *     End of DSTEDC
 482 *
 483       END