LAPACK: SRC/zlaed0.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZLAED0( QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS, RWORK,
00002      $                   IWORK, INFO )
00003 *
00004 *  -- LAPACK routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       INTEGER            INFO, LDQ, LDQS, N, QSIZ
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IWORK( * )
00014       DOUBLE PRECISION   D( * ), E( * ), RWORK( * )
00015       COMPLEX*16         Q( LDQ, * ), QSTORE( LDQS, * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  Using the divide and conquer method, ZLAED0 computes all eigenvalues
00022 *  of a symmetric tridiagonal matrix which is one diagonal block of
00023 *  those from reducing a dense or band Hermitian matrix and
00024 *  corresponding eigenvectors of the dense or band matrix.
00025 *
00026 *  Arguments
00027 *  =========
00028 *
00029 *  QSIZ   (input) INTEGER
00030 *         The dimension of the unitary matrix used to reduce
00031 *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
00032 *
00033 *  N      (input) INTEGER
00034 *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00035 *
00036 *  D      (input/output) DOUBLE PRECISION array, dimension (N)
00037 *         On entry, the diagonal elements of the tridiagonal matrix.
00038 *         On exit, the eigenvalues in ascending order.
00039 *
00040 *  E      (input/output) DOUBLE PRECISION array, dimension (N-1)
00041 *         On entry, the off-diagonal elements of the tridiagonal matrix.
00042 *         On exit, E has been destroyed.
00043 *
00044 *  Q      (input/output) COMPLEX*16 array, dimension (LDQ,N)
00045 *         On entry, Q must contain an QSIZ x N matrix whose columns
00046 *         unitarily orthonormal. It is a part of the unitary matrix
00047 *         that reduces the full dense Hermitian matrix to a
00048 *         (reducible) symmetric tridiagonal matrix.
00049 *
00050 *  LDQ    (input) INTEGER
00051 *         The leading dimension of the array Q.  LDQ >= max(1,N).
00052 *
00053 *  IWORK  (workspace) INTEGER array,
00054 *         the dimension of IWORK must be at least
00055 *                      6 + 6*N + 5*N*lg N
00056 *                      ( lg( N ) = smallest integer k
00057 *                                  such that 2^k >= N )
00058 *
00059 *  RWORK  (workspace) DOUBLE PRECISION array,
00060 *                               dimension (1 + 3*N + 2*N*lg N + 3*N**2)
00061 *                        ( lg( N ) = smallest integer k
00062 *                                    such that 2^k >= N )
00063 *
00064 *  QSTORE (workspace) COMPLEX*16 array, dimension (LDQS, N)
00065 *         Used to store parts of
00066 *         the eigenvector matrix when the updating matrix multiplies
00067 *         take place.
00068 *
00069 *  LDQS   (input) INTEGER
00070 *         The leading dimension of the array QSTORE.
00071 *         LDQS >= max(1,N).
00072 *
00073 *  INFO   (output) INTEGER
00074 *          = 0:  successful exit.
00075 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00076 *          > 0:  The algorithm failed to compute an eigenvalue while
00077 *                working on the submatrix lying in rows and columns
00078 *                INFO/(N+1) through mod(INFO,N+1).
00079 *
00080 *  =====================================================================
00081 *
00082 *  Warning:      N could be as big as QSIZ!
00083 *
00084 *     .. Parameters ..
00085       DOUBLE PRECISION   TWO
00086       PARAMETER          ( TWO = 2.D+0 )
00087 *     ..
00088 *     .. Local Scalars ..
00089       INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
00090      $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
00091      $                   J, K, LGN, LL, MATSIZ, MSD2, SMLSIZ, SMM1,
00092      $                   SPM1, SPM2, SUBMAT, SUBPBS, TLVLS
00093       DOUBLE PRECISION   TEMP
00094 *     ..
00095 *     .. External Subroutines ..
00096       EXTERNAL           DCOPY, DSTEQR, XERBLA, ZCOPY, ZLACRM, ZLAED7
00097 *     ..
00098 *     .. External Functions ..
00099       INTEGER            ILAENV
00100       EXTERNAL           ILAENV
00101 *     ..
00102 *     .. Intrinsic Functions ..
00103       INTRINSIC          ABS, DBLE, INT, LOG, MAX
00104 *     ..
00105 *     .. Executable Statements ..
00106 *
00107 *     Test the input parameters.
00108 *
00109       INFO = 0
00110 *
00111 *     IF( ICOMPQ .LT. 0 .OR. ICOMPQ .GT. 2 ) THEN
00112 *        INFO = -1
00113 *     ELSE IF( ( ICOMPQ .EQ. 1 ) .AND. ( QSIZ .LT. MAX( 0, N ) ) )
00114 *    $        THEN
00115       IF( QSIZ.LT.MAX( 0, N ) ) THEN
00116          INFO = -1
00117       ELSE IF( N.LT.0 ) THEN
00118          INFO = -2
00119       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00120          INFO = -6
00121       ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
00122          INFO = -8
00123       END IF
00124       IF( INFO.NE.0 ) THEN
00125          CALL XERBLA( 'ZLAED0', -INFO )
00126          RETURN
00127       END IF
00128 *
00129 *     Quick return if possible
00130 *
00131       IF( N.EQ.0 )
00132      $   RETURN
00133 *
00134       SMLSIZ = ILAENV( 9, 'ZLAED0', ' ', 0, 0, 0, 0 )
00135 *
00136 *     Determine the size and placement of the submatrices, and save in
00137 *     the leading elements of IWORK.
00138 *
00139       IWORK( 1 ) = N
00140       SUBPBS = 1
00141       TLVLS = 0
00142    10 CONTINUE
00143       IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
00144          DO 20 J = SUBPBS, 1, -1
00145             IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
00146             IWORK( 2*J-1 ) = IWORK( J ) / 2
00147    20    CONTINUE
00148          TLVLS = TLVLS + 1
00149          SUBPBS = 2*SUBPBS
00150          GO TO 10
00151       END IF
00152       DO 30 J = 2, SUBPBS
00153          IWORK( J ) = IWORK( J ) + IWORK( J-1 )
00154    30 CONTINUE
00155 *
00156 *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
00157 *     using rank-1 modifications (cuts).
00158 *
00159       SPM1 = SUBPBS - 1
00160       DO 40 I = 1, SPM1
00161          SUBMAT = IWORK( I ) + 1
00162          SMM1 = SUBMAT - 1
00163          D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
00164          D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
00165    40 CONTINUE
00166 *
00167       INDXQ = 4*N + 3
00168 *
00169 *     Set up workspaces for eigenvalues only/accumulate new vectors
00170 *     routine
00171 *
00172       TEMP = LOG( DBLE( N ) ) / LOG( TWO )
00173       LGN = INT( TEMP )
00174       IF( 2**LGN.LT.N )
00175      $   LGN = LGN + 1
00176       IF( 2**LGN.LT.N )
00177      $   LGN = LGN + 1
00178       IPRMPT = INDXQ + N + 1
00179       IPERM = IPRMPT + N*LGN
00180       IQPTR = IPERM + N*LGN
00181       IGIVPT = IQPTR + N + 2
00182       IGIVCL = IGIVPT + N*LGN
00183 *
00184       IGIVNM = 1
00185       IQ = IGIVNM + 2*N*LGN
00186       IWREM = IQ + N**2 + 1
00187 *     Initialize pointers
00188       DO 50 I = 0, SUBPBS
00189          IWORK( IPRMPT+I ) = 1
00190          IWORK( IGIVPT+I ) = 1
00191    50 CONTINUE
00192       IWORK( IQPTR ) = 1
00193 *
00194 *     Solve each submatrix eigenproblem at the bottom of the divide and
00195 *     conquer tree.
00196 *
00197       CURR = 0
00198       DO 70 I = 0, SPM1
00199          IF( I.EQ.0 ) THEN
00200             SUBMAT = 1
00201             MATSIZ = IWORK( 1 )
00202          ELSE
00203             SUBMAT = IWORK( I ) + 1
00204             MATSIZ = IWORK( I+1 ) - IWORK( I )
00205          END IF
00206          LL = IQ - 1 + IWORK( IQPTR+CURR )
00207          CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
00208      $                RWORK( LL ), MATSIZ, RWORK, INFO )
00209          CALL ZLACRM( QSIZ, MATSIZ, Q( 1, SUBMAT ), LDQ, RWORK( LL ),
00210      $                MATSIZ, QSTORE( 1, SUBMAT ), LDQS,
00211      $                RWORK( IWREM ) )
00212          IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
00213          CURR = CURR + 1
00214          IF( INFO.GT.0 ) THEN
00215             INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
00216             RETURN
00217          END IF
00218          K = 1
00219          DO 60 J = SUBMAT, IWORK( I+1 )
00220             IWORK( INDXQ+J ) = K
00221             K = K + 1
00222    60    CONTINUE
00223    70 CONTINUE
00224 *
00225 *     Successively merge eigensystems of adjacent submatrices
00226 *     into eigensystem for the corresponding larger matrix.
00227 *
00228 *     while ( SUBPBS > 1 )
00229 *
00230       CURLVL = 1
00231    80 CONTINUE
00232       IF( SUBPBS.GT.1 ) THEN
00233          SPM2 = SUBPBS - 2
00234          DO 90 I = 0, SPM2, 2
00235             IF( I.EQ.0 ) THEN
00236                SUBMAT = 1
00237                MATSIZ = IWORK( 2 )
00238                MSD2 = IWORK( 1 )
00239                CURPRB = 0
00240             ELSE
00241                SUBMAT = IWORK( I ) + 1
00242                MATSIZ = IWORK( I+2 ) - IWORK( I )
00243                MSD2 = MATSIZ / 2
00244                CURPRB = CURPRB + 1
00245             END IF
00246 *
00247 *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
00248 *     into an eigensystem of size MATSIZ.  ZLAED7 handles the case
00249 *     when the eigenvectors of a full or band Hermitian matrix (which
00250 *     was reduced to tridiagonal form) are desired.
00251 *
00252 *     I am free to use Q as a valuable working space until Loop 150.
00253 *
00254             CALL ZLAED7( MATSIZ, MSD2, QSIZ, TLVLS, CURLVL, CURPRB,
00255      $                   D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
00256      $                   E( SUBMAT+MSD2-1 ), IWORK( INDXQ+SUBMAT ),
00257      $                   RWORK( IQ ), IWORK( IQPTR ), IWORK( IPRMPT ),
00258      $                   IWORK( IPERM ), IWORK( IGIVPT ),
00259      $                   IWORK( IGIVCL ), RWORK( IGIVNM ),
00260      $                   Q( 1, SUBMAT ), RWORK( IWREM ),
00261      $                   IWORK( SUBPBS+1 ), INFO )
00262             IF( INFO.GT.0 ) THEN
00263                INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
00264                RETURN
00265             END IF
00266             IWORK( I / 2+1 ) = IWORK( I+2 )
00267    90    CONTINUE
00268          SUBPBS = SUBPBS / 2
00269          CURLVL = CURLVL + 1
00270          GO TO 80
00271       END IF
00272 *
00273 *     end while
00274 *
00275 *     Re-merge the eigenvalues/vectors which were deflated at the final
00276 *     merge step.
00277 *
00278       DO 100 I = 1, N
00279          J = IWORK( INDXQ+I )
00280          RWORK( I ) = D( J )
00281          CALL ZCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
00282   100 CONTINUE
00283       CALL DCOPY( N, RWORK, 1, D, 1 )
00284 *
00285       RETURN
00286 *
00287 *     End of ZLAED0
00288 *
00289       END