Annotation of rpl/lapack/lapack/dlasd7.f, revision 1.19

1.11      bertrand    1: *> \brief \b DLASD7 merges the two sets of singular values together into a single sorted set. Then it tries to deflate the size of the problem. Used by sbdsdc.
1.8       bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
1.15      bertrand    5: * Online html documentation available at
                      6: *            http://www.netlib.org/lapack/explore-html/
1.8       bertrand    7: *
                      8: *> \htmlonly
1.15      bertrand    9: *> Download DLASD7 + dependencies
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd7.f">
                     11: *> [TGZ]</a>
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd7.f">
                     13: *> [ZIP]</a>
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd7.f">
1.8       bertrand   15: *> [TXT]</a>
1.15      bertrand   16: *> \endhtmlonly
1.8       bertrand   17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
                     22: *                          VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
                     23: *                          PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
                     24: *                          C, S, INFO )
1.15      bertrand   25: *
1.8       bertrand   26: *       .. Scalar Arguments ..
                     27: *       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
                     28: *      $                   NR, SQRE
                     29: *       DOUBLE PRECISION   ALPHA, BETA, C, S
                     30: *       ..
                     31: *       .. Array Arguments ..
                     32: *       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
                     33: *      $                   IDXQ( * ), PERM( * )
                     34: *       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
                     35: *      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
                     36: *      $                   ZW( * )
                     37: *       ..
