LAPACK: SRC/slaed0.f Source File

Go to the documentation of this file.
00001       SUBROUTINE SLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
00002      $                   WORK, 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            ICOMPQ, INFO, LDQ, LDQS, N, QSIZ
00011 *     ..
00012 *     .. Array Arguments ..
00013       INTEGER            IWORK( * )
00014       REAL               D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
00015      $                   WORK( * )
00016 *     ..
00017 *
00018 *  Purpose
00019 *  =======
00020 *
00021 *  SLAED0 computes all eigenvalues and corresponding eigenvectors of a
00022 *  symmetric tridiagonal matrix using the divide and conquer method.
00023 *
00024 *  Arguments
00025 *  =========
00026 *
00027 *  ICOMPQ  (input) INTEGER
00028 *          = 0:  Compute eigenvalues only.
00029 *          = 1:  Compute eigenvectors of original dense symmetric matrix
00030 *                also.  On entry, Q contains the orthogonal matrix used
00031 *                to reduce the original matrix to tridiagonal form.
00032 *          = 2:  Compute eigenvalues and eigenvectors of tridiagonal
00033 *                matrix.
00034 *
00035 *  QSIZ   (input) INTEGER
00036 *         The dimension of the orthogonal matrix used to reduce
00037 *         the full matrix to tridiagonal form.  QSIZ >= N if ICOMPQ = 1.
00038 *
00039 *  N      (input) INTEGER
00040 *         The dimension of the symmetric tridiagonal matrix.  N >= 0.
00041 *
00042 *  D      (input/output) REAL array, dimension (N)
00043 *         On entry, the main diagonal of the tridiagonal matrix.
00044 *         On exit, its eigenvalues.
00045 *
00046 *  E      (input) REAL array, dimension (N-1)
00047 *         The off-diagonal elements of the tridiagonal matrix.
00048 *         On exit, E has been destroyed.
00049 *
00050 *  Q      (input/output) REAL array, dimension (LDQ, N)
00051 *         On entry, Q must contain an N-by-N orthogonal matrix.
00052 *         If ICOMPQ = 0    Q is not referenced.
00053 *         If ICOMPQ = 1    On entry, Q is a subset of the columns of the
00054 *                          orthogonal matrix used to reduce the full
00055 *                          matrix to tridiagonal form corresponding to
00056 *                          the subset of the full matrix which is being
00057 *                          decomposed at this time.
00058 *         If ICOMPQ = 2    On entry, Q will be the identity matrix.
00059 *                          On exit, Q contains the eigenvectors of the
00060 *                          tridiagonal matrix.
00061 *
00062 *  LDQ    (input) INTEGER
00063 *         The leading dimension of the array Q.  If eigenvectors are
00064 *         desired, then  LDQ >= max(1,N).  In any case,  LDQ >= 1.
00065 *
00066 *  QSTORE (workspace) REAL array, dimension (LDQS, N)
00067 *         Referenced only when ICOMPQ = 1.  Used to store parts of
00068 *         the eigenvector matrix when the updating matrix multiplies
00069 *         take place.
00070 *
00071 *  LDQS   (input) INTEGER
00072 *         The leading dimension of the array QSTORE.  If ICOMPQ = 1,
00073 *         then  LDQS >= max(1,N).  In any case,  LDQS >= 1.
00074 *
00075 *  WORK   (workspace) REAL array,
00076 *         If ICOMPQ = 0 or 1, the dimension of WORK must be at least
00077 *                     1 + 3*N + 2*N*lg N + 2*N**2
00078 *                     ( lg( N ) = smallest integer k
00079 *                                 such that 2^k >= N )
00080 *         If ICOMPQ = 2, the dimension of WORK must be at least
00081 *                     4*N + N**2.
00082 *
00083 *  IWORK  (workspace) INTEGER array,
00084 *         If ICOMPQ = 0 or 1, the dimension of IWORK must be at least
00085 *                        6 + 6*N + 5*N*lg N.
00086 *                        ( lg( N ) = smallest integer k
00087 *                                    such that 2^k >= N )
00088 *         If ICOMPQ = 2, the dimension of IWORK must be at least
00089 *                        3 + 5*N.
00090 *
00091 *  INFO   (output) INTEGER
00092 *          = 0:  successful exit.
00093 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00094 *          > 0:  The algorithm failed to compute an eigenvalue while
00095 *                working on the submatrix lying in rows and columns
00096 *                INFO/(N+1) through mod(INFO,N+1).
00097 *
00098 *  Further Details
00099 *  ===============
00100 *
00101 *  Based on contributions by
00102 *     Jeff Rutter, Computer Science Division, University of California
00103 *     at Berkeley, USA
00104 *
00105 *  =====================================================================
00106 *
00107 *     .. Parameters ..
00108       REAL               ZERO, ONE, TWO
00109       PARAMETER          ( ZERO = 0.E0, ONE = 1.E0, TWO = 2.E0 )
00110 *     ..
00111 *     .. Local Scalars ..
00112       INTEGER            CURLVL, CURPRB, CURR, I, IGIVCL, IGIVNM,
00113      $                   IGIVPT, INDXQ, IPERM, IPRMPT, IQ, IQPTR, IWREM,
00114      $                   J, K, LGN, MATSIZ, MSD2, SMLSIZ, SMM1, SPM1,
00115      $                   SPM2, SUBMAT, SUBPBS, TLVLS
00116       REAL               TEMP
00117 *     ..
00118 *     .. External Subroutines ..
00119       EXTERNAL           SCOPY, SGEMM, SLACPY, SLAED1, SLAED7, SSTEQR,
00120      $                   XERBLA
00121 *     ..
00122 *     .. External Functions ..
00123       INTEGER            ILAENV
00124       EXTERNAL           ILAENV
00125 *     ..
00126 *     .. Intrinsic Functions ..
00127       INTRINSIC          ABS, INT, LOG, MAX, REAL
00128 *     ..
00129 *     .. Executable Statements ..
00130 *
00131 *     Test the input parameters.
00132 *
00133       INFO = 0
00134 *
00135       IF( ICOMPQ.LT.0 .OR. ICOMPQ.GT.2 ) THEN
00136          INFO = -1
00137       ELSE IF( ( ICOMPQ.EQ.1 ) .AND. ( QSIZ.LT.MAX( 0, N ) ) ) THEN
00138          INFO = -2
00139       ELSE IF( N.LT.0 ) THEN
00140          INFO = -3
00141       ELSE IF( LDQ.LT.MAX( 1, N ) ) THEN
00142          INFO = -7
00143       ELSE IF( LDQS.LT.MAX( 1, N ) ) THEN
00144          INFO = -9
00145       END IF
00146       IF( INFO.NE.0 ) THEN
00147          CALL XERBLA( 'SLAED0', -INFO )
00148          RETURN
00149       END IF
00150 *
00151 *     Quick return if possible
00152 *
00153       IF( N.EQ.0 )
00154      $   RETURN
00155 *
00156       SMLSIZ = ILAENV( 9, 'SLAED0', ' ', 0, 0, 0, 0 )
00157 *
00158 *     Determine the size and placement of the submatrices, and save in
00159 *     the leading elements of IWORK.
00160 *
00161       IWORK( 1 ) = N
00162       SUBPBS = 1
00163       TLVLS = 0
00164    10 CONTINUE
00165       IF( IWORK( SUBPBS ).GT.SMLSIZ ) THEN
00166          DO 20 J = SUBPBS, 1, -1
00167             IWORK( 2*J ) = ( IWORK( J )+1 ) / 2
00168             IWORK( 2*J-1 ) = IWORK( J ) / 2
00169    20    CONTINUE
00170          TLVLS = TLVLS + 1
00171          SUBPBS = 2*SUBPBS
00172          GO TO 10
00173       END IF
00174       DO 30 J = 2, SUBPBS
00175          IWORK( J ) = IWORK( J ) + IWORK( J-1 )
00176    30 CONTINUE
00177 *
00178 *     Divide the matrix into SUBPBS submatrices of size at most SMLSIZ+1
00179 *     using rank-1 modifications (cuts).
00180 *
00181       SPM1 = SUBPBS - 1
00182       DO 40 I = 1, SPM1
00183          SUBMAT = IWORK( I ) + 1
00184          SMM1 = SUBMAT - 1
00185          D( SMM1 ) = D( SMM1 ) - ABS( E( SMM1 ) )
00186          D( SUBMAT ) = D( SUBMAT ) - ABS( E( SMM1 ) )
00187    40 CONTINUE
00188 *
00189       INDXQ = 4*N + 3
00190       IF( ICOMPQ.NE.2 ) THEN
00191 *
00192 *        Set up workspaces for eigenvalues only/accumulate new vectors
00193 *        routine
00194 *
00195          TEMP = LOG( REAL( N ) ) / LOG( TWO )
00196          LGN = INT( TEMP )
00197          IF( 2**LGN.LT.N )
00198      $      LGN = LGN + 1
00199          IF( 2**LGN.LT.N )
00200      $      LGN = LGN + 1
00201          IPRMPT = INDXQ + N + 1
00202          IPERM = IPRMPT + N*LGN
00203          IQPTR = IPERM + N*LGN
00204          IGIVPT = IQPTR + N + 2
00205          IGIVCL = IGIVPT + N*LGN
00206 *
00207          IGIVNM = 1
00208          IQ = IGIVNM + 2*N*LGN
00209          IWREM = IQ + N**2 + 1
00210 *
00211 *        Initialize pointers
00212 *
00213          DO 50 I = 0, SUBPBS
00214             IWORK( IPRMPT+I ) = 1
00215             IWORK( IGIVPT+I ) = 1
00216    50    CONTINUE
00217          IWORK( IQPTR ) = 1
00218       END IF
00219 *
00220 *     Solve each submatrix eigenproblem at the bottom of the divide and
00221 *     conquer tree.
00222 *
00223       CURR = 0
00224       DO 70 I = 0, SPM1
00225          IF( I.EQ.0 ) THEN
00226             SUBMAT = 1
00227             MATSIZ = IWORK( 1 )
00228          ELSE
00229             SUBMAT = IWORK( I ) + 1
00230             MATSIZ = IWORK( I+1 ) - IWORK( I )
00231          END IF
00232          IF( ICOMPQ.EQ.2 ) THEN
00233             CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
00234      $                   Q( SUBMAT, SUBMAT ), LDQ, WORK, INFO )
00235             IF( INFO.NE.0 )
00236      $         GO TO 130
00237          ELSE
00238             CALL SSTEQR( 'I', MATSIZ, D( SUBMAT ), E( SUBMAT ),
00239      $                   WORK( IQ-1+IWORK( IQPTR+CURR ) ), MATSIZ, WORK,
00240      $                   INFO )
00241             IF( INFO.NE.0 )
00242      $         GO TO 130
00243             IF( ICOMPQ.EQ.1 ) THEN
00244                CALL SGEMM( 'N', 'N', QSIZ, MATSIZ, MATSIZ, ONE,
00245      $                     Q( 1, SUBMAT ), LDQ, WORK( IQ-1+IWORK( IQPTR+
00246      $                     CURR ) ), MATSIZ, ZERO, QSTORE( 1, SUBMAT ),
00247      $                     LDQS )
00248             END IF
00249             IWORK( IQPTR+CURR+1 ) = IWORK( IQPTR+CURR ) + MATSIZ**2
00250             CURR = CURR + 1
00251          END IF
00252          K = 1
00253          DO 60 J = SUBMAT, IWORK( I+1 )
00254             IWORK( INDXQ+J ) = K
00255             K = K + 1
00256    60    CONTINUE
00257    70 CONTINUE
00258 *
00259 *     Successively merge eigensystems of adjacent submatrices
00260 *     into eigensystem for the corresponding larger matrix.
00261 *
00262 *     while ( SUBPBS > 1 )
00263 *
00264       CURLVL = 1
00265    80 CONTINUE
00266       IF( SUBPBS.GT.1 ) THEN
00267          SPM2 = SUBPBS - 2
00268          DO 90 I = 0, SPM2, 2
00269             IF( I.EQ.0 ) THEN
00270                SUBMAT = 1
00271                MATSIZ = IWORK( 2 )
00272                MSD2 = IWORK( 1 )
00273                CURPRB = 0
00274             ELSE
00275                SUBMAT = IWORK( I ) + 1
00276                MATSIZ = IWORK( I+2 ) - IWORK( I )
00277                MSD2 = MATSIZ / 2
00278                CURPRB = CURPRB + 1
00279             END IF
00280 *
00281 *     Merge lower order eigensystems (of size MSD2 and MATSIZ - MSD2)
00282 *     into an eigensystem of size MATSIZ.
00283 *     SLAED1 is used only for the full eigensystem of a tridiagonal
00284 *     matrix.
00285 *     SLAED7 handles the cases in which eigenvalues only or eigenvalues
00286 *     and eigenvectors of a full symmetric matrix (which was reduced to
00287 *     tridiagonal form) are desired.
00288 *
00289             IF( ICOMPQ.EQ.2 ) THEN
00290                CALL SLAED1( MATSIZ, D( SUBMAT ), Q( SUBMAT, SUBMAT ),
00291      $                      LDQ, IWORK( INDXQ+SUBMAT ),
00292      $                      E( SUBMAT+MSD2-1 ), MSD2, WORK,
00293      $                      IWORK( SUBPBS+1 ), INFO )
00294             ELSE
00295                CALL SLAED7( ICOMPQ, MATSIZ, QSIZ, TLVLS, CURLVL, CURPRB,
00296      $                      D( SUBMAT ), QSTORE( 1, SUBMAT ), LDQS,
00297      $                      IWORK( INDXQ+SUBMAT ), E( SUBMAT+MSD2-1 ),
00298      $                      MSD2, WORK( IQ ), IWORK( IQPTR ),
00299      $                      IWORK( IPRMPT ), IWORK( IPERM ),
00300      $                      IWORK( IGIVPT ), IWORK( IGIVCL ),
00301      $                      WORK( IGIVNM ), WORK( IWREM ),
00302      $                      IWORK( SUBPBS+1 ), INFO )
00303             END IF
00304             IF( INFO.NE.0 )
00305      $         GO TO 130
00306             IWORK( I / 2+1 ) = IWORK( I+2 )
00307    90    CONTINUE
00308          SUBPBS = SUBPBS / 2
00309          CURLVL = CURLVL + 1
00310          GO TO 80
00311       END IF
00312 *
00313 *     end while
00314 *
00315 *     Re-merge the eigenvalues/vectors which were deflated at the final
00316 *     merge step.
00317 *
00318       IF( ICOMPQ.EQ.1 ) THEN
00319          DO 100 I = 1, N
00320             J = IWORK( INDXQ+I )
00321             WORK( I ) = D( J )
00322             CALL SCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
00323   100    CONTINUE
00324          CALL SCOPY( N, WORK, 1, D, 1 )
00325       ELSE IF( ICOMPQ.EQ.2 ) THEN
00326          DO 110 I = 1, N
00327             J = IWORK( INDXQ+I )
00328             WORK( I ) = D( J )
00329             CALL SCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
00330   110    CONTINUE
00331          CALL SCOPY( N, WORK, 1, D, 1 )
00332          CALL SLACPY( 'A', N, N, WORK( N+1 ), N, Q, LDQ )
00333       ELSE
00334          DO 120 I = 1, N
00335             J = IWORK( INDXQ+I )
00336             WORK( I ) = D( J )
00337   120    CONTINUE
00338          CALL SCOPY( N, WORK, 1, D, 1 )
00339       END IF
00340       GO TO 140
00341 *
00342   130 CONTINUE
00343       INFO = SUBMAT*( N+1 ) + SUBMAT + MATSIZ - 1
00344 *
00345   140 CONTINUE
00346       RETURN
00347 *
00348 *     End of SLAED0
00349 *
00350       END