SRC/zlaed0.f

   1 *> \brief \b ZLAED0 used by sstedc. Computes all eigenvalues and corresponding eigenvectors of an unreduced symmetric tridiagonal matrix using the divide and conquer method.
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *> \htmlonly
   9 *> Download ZLAED0 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlaed0.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlaed0.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlaed0.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
  22 *                          IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            INFO, LDQ, LDQS, N, QSIZ
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IWORK( * )
  29 *       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
  30 *       COMPLEX*16         Q( LDQ, * ), QSTORE( LDQS, * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> Using the divide and conquer method, ZLAED0 computes all eigenvalues
  40 *> of a symmetric tridiagonal matrix which is one diagonal block of
  41 *> those from reducing a dense or band Hermitian matrix and
  42 *> corresponding eigenvectors of the dense or band matrix.
  43 *> \endverbatim
  44 *
  45 *  Arguments:
  46 *  ==========
  47 *
  48 *> \param[in] QSIZ
  49 *> \verbatim
  50 *>          QSIZ is INTEGER
  51 *>         The dimension of the unitary matrix used to reduce
  52 *>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
  53 *> \endverbatim
  54 *>
  55 *> \param[in] N
  56 *> \verbatim
  57 *>          N is INTEGER
  58 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
  59 *> \endverbatim
  60 *>
  61 *> \param[in,out] D
  62 *> \verbatim
  63 *>          D is DOUBLE PRECISION array, dimension (N)
  64 *>         On entry, the diagonal elements of the tridiagonal matrix.
  65 *>         On exit, the eigenvalues in ascending order.
  66 *> \endverbatim
  67 *>
  68 *> \param[in,out] E
  69 *> \verbatim
  70 *>          E is DOUBLE PRECISION array, dimension (N-1)
  71 *>         On entry, the off-diagonal elements of the tridiagonal matrix.
  72 *>         On exit, E has been destroyed.
  73 *> \endverbatim
  74 *>
  75 *> \param[in,out] Q
  76 *> \verbatim
  77 *>          Q is COMPLEX*16 array, dimension (LDQ,N)
  78 *>         On entry, Q must contain an QSIZ x N matrix whose columns
  79 *>         unitarily orthonormal. It is a part of the unitary matrix
  80 *>         that reduces the full dense Hermitian matrix to a
  81 *>         (reducible) symmetric tridiagonal matrix.
  82 *> \endverbatim
  83 *>
  84 *> \param[in] LDQ
  85 *> \verbatim
  86 *>          LDQ is INTEGER
  87 *>         The leading dimension of the array Q.  LDQ >= max(1,N).
  88 *> \endverbatim
  89 *>
  90 *> \param[out] IWORK
  91 *> \verbatim
  92 *>          IWORK is INTEGER array,
  93 *>         the dimension of IWORK must be at least
  94 *>                      6 + 6*N + 5*N*lg N
  95 *>                      ( lg( N ) = smallest integer k
  96 *>                                  such that 2^k >= N )
  97 *> \endverbatim
  98 *>
  99 *> \param[out] RWORK
 100 *> \verbatim
 101 *>          RWORK is DOUBLE PRECISION array,
 102 *>                               dimension (1 + 3*N + 2*N*lg N + 3*N**2)
 103 *>                        ( lg( N ) = smallest integer k
 104 *>                                    such that 2^k >= N )
 105 *> \endverbatim
 106 *>
 107 *> \param[out] QSTORE
 108 *> \verbatim
 109 *>          QSTORE is COMPLEX*16 array, dimension (LDQS, N)
 110 *>         Used to store parts of
 111 *>         the eigenvector matrix when the updating matrix multiplies
 112 *>         take place.
 113 *> \endverbatim
 114 *>
 115 *> \param[in] LDQS
 116 *> \verbatim
 117 *>          LDQS is INTEGER
 118 *>         The leading dimension of the array QSTORE.
 119 *>         LDQS >= max(1,N).
 120 *> \endverbatim
 121 *>
 122 *> \param[out] INFO
 123 *> \verbatim
 124 *>          INFO is INTEGER
 125 *>          = 0:  successful exit.
 126 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 127 *>          > 0:  The algorithm failed to compute an eigenvalue while
 128 *>                working on the submatrix lying in rows and columns
 129 *>                INFO/(N+1) through mod(INFO,N+1).
 130 *> \endverbatim
 131 *
 132 *  Authors:
 133 *  ========
 134 *
 135 *> \author Univ. of Tennessee
 136 *> \author Univ. of California Berkeley
 137 *> \author Univ. of Colorado Denver
 138 *> \author NAG Ltd.
 139 *
 140 *> \date September 2012
 141 *
 142 *> \ingroup complex16OTHERcomputational
 143 *
 144 *  =====================================================================
 145       SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
 146      $                   IWORK, INFO )
 147 *
 148 *  -- LAPACK computational routine (version 3.4.2) --
 149 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 150 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 151 *     September 2012
 152 *
 153 *     .. Scalar Arguments ..
 154       INTEGER            INFO, LDQ, LDQS, N, QSIZ
 155 *     ..
 156 *     .. Array Arguments ..
 157       INTEGER            IWORK( * )
 158       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
 159       COMPLEX*16         Q( LDQ, * ), QSTORE( LDQS, * )
 160 *     ..
 161 *
 162 *  =====================================================================
 163 *
 164 *  Warning:      N could be as big as QSIZ!
 165 *
 166 *     .. Parameters ..
 167       DOUBLE PRECISION   TWO
 168       PARAMETER          ( TWO = 2.D+0 )
 169 *     ..
 170 *     .. Local Scalars ..
 171       INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
 172      $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
 173      $                   J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
 174      $                   SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
 175       DOUBLE PRECISION   TEMP
 176 *     ..
 177 *     .. External Subroutines ..
 178       EXTERNAL           DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
 179 *     ..
 180 *     .. External Functions ..
 181       INTEGER            ILAENV
 182       EXTERNAL           ILAENV
 183 *     ..
 184 *     .. Intrinsic Functions ..
 185       INTRINSIC          ABS, DBLE, INT, LOG, MAX
 186 *     ..
 187 *     .. Executable Statements ..
 188 *
 189 *     Test the input parameters.
 190 *
 191       INFO = 0
 192 *
 193 *     IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
 194 *        INFO = -1
 195 *     ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
 196 *    $        THEN
 197       IF( QSIZ.LT.MAX( 0, N ) ) THEN
 198          INFO = -1
 199       ELSE IF( N.LT.0 ) THEN
 200          INFO = -2
 201       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
 202          INFO = -6
 203       ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
 204          INFO = -8
 205       END IF
 206       IF( INFO.NE.0 ) THEN
 207          CALL XERBLA( 'ZLAED0', -INFO )
 208          RETURN
 209       END IF
 210 *
 211 *     Quick return if possible
 212 *
 213       IF( N.EQ.0 )
 214      $   RETURN
 215 *
 216       SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
 217 *
 218 *     Determine the size and placement of the submatrices, and save in
 219 *     the leading elements of IWORK.
 220 *
 221       IWORK( 1 ) = N
 222       SUBPBS = 1
 223       TLVLS = 0
 224    10 CONTINUE
 225       IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
 226          DO 20 J = SUBPBS, 1, -1
 227             IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
 228             IWORK( 2*J-1 ) = IWORK( J ) / 2
 229    20    CONTINUE
 230          TLVLS = TLVLS + 1
 231          SUBPBS = 2*SUBPBS
 232          GO TO 10
 233       END IF
 234       DO 30 J = 2, SUBPBS
 235          IWORK( J ) = IWORK( J ) + IWORK( J-1 )
 236    30 CONTINUE
 237 *
 238 *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
 239 *     using rank-1 modifications (cuts).
 240 *
 241       SPM1 = SUBPBS - 1
 242       DO 40 I = 1, SPM1
 243          SUBMAT = IWORK( I ) + 1
 244          SMM1 = SUBMAT - 1
 245          D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
 246          D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
 247    40 CONTINUE
 248 *
 249       INDXQ = 4*N + 3
 250 *
 251 *     Set up workspaces for eigenvalues only/accumulate new vectors
 252 *     routine
 253 *
 254       TEMP = LOG( DBLE( N ) ) / LOG( TWO )
 255       LGN = INT( TEMP )
 256       IF( 2**LGN.LT.N )
 257      $   LGN = LGN + 1
 258       IF( 2**LGN.LT.N )
 259      $   LGN = LGN + 1
 260       IPRMPT = INDXQ + N + 1
 261       IPERM = IPRMPT + N*LGN
 262       IQPTR = IPERM + N*LGN
 263       IGIVPT = IQPTR + N + 2
 264       IGIVCL = IGIVPT + N*LGN
 265 *
 266       IGIVNM = 1
 267       IQ = IGIVNM + 2*N*LGN
 268       IWREM = IQ + N**2 + 1
 269 *     Initialize pointers
 270       DO 50 I = 0, SUBPBS
 271          IWORK( IPRMPT+I ) = 1
 272          IWORK( IGIVPT+I ) = 1
 273    50 CONTINUE
 274       IWORK( IQPTR ) = 1
 275 *
 276 *     Solve each submatrix eigenproblem at the bottom of the divide and
 277 *     conquer tree.
 278 *
 279       CURR = 0
 280       DO 70 I = 0, SPM1
 281          IF( I.EQ.0 ) THEN
 282             SUBMAT = 1
 283             MATSIZ = IWORK( 1 )
 284          ELSE
 285             SUBMAT = IWORK( I ) + 1
 286             MATSIZ = IWORK( I+1 ) - IWORK( I )
 287          END IF
 288          LL = IQ - 1 + IWORK( IQPTR+CURR )
 289          CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
 290      $                RWORK( LL ), MATSIZ, RWORK, INFO )
 291          CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
 292      $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
 293      $                RWORK( IWREM ) )
 294          IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
 295          CURR = CURR + 1
 296          IF( INFO.GT.0 ) THEN
 297             INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
 298             RETURN
 299          END IF
 300          K = 1
 301          DO 60 J = SUBMAT, IWORK( I+1 )
 302             IWORK( INDXQ+J ) = K
 303             K = K + 1
 304    60    CONTINUE
 305    70 CONTINUE
 306 *
 307 *     Successively merge eigensystems of adjacent submatrices
 308 *     into eigensystem for the corresponding larger matrix.
 309 *
 310 *     while ( SUBPBS > 1 )
 311 *
 312       CURLVL = 1
 313    80 CONTINUE
 314       IF( SUBPBS.GT.1 ) THEN
 315          SPM2 = SUBPBS - 2
 316          DO 90 I = 0, SPM2, 2
 317             IF( I.EQ.0 ) THEN
 318                SUBMAT = 1
 319                MATSIZ = IWORK( 2 )
 320                MSD2 = IWORK( 1 )
 321                CURPRB = 0
 322             ELSE
 323                SUBMAT = IWORK( I ) + 1
 324                MATSIZ = IWORK( I+2 ) - IWORK( I )
 325                MSD2 = MATSIZ / 2
 326                CURPRB = CURPRB + 1
 327             END IF
 328 *
 329 *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
 330 *     into an eigensystem of size MATSIZ.  ZLAED7 handles the case
 331 *     when the eigenvectors of a full or band Hermitian matrix (which
 332 *     was reduced to tridiagonal form) are desired.
 333 *
 334 *     I am free to use Q as a valuable working space until Loop 150.
 335 *
 336             CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
 337      $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
 338      $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
 339      $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
 340      $                   IWORK( IPERM ), IWORK( IGIVPT ),
 341      $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
 342      $                   Q( 1, SUBMAT ), RWORK( IWREM ),
 343      $                   IWORK( SUBPBS+1 ), INFO )
 344             IF( INFO.GT.0 ) THEN
 345                INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
 346                RETURN
 347             END IF
 348             IWORK( I / 2+1 ) = IWORK( I+2 )
 349    90    CONTINUE
 350          SUBPBS = SUBPBS / 2
 351          CURLVL = CURLVL + 1
 352          GO TO 80
 353       END IF
 354 *
 355 *     end while
 356 *
 357 *     Re-merge the eigenvalues/vectors which were deflated at the final
 358 *     merge step.
 359 *
 360       DO 100 I = 1, N
 361          J = IWORK( INDXQ+I )
 362          RWORK( I ) = D( J )
 363          CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
 364   100 CONTINUE
 365       CALL DCOPY( N, RWORK, 1, D, 1 )
 366 *
 367       RETURN
 368 *
 369 *     End of ZLAED0
 370 *
 371       END