1.15      bertrand   38: *
1.8       bertrand   39: *
                     40: *> \par Purpose:
                     41: *  =============
                     42: *>
                     43: *> \verbatim
                     44: *>
                     45: *> DLASD7 merges the two sets of singular values together into a single
                     46: *> sorted set. Then it tries to deflate the size of the problem. There
                     47: *> are two ways in which deflation can occur:  when two or more singular
                     48: *> values are close together or if there is a tiny entry in the Z
                     49: *> vector. For each such occurrence the order of the related
                     50: *> secular equation problem is reduced by one.
                     51: *>
                     52: *> DLASD7 is called from DLASD6.
                     53: *> \endverbatim
                     54: *
                     55: *  Arguments:
                     56: *  ==========
                     57: *
                     58: *> \param[in] ICOMPQ
                     59: *> \verbatim
                     60: *>          ICOMPQ is INTEGER
                     61: *>          Specifies whether singular vectors are to be computed
                     62: *>          in compact form, as follows:
                     63: *>          = 0: Compute singular values only.
                     64: *>          = 1: Compute singular vectors of upper
                     65: *>               bidiagonal matrix in compact form.
                     66: *> \endverbatim
                     67: *>
                     68: *> \param[in] NL
                     69: *> \verbatim
                     70: *>          NL is INTEGER
                     71: *>         The row dimension of the upper block. NL >= 1.
                     72: *> \endverbatim
                     73: *>
                     74: *> \param[in] NR
                     75: *> \verbatim
                     76: *>          NR is INTEGER
                     77: *>         The row dimension of the lower block. NR >= 1.
                     78: *> \endverbatim
                     79: *>
                     80: *> \param[in] SQRE
                     81: *> \verbatim
                     82: *>          SQRE is INTEGER
                     83: *>         = 0: the lower block is an NR-by-NR square matrix.
                     84: *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
                     85: *>
                     86: *>         The bidiagonal matrix has
                     87: *>         N = NL + NR + 1 rows and
                     88: *>         M = N + SQRE >= N columns.
                     89: *> \endverbatim
                     90: *>
                     91: *> \param[out] K
                     92: *> \verbatim
                     93: *>          K is INTEGER
                     94: *>         Contains the dimension of the non-deflated matrix, this is
                     95: *>         the order of the related secular equation. 1 <= K <=N.
                     96: *> \endverbatim
                     97: *>
                     98: *> \param[in,out] D
                     99: *> \verbatim
                    100: *>          D is DOUBLE PRECISION array, dimension ( N )
                    101: *>         On entry D contains the singular values of the two submatrices
                    102: *>         to be combined. On exit D contains the trailing (N-K) updated
                    103: *>         singular values (those which were deflated) sorted into
                    104: *>         increasing order.
                    105: *> \endverbatim
                    106: *>
                    107: *> \param[out] Z
                    108: *> \verbatim
                    109: *>          Z is DOUBLE PRECISION array, dimension ( M )
                    110: *>         On exit Z contains the updating row vector in the secular
                    111: *>         equation.
                    112: *> \endverbatim
                    113: *>
                    114: *> \param[out] ZW
                    115: *> \verbatim
                    116: *>          ZW is DOUBLE PRECISION array, dimension ( M )
                    117: *>         Workspace for Z.
                    118: *> \endverbatim
                    119: *>
                    120: *> \param[in,out] VF
                    121: *> \verbatim
                    122: *>          VF is DOUBLE PRECISION array, dimension ( M )
                    123: *>         On entry, VF(1:NL+1) contains the first components of all
                    124: *>         right singular vectors of the upper block; and VF(NL+2:M)
                    125: *>         contains the first components of all right singular vectors
                    126: *>         of the lower block. On exit, VF contains the first components
                    127: *>         of all right singular vectors of the bidiagonal matrix.
                    128: *> \endverbatim
                    129: *>
                    130: *> \param[out] VFW
                    131: *> \verbatim
                    132: *>          VFW is DOUBLE PRECISION array, dimension ( M )
                    133: *>         Workspace for VF.
                    134: *> \endverbatim
                    135: *>
                    136: *> \param[in,out] VL
                    137: *> \verbatim
                    138: *>          VL is DOUBLE PRECISION array, dimension ( M )
                    139: *>         On entry, VL(1:NL+1) contains the  last components of all
                    140: *>         right singular vectors of the upper block; and VL(NL+2:M)
                    141: *>         contains the last components of all right singular vectors
                    142: *>         of the lower block. On exit, VL contains the last components
                    143: *>         of all right singular vectors of the bidiagonal matrix.
                    144: *> \endverbatim
                    145: *>
                    146: *> \param[out] VLW
                    147: *> \verbatim
                    148: *>          VLW is DOUBLE PRECISION array, dimension ( M )
                    149: *>         Workspace for VL.
                    150: *> \endverbatim
                    151: *>
                    152: *> \param[in] ALPHA
                    153: *> \verbatim
                    154: *>          ALPHA is DOUBLE PRECISION
                    155: *>         Contains the diagonal element associated with the added row.
                    156: *> \endverbatim
                    157: *>
                    158: *> \param[in] BETA
                    159: *> \verbatim
                    160: *>          BETA is DOUBLE PRECISION
                    161: *>         Contains the off-diagonal element associated with the added
                    162: *>         row.
                    163: *> \endverbatim
                    164: *>
                    165: *> \param[out] DSIGMA
                    166: *> \verbatim
                    167: *>          DSIGMA is DOUBLE PRECISION array, dimension ( N )
                    168: *>         Contains a copy of the diagonal elements (K-1 singular values
                    169: *>         and one zero) in the secular equation.
                    170: *> \endverbatim
                    171: *>
                    172: *> \param[out] IDX
                    173: *> \verbatim
                    174: *>          IDX is INTEGER array, dimension ( N )
                    175: *>         This will contain the permutation used to sort the contents of
                    176: *>         D into ascending order.
                    177: *> \endverbatim
                    178: *>
                    179: *> \param[out] IDXP
                    180: *> \verbatim
                    181: *>          IDXP is INTEGER array, dimension ( N )
                    182: *>         This will contain the permutation used to place deflated
                    183: *>         values of D at the end of the array. On output IDXP(2:K)
                    184: *>         points to the nondeflated D-values and IDXP(K+1:N)
                    185: *>         points to the deflated singular values.
                    186: *> \endverbatim
                    187: *>
                    188: *> \param[in] IDXQ
                    189: *> \verbatim
                    190: *>          IDXQ is INTEGER array, dimension ( N )
                    191: *>         This contains the permutation which separately sorts the two
                    192: *>         sub-problems in D into ascending order.  Note that entries in
                    193: *>         the first half of this permutation must first be moved one
                    194: *>         position backward; and entries in the second half
                    195: *>         must first have NL+1 added to their values.
                    196: *> \endverbatim
                    197: *>
                    198: *> \param[out] PERM
                    199: *> \verbatim
                    200: *>          PERM is INTEGER array, dimension ( N )
                    201: *>         The permutations (from deflation and sorting) to be applied
                    202: *>         to each singular block. Not referenced if ICOMPQ = 0.
                    203: *> \endverbatim
                    204: *>
                    205: *> \param[out] GIVPTR
                    206: *> \verbatim
                    207: *>          GIVPTR is INTEGER
                    208: *>         The number of Givens rotations which took place in this
                    209: *>         subproblem. Not referenced if ICOMPQ = 0.
                    210: *> \endverbatim
                    211: *>
                    212: *> \param[out] GIVCOL
                    213: *> \verbatim
                    214: *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
                    215: *>         Each pair of numbers indicates a pair of columns to take place
                    216: *>         in a Givens rotation. Not referenced if ICOMPQ = 0.
                    217: *> \endverbatim
                    218: *>
                    219: *> \param[in] LDGCOL
                    220: *> \verbatim
                    221: *>          LDGCOL is INTEGER
                    222: *>         The leading dimension of GIVCOL, must be at least N.
                    223: *> \endverbatim
                    224: *>
                    225: *> \param[out] GIVNUM
                    226: *> \verbatim
                    227: *>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
                    228: *>         Each number indicates the C or S value to be used in the
                    229: *>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
                    230: *> \endverbatim
                    231: *>
                    232: *> \param[in] LDGNUM
                    233: *> \verbatim
                    234: *>          LDGNUM is INTEGER
                    235: *>         The leading dimension of GIVNUM, must be at least N.
                    236: *> \endverbatim
                    237: *>
                    238: *> \param[out] C
                    239: *> \verbatim
                    240: *>          C is DOUBLE PRECISION
                    241: *>         C contains garbage if SQRE =0 and the C-value of a Givens
                    242: *>         rotation related to the right null space if SQRE = 1.
                    243: *> \endverbatim
                    244: *>
                    245: *> \param[out] S
                    246: *> \verbatim
                    247: *>          S is DOUBLE PRECISION
                    248: *>         S contains garbage if SQRE =0 and the S-value of a Givens
                    249: *>         rotation related to the right null space if SQRE = 1.
                    250: *> \endverbatim
                    251: *>
                    252: *> \param[out] INFO
                    253: *> \verbatim
                    254: *>          INFO is INTEGER
                    255: *>         = 0:  successful exit.
                    256: *>         < 0:  if INFO = -i, the i-th argument had an illegal value.
                    257: *> \endverbatim
                    258: *
                    259: *  Authors:
                    260: *  ========
                    261: *
