LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaed0.f
Go to the documentation of this file.
00001 *> \brief \b DLAED0
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download DLAED0 + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
00022 *                          WORK, IWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
00026 *       ..
00027 *       .. Array Arguments ..
00028 *       INTEGER            IWORK( * )
00029 *       DOUBLE PRECISION   D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
00030 *      $                   WORK( * )
00031 *       ..
00032 *  
00033 *
00034 *> \par Purpose:
00035 *  =============
00036 *>
00037 *> \verbatim
00038 *>
00039 *> DLAED0 computes all eigenvalues and corresponding eigenvectors of a
00040 *> symmetric tridiagonal matrix using the divide and conquer method.
00041 *> \endverbatim
00042 *
00043 *  Arguments:
00044 *  ==========
00045 *
00046 *> \param[in] ICOMPQ
00047 *> \verbatim
00048 *>          ICOMPQ is INTEGER
00049 *>          = 0:  Compute eigenvalues only.
00050 *>          = 1:  Compute eigenvectors of original dense symmetric matrix
00051 *>                also.  On entry, Q contains the orthogonal matrix used
00052 *>                to reduce the original matrix to tridiagonal form.
00053 *>          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
00054 *>                matrix.
00055 *> \endverbatim
00056 *>
00057 *> \param[in] QSIZ
00058 *> \verbatim
00059 *>          QSIZ is INTEGER
00060 *>         The dimension of the orthogonal matrix used to reduce
00061 *>         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] N
00065 *> \verbatim
00066 *>          N is INTEGER
00067 *>         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00068 *> \endverbatim
00069 *>
00070 *> \param[in,out] D
00071 *> \verbatim
00072 *>          D is DOUBLE PRECISION array, dimension (N)
00073 *>         On entry, the main diagonal of the tridiagonal matrix.
00074 *>         On exit, its eigenvalues.
00075 *> \endverbatim
00076 *>
00077 *> \param[in] E
00078 *> \verbatim
00079 *>          E is DOUBLE PRECISION array, dimension (N-1)
00080 *>         The off-diagonal elements of the tridiagonal matrix.
00081 *>         On exit, E has been destroyed.
00082 *> \endverbatim
00083 *>
00084 *> \param[in,out] Q
00085 *> \verbatim
00086 *>          Q is DOUBLE PRECISION array, dimension (LDQ, N)
00087 *>         On entry, Q must contain an N-by-N orthogonal matrix.
00088 *>         If ICOMPQ = 0    Q is not referenced.
00089 *>         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
00090 *>                          orthogonal matrix used to reduce the full
00091 *>                          matrix to tridiagonal form corresponding to
00092 *>                          the subset of the full matrix which is being
00093 *>                          decomposed at this time.
00094 *>         If ICOMPQ = 2    On entry, Q will be the identity matrix.
00095 *>                          On exit, Q contains the eigenvectors of the
00096 *>                          tridiagonal matrix.
00097 *> \endverbatim
00098 *>
00099 *> \param[in] LDQ
00100 *> \verbatim
00101 *>          LDQ is INTEGER
00102 *>         The leading dimension of the array Q.  If eigenvectors are
00103 *>         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
00104 *> \endverbatim
00105 *>
00106 *> \param[out] QSTORE
00107 *> \verbatim
00108 *>          QSTORE is DOUBLE PRECISION array, dimension (LDQS, N)
00109 *>         Referenced only when ICOMPQ = 1.  Used to store parts of
00110 *>         the eigenvector matrix when the updating matrix multiplies
00111 *>         take place.
00112 *> \endverbatim
00113 *>
00114 *> \param[in] LDQS
00115 *> \verbatim
00116 *>          LDQS is INTEGER
00117 *>         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
00118 *>         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
00119 *> \endverbatim
00120 *>
00121 *> \param[out] WORK
00122 *> \verbatim
00123 *>          WORK is DOUBLE PRECISION array,
00124 *>         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
00125 *>                     1 + 3*N + 2*N*lg N + 3*N**2
00126 *>                     ( lg( N ) = smallest integer k
00127 *>                                 such that 2^k >= N )
00128 *>         If ICOMPQ = 2, the dimension of WORK must be at least
00129 *>                     4*N + N**2.
00130 *> \endverbatim
00131 *>
00132 *> \param[out] IWORK
00133 *> \verbatim
00134 *>          IWORK is INTEGER array,
00135 *>         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
00136 *>                        6 + 6*N + 5*N*lg N.
00137 *>                        ( lg( N ) = smallest integer k
00138 *>                                    such that 2^k >= N )
00139 *>         If ICOMPQ = 2, the dimension of IWORK must be at least
00140 *>                        3 + 5*N.
00141 *> \endverbatim
00142 *>
00143 *> \param[out] INFO
00144 *> \verbatim
00145 *>          INFO is INTEGER
00146 *>          = 0:  successful exit.
00147 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00148 *>          > 0:  The algorithm failed to compute an eigenvalue while
00149 *>                working on the submatrix lying in rows and columns
00150 *>                INFO/(N+1) through mod(INFO,N+1).
00151 *> \endverbatim
00152 *
00153 *  Authors:
00154 *  ========
00155 *
00156 *> \author Univ. of Tennessee 
00157 *> \author Univ. of California Berkeley 
00158 *> \author Univ. of Colorado Denver 
00159 *> \author NAG Ltd. 
00160 *
00161 *> \date November 2011
00162 *
00163 *> \ingroup auxOTHERcomputational
00164 *
00165 *> \par Contributors:
00166 *  ==================
00167 *>
00168 *> Jeff Rutter, Computer Science Division, University of California
00169 *> at Berkeley, USA
00170 *
00171 *  =====================================================================
00172       SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
00173      $                   WORK, IWORK, INFO )
00174 *
00175 *  -- LAPACK computational routine (version 3.4.0) --
00176 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00177 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00178 *     November 2011
00179 *
00180 *     .. Scalar Arguments ..
00181       INTEGER            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
00182 *     ..
00183 *     .. Array Arguments ..
00184       INTEGER            IWORK( * )
00185       DOUBLE PRECISION   D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
00186      $                   WORK( * )
00187 *     ..
00188 *
00189 *  =====================================================================
00190 *
00191 *     .. Parameters ..
00192       DOUBLE PRECISION   ZERO, ONE, TWO
00193       PARAMETER          ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
00194 *     ..
00195 *     .. Local Scalars ..
00196       INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
00197      $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
00198      $                   J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
00199      $                   SPM2, SUBMAT, SUBPBS, TLVLS
00200       DOUBLE PRECISION   TEMP
00201 *     ..
00202 *     .. External Subroutines ..
00203       EXTERNAL           DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
00204      $                   XERBLA
00205 *     ..
00206 *     .. External Functions ..
00207       INTEGER            ILAENV
00208       EXTERNAL           ILAENV
00209 *     ..
00210 *     .. Intrinsic Functions ..
00211       INTRINSIC          ABS, DBLE, INT, LOG, MAX
00212 *     ..
00213 *     .. Executable Statements ..
00214 *
00215 *     Test the input parameters.
00216 *
00217       INFO = 0
00218 *
00219       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
00220          INFO = -1
00221       ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
00222          INFO = -2
00223       ELSE IF( N.LT.0 ) THEN
00224          INFO = -3
00225       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00226          INFO = -7
00227       ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
00228          INFO = -9
00229       END IF
00230       IF( INFO.NE.0 ) THEN
00231          CALL XERBLA( 'DLAED0', -INFO )
00232          RETURN
00233       END IF
00234 *
00235 *     Quick return if possible
00236 *
00237       IF( N.EQ.0 )
00238      $   RETURN
00239 *
00240       SMLSIZ = ILAENV( 9, 'DLAED0', ' ', 0, 0, 0, 0 )
00241 *
00242 *     Determine the size and placement of the submatrices, and save in
00243 *     the leading elements of IWORK.
00244 *
00245       IWORK( 1 ) = N
00246       SUBPBS = 1
00247       TLVLS = 0
00248    10 CONTINUE
00249       IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
00250          DO 20 J = SUBPBS, 1, -1
00251             IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
00252             IWORK( 2*J-1 ) = IWORK( J ) / 2
00253    20    CONTINUE
00254          TLVLS = TLVLS + 1
00255          SUBPBS = 2*SUBPBS
00256          GO TO 10
00257       END IF
00258       DO 30 J = 2, SUBPBS
00259          IWORK( J ) = IWORK( J ) + IWORK( J-1 )
00260    30 CONTINUE
00261 *
00262 *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
00263 *     using rank-1 modifications (cuts).
00264 *
00265       SPM1 = SUBPBS - 1
00266       DO 40 I = 1, SPM1
00267          SUBMAT = IWORK( I ) + 1
00268          SMM1 = SUBMAT - 1
00269          D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
00270          D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
00271    40 CONTINUE
00272 *
00273       INDXQ = 4*N + 3
00274       IF( ICOMPQ.NE.2 ) THEN
00275 *
00276 *        Set up workspaces for eigenvalues only/accumulate new vectors
00277 *        routine
00278 *
00279          TEMP = LOG( DBLE( N ) ) / LOG( TWO )
00280          LGN = INT( TEMP )
00281          IF( 2**LGN.LT.N )
00282      $      LGN = LGN + 1
00283          IF( 2**LGN.LT.N )
00284      $      LGN = LGN + 1
00285          IPRMPT = INDXQ + N + 1
00286          IPERM = IPRMPT + N*LGN
00287          IQPTR = IPERM + N*LGN
00288          IGIVPT = IQPTR + N + 2
00289          IGIVCL = IGIVPT + N*LGN
00290 *
00291          IGIVNM = 1
00292          IQ = IGIVNM + 2*N*LGN
00293          IWREM = IQ + N**2 + 1
00294 *
00295 *        Initialize pointers
00296 *
00297          DO 50 I = 0, SUBPBS
00298             IWORK( IPRMPT+I ) = 1
00299             IWORK( IGIVPT+I ) = 1
00300    50    CONTINUE
00301          IWORK( IQPTR ) = 1
00302       END IF
00303 *
00304 *     Solve each submatrix eigenproblem at the bottom of the divide and
00305 *     conquer tree.
00306 *
00307       CURR = 0
00308       DO 70 I = 0, SPM1
00309          IF( I.EQ.0 ) THEN
00310             SUBMAT = 1
00311             MATSIZ = IWORK( 1 )
00312          ELSE
00313             SUBMAT = IWORK( I ) + 1
00314             MATSIZ = IWORK( I+1 ) - IWORK( I )
00315          END IF
00316          IF( ICOMPQ.EQ.2 ) THEN
00317             CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
00318      $                   Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
00319             IF( INFO.NE.0 )
00320      $         GO TO 130
00321          ELSE
00322             CALL DSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
00323      $                   WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
00324      $                   INFO )
00325             IF( INFO.NE.0 )
00326      $         GO TO 130
00327             IF( ICOMPQ.EQ.1 ) THEN
00328                CALL DGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
00329      $                     Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
00330      $                     CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
00331      $                     LDQS )
00332             END IF
00333             IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
00334             CURR = CURR + 1
00335          END IF
00336          K = 1
00337          DO 60 J = SUBMAT, IWORK( I+1 )
00338             IWORK( INDXQ+J ) = K
00339             K = K + 1
00340    60    CONTINUE
00341    70 CONTINUE
00342 *
00343 *     Successively merge eigensystems of adjacent submatrices
00344 *     into eigensystem for the corresponding larger matrix.
00345 *
00346 *     while ( SUBPBS > 1 )
00347 *
00348       CURLVL = 1
00349    80 CONTINUE
00350       IF( SUBPBS.GT.1 ) THEN
00351          SPM2 = SUBPBS - 2
00352          DO 90 I = 0, SPM2, 2
00353             IF( I.EQ.0 ) THEN
00354                SUBMAT = 1
00355                MATSIZ = IWORK( 2 )
00356                MSD2 = IWORK( 1 )
00357                CURPRB = 0
00358             ELSE
00359                SUBMAT = IWORK( I ) + 1
00360                MATSIZ = IWORK( I+2 ) - IWORK( I )
00361                MSD2 = MATSIZ / 2
00362                CURPRB = CURPRB + 1
00363             END IF
00364 *
00365 *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
00366 *     into an eigensystem of size MATSIZ.
00367 *     DLAED1 is used only for the full eigensystem of a tridiagonal
00368 *     matrix.
00369 *     DLAED7 handles the cases in which eigenvalues only or eigenvalues
00370 *     and eigenvectors of a full symmetric matrix (which was reduced to
00371 *     tridiagonal form) are desired.
00372 *
00373             IF( ICOMPQ.EQ.2 ) THEN
00374                CALL DLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
00375      $                      LDQ, IWORK( INDXQ+SUBMAT ),
00376      $                      E( SUBMAT+MSD2-1 ), MSD2, WORK,
00377      $                      IWORK( SUBPBS+1 ), INFO )
00378             ELSE
00379                CALL DLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
00380      $                      D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
00381      $                      IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
00382      $                      MSD2, WORK( IQ ), IWORK( IQPTR ),
00383      $                      IWORK( IPRMPT ), IWORK( IPERM ),
00384      $                      IWORK( IGIVPT ), IWORK( IGIVCL ),
00385      $                      WORK( IGIVNM ), WORK( IWREM ),
00386      $                      IWORK( SUBPBS+1 ), INFO )
00387             END IF
00388             IF( INFO.NE.0 )
00389      $         GO TO 130
00390             IWORK( I / 2+1 ) = IWORK( I+2 )
00391    90    CONTINUE
00392          SUBPBS = SUBPBS / 2
00393          CURLVL = CURLVL + 1
00394          GO TO 80
00395       END IF
00396 *
00397 *     end while
00398 *
00399 *     Re-merge the eigenvalues/vectors which were deflated at the final
00400 *     merge step.
00401 *
00402       IF( ICOMPQ.EQ.1 ) THEN
00403          DO 100 I = 1, N
00404             J = IWORK( INDXQ+I )
00405             WORK( I ) = D( J )
00406             CALL DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
00407   100    CONTINUE
00408          CALL DCOPY( N, WORK, 1, D, 1 )
00409       ELSE IF( ICOMPQ.EQ.2 ) THEN
00410          DO 110 I = 1, N
00411             J = IWORK( INDXQ+I )
00412             WORK( I ) = D( J )
00413             CALL DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
00414   110    CONTINUE
00415          CALL DCOPY( N, WORK, 1, D, 1 )
00416          CALL DLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
00417       ELSE
00418          DO 120 I = 1, N
00419             J = IWORK( INDXQ+I )
00420             WORK( I ) = D( J )
00421   120    CONTINUE
00422          CALL DCOPY( N, WORK, 1, D, 1 )
00423       END IF
00424       GO TO 140
00425 *
00426   130 CONTINUE
00427       INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
00428 *
00429   140 CONTINUE
00430       RETURN
00431 *
00432 *     End of DLAED0
00433 *
00434       END
 All Files Functions