SRC/slaed0.f

   1 *> \brief \b SLAED0 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 SLAED0 + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaed0.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaed0.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaed0.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
  22 *                          WORK, IWORK, INFO )
  23 *
  24 *       .. Scalar Arguments ..
  25 *       INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
  26 *       ..
  27 *       .. Array Arguments ..
  28 *       INTEGER            IWORK( * )
  29 *       REAL               D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
  30 *      $                   WORK( * )
  31 *       ..
  32 *
  33 *
  34 *> \par Purpose:
  35 *  =============
  36 *>
  37 *> \verbatim
  38 *>
  39 *> SLAED0 computes all eigenvalues and corresponding eigenvectors of a
  40 *> symmetric tridiagonal matrix using the divide and conquer method.
  41 *> \endverbatim
  42 *
  43 *  Arguments:
  44 *  ==========
  45 *
  46 *> \param[in] ICOMPQ
  47 *> \verbatim
  48 *>          ICOMPQ is INTEGER
  49 *>          = 0:  Compute eigenvalues only.
  50 *>          = 1:  Compute eigenvectors of original dense symmetric matrix
  51 *>                also.  On entry, Q contains the orthogonal matrix used
  52 *>                to reduce the original matrix to tridiagonal form.
  53 *>          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
  54 *>                matrix.
  55 *> \endverbatim
  56 *>
  57 *> \param[in] QSIZ
  58 *> \verbatim
  59 *>          QSIZ is INTEGER
  60 *>         The dimension of the orthogonal matrix used to reduce
  61 *>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
  62 *> \endverbatim
  63 *>
  64 *> \param[in] N
  65 *> \verbatim
  66 *>          N is INTEGER
  67 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
  68 *> \endverbatim
  69 *>
  70 *> \param[in,out] D
  71 *> \verbatim
  72 *>          D is REAL array, dimension (N)
  73 *>         On entry, the main diagonal of the tridiagonal matrix.
  74 *>         On exit, its eigenvalues.
  75 *> \endverbatim
  76 *>
  77 *> \param[in] E
  78 *> \verbatim
  79 *>          E is REAL array, dimension (N-1)
  80 *>         The off-diagonal elements of the tridiagonal matrix.
  81 *>         On exit, E has been destroyed.
  82 *> \endverbatim
  83 *>
  84 *> \param[in,out] Q
  85 *> \verbatim
  86 *>          Q is REAL array, dimension (LDQ, N)
  87 *>         On entry, Q must contain an N-by-N orthogonal matrix.
  88 *>         If ICOMPQ = 0    Q is not referenced.
  89 *>         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
  90 *>                          orthogonal matrix used to reduce the full
  91 *>                          matrix to tridiagonal form corresponding to
  92 *>                          the subset of the full matrix which is being
  93 *>                          decomposed at this time.
  94 *>         If ICOMPQ = 2    On entry, Q will be the identity matrix.
  95 *>                          On exit, Q contains the eigenvectors of the
  96 *>                          tridiagonal matrix.
  97 *> \endverbatim
  98 *>
  99 *> \param[in] LDQ
 100 *> \verbatim
 101 *>          LDQ is INTEGER
 102 *>         The leading dimension of the array Q.  If eigenvectors are
 103 *>         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
 104 *> \endverbatim
 105 *>
 106 *> \param[out] QSTORE
 107 *> \verbatim
 108 *>          QSTORE is REAL array, dimension (LDQS, N)
 109 *>         Referenced only when ICOMPQ = 1.  Used to store parts of
 110 *>         the eigenvector matrix when the updating matrix multiplies
 111 *>         take place.
 112 *> \endverbatim
 113 *>
 114 *> \param[in] LDQS
 115 *> \verbatim
 116 *>          LDQS is INTEGER
 117 *>         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
 118 *>         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
 119 *> \endverbatim
 120 *>
 121 *> \param[out] WORK
 122 *> \verbatim
 123 *>          WORK is REAL array,
 124 *>         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
 125 *>                     1 + 3*N + 2*N*lg N + 3*N**2
 126 *>                     ( lg( N ) = smallest integer k
 127 *>                                 such that 2^k >= N )
 128 *>         If ICOMPQ = 2, the dimension of WORK must be at least
 129 *>                     4*N + N**2.
 130 *> \endverbatim
 131 *>
 132 *> \param[out] IWORK
 133 *> \verbatim
 134 *>          IWORK is INTEGER array,
 135 *>         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
 136 *>                        6 + 6*N + 5*N*lg N.
 137 *>                        ( lg( N ) = smallest integer k
 138 *>                                    such that 2^k >= N )
 139 *>         If ICOMPQ = 2, the dimension of IWORK must be at least
 140 *>                        3 + 5*N.
 141 *> \endverbatim
 142 *>
 143 *> \param[out] INFO
 144 *> \verbatim
 145 *>          INFO is INTEGER
 146 *>          = 0:  successful exit.
 147 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
 148 *>          > 0:  The algorithm failed to compute an eigenvalue while
 149 *>                working on the submatrix lying in rows and columns
 150 *>                INFO/(N+1) through mod(INFO,N+1).
 151 *> \endverbatim
 152 *
 153 *  Authors:
 154 *  ========
 155 *
 156 *> \author Univ. of Tennessee
 157 *> \author Univ. of California Berkeley
 158 *> \author Univ. of Colorado Denver
 159 *> \author NAG Ltd.
 160 *
 161 *> \date September 2012
 162 *
 163 *> \ingroup auxOTHERcomputational
 164 *
 165 *> \par Contributors:
 166 *  ==================
 167 *>
 168 *> Jeff Rutter, Computer Science Division, University of California
 169 *> at Berkeley, USA
 170 *
 171 *  =====================================================================
 172       SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
 173      $                   WORK, IWORK, INFO )
 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 *
 434       END