1.15      bertrand  262: *> \author Univ. of Tennessee
                    263: *> \author Univ. of California Berkeley
                    264: *> \author Univ. of Colorado Denver
                    265: *> \author NAG Ltd.
1.8       bertrand  266: *
1.15      bertrand  267: *> \ingroup OTHERauxiliary
1.8       bertrand  268: *
                    269: *> \par Contributors:
                    270: *  ==================
                    271: *>
                    272: *>     Ming Gu and Huan Ren, Computer Science Division, University of
                    273: *>     California at Berkeley, USA
                    274: *>
                    275: *  =====================================================================
1.1       bertrand  276:       SUBROUTINE DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, ZW, VF, VFW, VL,
                    277:      $                   VLW, ALPHA, BETA, DSIGMA, IDX, IDXP, IDXQ,
                    278:      $                   PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM,
                    279:      $                   C, S, INFO )
                    280: *
1.19    ! bertrand  281: *  -- LAPACK auxiliary routine --
1.1       bertrand  282: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    283: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
                    284: *
                    285: *     .. Scalar Arguments ..
                    286:       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
                    287:      $                   NR, SQRE
                    288:       DOUBLE PRECISION   ALPHA, BETA, C, S
                    289: *     ..
                    290: *     .. Array Arguments ..
                    291:       INTEGER            GIVCOL( LDGCOL, * ), IDX( * ), IDXP( * ),
                    292:      $                   IDXQ( * ), PERM( * )
                    293:       DOUBLE PRECISION   D( * ), DSIGMA( * ), GIVNUM( LDGNUM, * ),
                    294:      $                   VF( * ), VFW( * ), VL( * ), VLW( * ), Z( * ),
                    295:      $                   ZW( * )
                    296: *     ..
                    297: *
                    298: *  =====================================================================
                    299: *
                    300: *     .. Parameters ..
                    301:       DOUBLE PRECISION   ZERO, ONE, TWO, EIGHT
                    302:       PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, TWO = 2.0D+0,
                    303:      $                   EIGHT = 8.0D+0 )
                    304: *     ..
                    305: *     .. Local Scalars ..
                    306: *
                    307:       INTEGER            I, IDXI, IDXJ, IDXJP, J, JP, JPREV, K2, M, N,
                    308:      $                   NLP1, NLP2
                    309:       DOUBLE PRECISION   EPS, HLFTOL, TAU, TOL, Z1
                    310: *     ..
                    311: *     .. External Subroutines ..
                    312:       EXTERNAL           DCOPY, DLAMRG, DROT, XERBLA
                    313: *     ..
                    314: *     .. External Functions ..
                    315:       DOUBLE PRECISION   DLAMCH, DLAPY2
                    316:       EXTERNAL           DLAMCH, DLAPY2
                    317: *     ..
                    318: *     .. Intrinsic Functions ..
                    319:       INTRINSIC          ABS, MAX
                    320: *     ..
                    321: *     .. Executable Statements ..
                    322: *
                    323: *     Test the input parameters.
                    324: *
                    325:       INFO = 0
                    326:       N = NL + NR + 1
                    327:       M = N + SQRE
                    328: *
                    329:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
                    330:          INFO = -1
                    331:       ELSE IF( NL.LT.1 ) THEN
                    332:          INFO = -2
                    333:       ELSE IF( NR.LT.1 ) THEN
                    334:          INFO = -3
                    335:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
                    336:          INFO = -4
                    337:       ELSE IF( LDGCOL.LT.N ) THEN
                    338:          INFO = -22
                    339:       ELSE IF( LDGNUM.LT.N ) THEN
                    340:          INFO = -24
                    341:       END IF
                    342:       IF( INFO.NE.0 ) THEN
                    343:          CALL XERBLA( 'DLASD7', -INFO )
                    344:          RETURN
                    345:       END IF
                    346: *
                    347:       NLP1 = NL + 1
                    348:       NLP2 = NL + 2
                    349:       IF( ICOMPQ.EQ.1 ) THEN
                    350:          GIVPTR = 0
                    351:       END IF
                    352: *
                    353: *     Generate the first part of the vector Z and move the singular
                    354: *     values in the first part of D one position backward.
                    355: *
                    356:       Z1 = ALPHA*VL( NLP1 )
                    357:       VL( NLP1 ) = ZERO
                    358:       TAU = VF( NLP1 )
                    359:       DO 10 I = NL, 1, -1
                    360:          Z( I+1 ) = ALPHA*VL( I )
                    361:          VL( I ) = ZERO
                    362:          VF( I+1 ) = VF( I )
                    363:          D( I+1 ) = D( I )
                    364:          IDXQ( I+1 ) = IDXQ( I ) + 1
                    365:    10 CONTINUE
                    366:       VF( 1 ) = TAU
                    367: *
                    368: *     Generate the second part of the vector Z.
                    369: *
                    370:       DO 20 I = NLP2, M
                    371:          Z( I ) = BETA*VF( I )
                    372:          VF( I ) = ZERO
                    373:    20 CONTINUE
                    374: *
                    375: *     Sort the singular values into increasing order
                    376: *
                    377:       DO 30 I = NLP2, N
                    378:          IDXQ( I ) = IDXQ( I ) + NLP1
                    379:    30 CONTINUE
                    380: *
                    381: *     DSIGMA, IDXC, IDXC, and ZW are used as storage space.
                    382: *
                    383:       DO 40 I = 2, N
                    384:          DSIGMA( I ) = D( IDXQ( I ) )
                    385:          ZW( I ) = Z( IDXQ( I ) )
                    386:          VFW( I ) = VF( IDXQ( I ) )
                    387:          VLW( I ) = VL( IDXQ( I ) )
                    388:    40 CONTINUE
                    389: *
                    390:       CALL DLAMRG( NL, NR, DSIGMA( 2 ), 1, 1, IDX( 2 ) )
                    391: *
                    392:       DO 50 I = 2, N
                    393:          IDXI = 1 + IDX( I )
                    394:          D( I ) = DSIGMA( IDXI )
                    395:          Z( I ) = ZW( IDXI )
                    396:          VF( I ) = VFW( IDXI )
                    397:          VL( I ) = VLW( IDXI )
                    398:    50 CONTINUE
                    399: *
