Return to zlaed0.f CVS log

Up to [local] / rpl / lapack / lapack

Patches pour OS/2

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