Annotation of rpl/lapack/lapack/dlaed0.f, revision 1.11

1.11    ! bertrand    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.
1.8       bertrand    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: *
1.11    ! bertrand  161: *> \date September 2012
1.8       bertrand  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: *  =====================================================================
1.1       bertrand  172:       SUBROUTINE DLAED0( ICOMPQ, QSIZ, N, D, E, Q, LDQ, QSTORE, LDQS,
                    173:      $                   WORK, IWORK, INFO )
                    174: *
1.11    ! bertrand  175: *  -- LAPACK computational routine (version 3.4.2) --
1.1       bertrand  176: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    177: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.11    ! bertrand  178: *     September 2012
1.1       bertrand  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

CVSweb interface <joel.bertrand@systella.fr>