Return to dlaed0.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    1: *> \brief \b DLAED0 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 DLAED0 + dependencies 
   10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaed0.f"> 
   11: *> [TGZ]</a> 
   12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaed0.f"> 
   13: *> [ZIP]</a> 
   14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaed0.f"> 
   15: *> [TXT]</a>
   16: *> \endhtmlonly 
   17: *
   18: *  Definition:
   19: *  ===========
   20: *
   21: *       SUBROUTINE DLAED0( 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: *       DOUBLE PRECISION   D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
   30: *      $                   WORK( * )
   31: *       ..
   32: *  
   33: *
   34: *> \par Purpose:
   35: *  =============
   36: *>
   37: *> \verbatim
   38: *>
   39: *> DLAED0 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DLAED0( 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:       DOUBLE PRECISION   D( * ), E( * ), Q( LDQ, * ), QSTORE( LDQS, * ),
  186:      $                   WORK( * )
  187: *     ..
  188: *
  189: *  =====================================================================
  190: *
  191: *     .. Parameters ..
  192:       DOUBLE PRECISION   ZERO, ONE, TWO
  193:       PARAMETER          ( ZERO = 0.D0, ONE = 1.D0, TWO = 2.D0 )
  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:       DOUBLE PRECISION   TEMP
  201: *     ..
  202: *     .. External Subroutines ..
  203:       EXTERNAL           DCOPY, DGEMM, DLACPY, DLAED1, DLAED7, DSTEQR,
  204:      $                   XERBLA
  205: *     ..
  206: *     .. External Functions ..
  207:       INTEGER            ILAENV
  208:       EXTERNAL           ILAENV
  209: *     ..
  210: *     .. Intrinsic Functions ..
  211:       INTRINSIC          ABS, DBLE, INT, LOG, MAX
  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( 'DLAED0', -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, 'DLAED0', ' ', 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( DBLE( 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 DSTEQR( '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 DSTEQR( '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 DGEMM( '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: *     DLAED1 is used only for the full eigensystem of a tridiagonal
  368: *     matrix.
  369: *     DLAED7 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 DLAED1( 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 DLAED7( 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 DCOPY( QSIZ, QSTORE( 1, J ), 1, Q( 1, I ), 1 )
  407:   100    CONTINUE
  408:          CALL DCOPY( 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 DCOPY( N, Q( 1, J ), 1, WORK( N*I+1 ), 1 )
  414:   110    CONTINUE
  415:          CALL DCOPY( N, WORK, 1, D, 1 )
  416:          CALL DLACPY( '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 DCOPY( 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 DLAED0
  433: *
  434:       END