SRC/cstedc.f

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