1.18      bertrand  400: *     Calculate the allowable deflation tolerance
1.1       bertrand  401: *
                    402:       EPS = DLAMCH( 'Epsilon' )
                    403:       TOL = MAX( ABS( ALPHA ), ABS( BETA ) )
                    404:       TOL = EIGHT*EIGHT*EPS*MAX( ABS( D( N ) ), TOL )
                    405: *
                    406: *     There are 2 kinds of deflation -- first a value in the z-vector
                    407: *     is small, second two (or more) singular values are very close
                    408: *     together (their difference is small).
                    409: *
                    410: *     If the value in the z-vector is small, we simply permute the
                    411: *     array so that the corresponding singular value is moved to the
                    412: *     end.
                    413: *
                    414: *     If two values in the D-vector are close, we perform a two-sided
                    415: *     rotation designed to make one of the corresponding z-vector
                    416: *     entries zero, and then permute the array so that the deflated
                    417: *     singular value is moved to the end.
                    418: *
                    419: *     If there are multiple singular values then the problem deflates.
                    420: *     Here the number of equal singular values are found.  As each equal
                    421: *     singular value is found, an elementary reflector is computed to
                    422: *     rotate the corresponding singular subspace so that the
                    423: *     corresponding components of Z are zero in this new basis.
                    424: *
                    425:       K = 1
                    426:       K2 = N + 1
                    427:       DO 60 J = 2, N
                    428:          IF( ABS( Z( J ) ).LE.TOL ) THEN
                    429: *
                    430: *           Deflate due to small z component.
                    431: *
                    432:             K2 = K2 - 1
                    433:             IDXP( K2 ) = J
                    434:             IF( J.EQ.N )
                    435:      $         GO TO 100
                    436:          ELSE
                    437:             JPREV = J
                    438:             GO TO 70
                    439:          END IF
                    440:    60 CONTINUE
                    441:    70 CONTINUE
                    442:       J = JPREV
                    443:    80 CONTINUE
                    444:       J = J + 1
                    445:       IF( J.GT.N )
                    446:      $   GO TO 90
                    447:       IF( ABS( Z( J ) ).LE.TOL ) THEN
                    448: *
                    449: *        Deflate due to small z component.
                    450: *
                    451:          K2 = K2 - 1
                    452:          IDXP( K2 ) = J
                    453:       ELSE
                    454: *
                    455: *        Check if singular values are close enough to allow deflation.
                    456: *
                    457:          IF( ABS( D( J )-D( JPREV ) ).LE.TOL ) THEN
                    458: *
                    459: *           Deflation is possible.
                    460: *
                    461:             S = Z( JPREV )
                    462:             C = Z( J )
                    463: *
                    464: *           Find sqrt(a**2+b**2) without overflow or
                    465: *           destructive underflow.
                    466: *
                    467:             TAU = DLAPY2( C, S )
                    468:             Z( J ) = TAU
                    469:             Z( JPREV ) = ZERO
                    470:             C = C / TAU
                    471:             S = -S / TAU
                    472: *
                    473: *           Record the appropriate Givens rotation
                    474: *
                    475:             IF( ICOMPQ.EQ.1 ) THEN
                    476:                GIVPTR = GIVPTR + 1
                    477:                IDXJP = IDXQ( IDX( JPREV )+1 )
                    478:                IDXJ = IDXQ( IDX( J )+1 )
                    479:                IF( IDXJP.LE.NLP1 ) THEN
                    480:                   IDXJP = IDXJP - 1
                    481:                END IF
                    482:                IF( IDXJ.LE.NLP1 ) THEN
                    483:                   IDXJ = IDXJ - 1
                    484:                END IF
                    485:                GIVCOL( GIVPTR, 2 ) = IDXJP
                    486:                GIVCOL( GIVPTR, 1 ) = IDXJ
                    487:                GIVNUM( GIVPTR, 2 ) = C
                    488:                GIVNUM( GIVPTR, 1 ) = S
                    489:             END IF
                    490:             CALL DROT( 1, VF( JPREV ), 1, VF( J ), 1, C, S )
                    491:             CALL DROT( 1, VL( JPREV ), 1, VL( J ), 1, C, S )
                    492:             K2 = K2 - 1
                    493:             IDXP( K2 ) = JPREV
                    494:             JPREV = J
                    495:          ELSE
                    496:             K = K + 1
                    497:             ZW( K ) = Z( JPREV )
                    498:             DSIGMA( K ) = D( JPREV )
                    499:             IDXP( K ) = JPREV
                    500:             JPREV = J
                    501:          END IF
                    502:       END IF
                    503:       GO TO 80
                    504:    90 CONTINUE
                    505: *
                    506: *     Record the last singular value.
                    507: *
                    508:       K = K + 1
                    509:       ZW( K ) = Z( JPREV )
                    510:       DSIGMA( K ) = D( JPREV )
                    511:       IDXP( K ) = JPREV
                    512: *
                    513:   100 CONTINUE
                    514: *
                    515: *     Sort the singular values into DSIGMA. The singular values which
                    516: *     were not deflated go into the first K slots of DSIGMA, except
                    517: *     that DSIGMA(1) is treated separately.
                    518: *
                    519:       DO 110 J = 2, N
                    520:          JP = IDXP( J )
                    521:          DSIGMA( J ) = D( JP )
                    522:          VFW( J ) = VF( JP )
                    523:          VLW( J ) = VL( JP )
                    524:   110 CONTINUE
                    525:       IF( ICOMPQ.EQ.1 ) THEN
                    526:          DO 120 J = 2, N
                    527:             JP = IDXP( J )
                    528:             PERM( J ) = IDXQ( IDX( JP )+1 )
                    529:             IF( PERM( J ).LE.NLP1 ) THEN
                    530:                PERM( J ) = PERM( J ) - 1
                    531:             END IF
                    532:   120    CONTINUE
                    533:       END IF
                    534: *
                    535: *     The deflated singular values go back into the last N - K slots of
                    536: *     D.
                    537: *
                    538:       CALL DCOPY( N-K, DSIGMA( K+1 ), 1, D( K+1 ), 1 )
                    539: *
                    540: *     Determine DSIGMA(1), DSIGMA(2), Z(1), VF(1), VL(1), VF(M), and
                    541: *     VL(M).
                    542: *
                    543:       DSIGMA( 1 ) = ZERO
                    544:       HLFTOL = TOL / TWO
                    545:       IF( ABS( DSIGMA( 2 ) ).LE.HLFTOL )
                    546:      $   DSIGMA( 2 ) = HLFTOL
                    547:       IF( M.GT.N ) THEN
                    548:          Z( 1 ) = DLAPY2( Z1, Z( M ) )
                    549:          IF( Z( 1 ).LE.TOL ) THEN
                    550:             C = ONE
                    551:             S = ZERO
                    552:             Z( 1 ) = TOL
                    553:          ELSE
                    554:             C = Z1 / Z( 1 )
                    555:             S = -Z( M ) / Z( 1 )
                    556:          END IF
                    557:          CALL DROT( 1, VF( M ), 1, VF( 1 ), 1, C, S )
                    558:          CALL DROT( 1, VL( M ), 1, VL( 1 ), 1, C, S )
                    559:       ELSE
                    560:          IF( ABS( Z1 ).LE.TOL ) THEN
                    561:             Z( 1 ) = TOL
                    562:          ELSE
                    563:             Z( 1 ) = Z1
                    564:          END IF
                    565:       END IF
                    566: *
                    567: *     Restore Z, VF, and VL.
                    568: *
                    569:       CALL DCOPY( K-1, ZW( 2 ), 1, Z( 2 ), 1 )
                    570:       CALL DCOPY( N-1, VFW( 2 ), 1, VF( 2 ), 1 )
                    571:       CALL DCOPY( N-1, VLW( 2 ), 1, VL( 2 ), 1 )
                    572: *
                    573:       RETURN
                    574: *
                    575: *     End of DLASD7
                    576: *
                    577:       END

CVSweb interface <joel.bertrand@systella.fr>