Annotation of rpl/lapack/lapack/dlasd6.f, revision 1.17

1.14      bertrand    1: *> \brief \b DLASD6 computes the SVD of an updated upper bidiagonal matrix obtained by merging two smaller ones by appending a row. Used by sbdsdc.
1.11      bertrand    2: *
                      3: *  =========== DOCUMENTATION ===========
                      4: *
                      5: * Online html documentation available at 
                      6: *            http://www.netlib.org/lapack/explore-html/ 
                      7: *
                      8: *> \htmlonly
                      9: *> Download DLASD6 + dependencies 
                     10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd6.f"> 
                     11: *> [TGZ]</a> 
                     12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd6.f"> 
                     13: *> [ZIP]</a> 
                     14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd6.f"> 
                     15: *> [TXT]</a>
                     16: *> \endhtmlonly 
                     17: *
                     18: *  Definition:
                     19: *  ===========
                     20: *
                     21: *       SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
                     22: *                          IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
                     23: *                          LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
                     24: *                          IWORK, INFO )
                     25: * 
                     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, * ), IDXQ( * ), IWORK( * ),
                     33: *      $                   PERM( * )
                     34: *       DOUBLE PRECISION   D( * ), DIFL( * ), DIFR( * ),
                     35: *      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
                     36: *      $                   VF( * ), VL( * ), WORK( * ), Z( * )
                     37: *       ..
                     38: *  
                     39: *
                     40: *> \par Purpose:
                     41: *  =============
                     42: *>
                     43: *> \verbatim
                     44: *>
                     45: *> DLASD6 computes the SVD of an updated upper bidiagonal matrix B
                     46: *> obtained by merging two smaller ones by appending a row. This
                     47: *> routine is used only for the problem which requires all singular
                     48: *> values and optionally singular vector matrices in factored form.
                     49: *> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE.
                     50: *> A related subroutine, DLASD1, handles the case in which all singular
                     51: *> values and singular vectors of the bidiagonal matrix are desired.
                     52: *>
                     53: *> DLASD6 computes the SVD as follows:
                     54: *>
                     55: *>               ( D1(in)    0    0       0 )
                     56: *>   B = U(in) * (   Z1**T   a   Z2**T    b ) * VT(in)
                     57: *>               (   0       0   D2(in)   0 )
                     58: *>
                     59: *>     = U(out) * ( D(out) 0) * VT(out)
                     60: *>
                     61: *> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M
                     62: *> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros
                     63: *> elsewhere; and the entry b is empty if SQRE = 0.
                     64: *>
                     65: *> The singular values of B can be computed using D1, D2, the first
                     66: *> components of all the right singular vectors of the lower block, and
                     67: *> the last components of all the right singular vectors of the upper
                     68: *> block. These components are stored and updated in VF and VL,
                     69: *> respectively, in DLASD6. Hence U and VT are not explicitly
                     70: *> referenced.
                     71: *>
                     72: *> The singular values are stored in D. The algorithm consists of two
                     73: *> stages:
                     74: *>
                     75: *>       The first stage consists of deflating the size of the problem
                     76: *>       when there are multiple singular values or if there is a zero
                     77: *>       in the Z vector. For each such occurence the dimension of the
                     78: *>       secular equation problem is reduced by one. This stage is
                     79: *>       performed by the routine DLASD7.
                     80: *>
                     81: *>       The second stage consists of calculating the updated
                     82: *>       singular values. This is done by finding the roots of the
                     83: *>       secular equation via the routine DLASD4 (as called by DLASD8).
                     84: *>       This routine also updates VF and VL and computes the distances
                     85: *>       between the updated singular values and the old singular
                     86: *>       values.
                     87: *>
                     88: *> DLASD6 is called from DLASDA.
                     89: *> \endverbatim
                     90: *
                     91: *  Arguments:
                     92: *  ==========
                     93: *
                     94: *> \param[in] ICOMPQ
                     95: *> \verbatim
                     96: *>          ICOMPQ is INTEGER
                     97: *>         Specifies whether singular vectors are to be computed in
                     98: *>         factored form:
                     99: *>         = 0: Compute singular values only.
                    100: *>         = 1: Compute singular vectors in factored form as well.
                    101: *> \endverbatim
                    102: *>
                    103: *> \param[in] NL
                    104: *> \verbatim
                    105: *>          NL is INTEGER
                    106: *>         The row dimension of the upper block.  NL >= 1.
                    107: *> \endverbatim
                    108: *>
                    109: *> \param[in] NR
                    110: *> \verbatim
                    111: *>          NR is INTEGER
                    112: *>         The row dimension of the lower block.  NR >= 1.
                    113: *> \endverbatim
                    114: *>
                    115: *> \param[in] SQRE
                    116: *> \verbatim
                    117: *>          SQRE is INTEGER
                    118: *>         = 0: the lower block is an NR-by-NR square matrix.
                    119: *>         = 1: the lower block is an NR-by-(NR+1) rectangular matrix.
                    120: *>
                    121: *>         The bidiagonal matrix has row dimension N = NL + NR + 1,
                    122: *>         and column dimension M = N + SQRE.
                    123: *> \endverbatim
                    124: *>
                    125: *> \param[in,out] D
                    126: *> \verbatim
                    127: *>          D is DOUBLE PRECISION array, dimension ( NL+NR+1 ).
                    128: *>         On entry D(1:NL,1:NL) contains the singular values of the
                    129: *>         upper block, and D(NL+2:N) contains the singular values
                    130: *>         of the lower block. On exit D(1:N) contains the singular
                    131: *>         values of the modified matrix.
                    132: *> \endverbatim
                    133: *>
                    134: *> \param[in,out] VF
                    135: *> \verbatim
                    136: *>          VF is DOUBLE PRECISION array, dimension ( M )
                    137: *>         On entry, VF(1:NL+1) contains the first components of all
                    138: *>         right singular vectors of the upper block; and VF(NL+2:M)
                    139: *>         contains the first components of all right singular vectors
                    140: *>         of the lower block. On exit, VF contains the first components
                    141: *>         of all right singular vectors of the bidiagonal matrix.
                    142: *> \endverbatim
                    143: *>
                    144: *> \param[in,out] VL
                    145: *> \verbatim
                    146: *>          VL is DOUBLE PRECISION array, dimension ( M )
                    147: *>         On entry, VL(1:NL+1) contains the  last components of all
                    148: *>         right singular vectors of the upper block; and VL(NL+2:M)
                    149: *>         contains the last components of all right singular vectors of
                    150: *>         the lower block. On exit, VL contains the last components of
                    151: *>         all right singular vectors of the bidiagonal matrix.
                    152: *> \endverbatim
                    153: *>
                    154: *> \param[in,out] ALPHA
                    155: *> \verbatim
                    156: *>          ALPHA is DOUBLE PRECISION
                    157: *>         Contains the diagonal element associated with the added row.
                    158: *> \endverbatim
                    159: *>
                    160: *> \param[in,out] BETA
                    161: *> \verbatim
                    162: *>          BETA is DOUBLE PRECISION
                    163: *>         Contains the off-diagonal element associated with the added
                    164: *>         row.
                    165: *> \endverbatim
                    166: *>
1.17    ! bertrand  167: *> \param[in,out] IDXQ
1.11      bertrand  168: *> \verbatim
                    169: *>          IDXQ is INTEGER array, dimension ( N )
                    170: *>         This contains the permutation which will reintegrate the
                    171: *>         subproblem just solved back into sorted order, i.e.
                    172: *>         D( IDXQ( I = 1, N ) ) will be in ascending order.
                    173: *> \endverbatim
                    174: *>
                    175: *> \param[out] PERM
                    176: *> \verbatim
                    177: *>          PERM is INTEGER array, dimension ( N )
                    178: *>         The permutations (from deflation and sorting) to be applied
                    179: *>         to each block. Not referenced if ICOMPQ = 0.
                    180: *> \endverbatim
                    181: *>
                    182: *> \param[out] GIVPTR
                    183: *> \verbatim
                    184: *>          GIVPTR is INTEGER
                    185: *>         The number of Givens rotations which took place in this
                    186: *>         subproblem. Not referenced if ICOMPQ = 0.
                    187: *> \endverbatim
                    188: *>
                    189: *> \param[out] GIVCOL
                    190: *> \verbatim
                    191: *>          GIVCOL is INTEGER array, dimension ( LDGCOL, 2 )
                    192: *>         Each pair of numbers indicates a pair of columns to take place
                    193: *>         in a Givens rotation. Not referenced if ICOMPQ = 0.
                    194: *> \endverbatim
                    195: *>
                    196: *> \param[in] LDGCOL
                    197: *> \verbatim
                    198: *>          LDGCOL is INTEGER
                    199: *>         leading dimension of GIVCOL, must be at least N.
                    200: *> \endverbatim
                    201: *>
                    202: *> \param[out] GIVNUM
                    203: *> \verbatim
                    204: *>          GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
                    205: *>         Each number indicates the C or S value to be used in the
                    206: *>         corresponding Givens rotation. Not referenced if ICOMPQ = 0.
                    207: *> \endverbatim
                    208: *>
                    209: *> \param[in] LDGNUM
                    210: *> \verbatim
                    211: *>          LDGNUM is INTEGER
                    212: *>         The leading dimension of GIVNUM and POLES, must be at least N.
                    213: *> \endverbatim
                    214: *>
                    215: *> \param[out] POLES
                    216: *> \verbatim
                    217: *>          POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 )
                    218: *>         On exit, POLES(1,*) is an array containing the new singular
                    219: *>         values obtained from solving the secular equation, and
                    220: *>         POLES(2,*) is an array containing the poles in the secular
                    221: *>         equation. Not referenced if ICOMPQ = 0.
                    222: *> \endverbatim
                    223: *>
                    224: *> \param[out] DIFL
                    225: *> \verbatim
                    226: *>          DIFL is DOUBLE PRECISION array, dimension ( N )
                    227: *>         On exit, DIFL(I) is the distance between I-th updated
                    228: *>         (undeflated) singular value and the I-th (undeflated) old
                    229: *>         singular value.
                    230: *> \endverbatim
                    231: *>
                    232: *> \param[out] DIFR
                    233: *> \verbatim
                    234: *>          DIFR is DOUBLE PRECISION array,
                    235: *>                  dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and
                    236: *>                  dimension ( N ) if ICOMPQ = 0.
                    237: *>         On exit, DIFR(I, 1) is the distance between I-th updated
                    238: *>         (undeflated) singular value and the I+1-th (undeflated) old
                    239: *>         singular value.
                    240: *>
                    241: *>         If ICOMPQ = 1, DIFR(1:K,2) is an array containing the
                    242: *>         normalizing factors for the right singular vector matrix.
                    243: *>
                    244: *>         See DLASD8 for details on DIFL and DIFR.
                    245: *> \endverbatim
                    246: *>
                    247: *> \param[out] Z
                    248: *> \verbatim
                    249: *>          Z is DOUBLE PRECISION array, dimension ( M )
                    250: *>         The first elements of this array contain the components
                    251: *>         of the deflation-adjusted updating row vector.
                    252: *> \endverbatim
                    253: *>
                    254: *> \param[out] K
                    255: *> \verbatim
                    256: *>          K is INTEGER
                    257: *>         Contains the dimension of the non-deflated matrix,
                    258: *>         This is the order of the related secular equation. 1 <= K <=N.
                    259: *> \endverbatim
                    260: *>
                    261: *> \param[out] C
                    262: *> \verbatim
                    263: *>          C is DOUBLE PRECISION
                    264: *>         C contains garbage if SQRE =0 and the C-value of a Givens
                    265: *>         rotation related to the right null space if SQRE = 1.
                    266: *> \endverbatim
                    267: *>
                    268: *> \param[out] S
                    269: *> \verbatim
                    270: *>          S is DOUBLE PRECISION
                    271: *>         S contains garbage if SQRE =0 and the S-value of a Givens
                    272: *>         rotation related to the right null space if SQRE = 1.
                    273: *> \endverbatim
                    274: *>
                    275: *> \param[out] WORK
                    276: *> \verbatim
                    277: *>          WORK is DOUBLE PRECISION array, dimension ( 4 * M )
                    278: *> \endverbatim
                    279: *>
                    280: *> \param[out] IWORK
                    281: *> \verbatim
                    282: *>          IWORK is INTEGER array, dimension ( 3 * N )
                    283: *> \endverbatim
                    284: *>
                    285: *> \param[out] INFO
                    286: *> \verbatim
                    287: *>          INFO is INTEGER
                    288: *>          = 0:  successful exit.
                    289: *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
                    290: *>          > 0:  if INFO = 1, a singular value did not converge
                    291: *> \endverbatim
                    292: *
                    293: *  Authors:
                    294: *  ========
                    295: *
                    296: *> \author Univ. of Tennessee 
                    297: *> \author Univ. of California Berkeley 
                    298: *> \author Univ. of Colorado Denver 
                    299: *> \author NAG Ltd. 
                    300: *
1.17    ! bertrand  301: *> \date November 2015
1.11      bertrand  302: *
                    303: *> \ingroup auxOTHERauxiliary
                    304: *
                    305: *> \par Contributors:
                    306: *  ==================
                    307: *>
                    308: *>     Ming Gu and Huan Ren, Computer Science Division, University of
                    309: *>     California at Berkeley, USA
                    310: *>
                    311: *  =====================================================================
1.1       bertrand  312:       SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA,
                    313:      $                   IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM,
                    314:      $                   LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK,
                    315:      $                   IWORK, INFO )
                    316: *
1.17    ! bertrand  317: *  -- LAPACK auxiliary routine (version 3.6.0) --
1.1       bertrand  318: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
                    319: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.17    ! bertrand  320: *     November 2015
1.1       bertrand  321: *
                    322: *     .. Scalar Arguments ..
                    323:       INTEGER            GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL,
                    324:      $                   NR, SQRE
                    325:       DOUBLE PRECISION   ALPHA, BETA, C, S
                    326: *     ..
                    327: *     .. Array Arguments ..
                    328:       INTEGER            GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ),
                    329:      $                   PERM( * )
                    330:       DOUBLE PRECISION   D( * ), DIFL( * ), DIFR( * ),
                    331:      $                   GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ),
                    332:      $                   VF( * ), VL( * ), WORK( * ), Z( * )
                    333: *     ..
                    334: *
                    335: *  =====================================================================
                    336: *
                    337: *     .. Parameters ..
                    338:       DOUBLE PRECISION   ONE, ZERO
                    339:       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
                    340: *     ..
                    341: *     .. Local Scalars ..
                    342:       INTEGER            I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M,
                    343:      $                   N, N1, N2
                    344:       DOUBLE PRECISION   ORGNRM
                    345: *     ..
                    346: *     .. External Subroutines ..
                    347:       EXTERNAL           DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA
                    348: *     ..
                    349: *     .. Intrinsic Functions ..
                    350:       INTRINSIC          ABS, MAX
                    351: *     ..
                    352: *     .. Executable Statements ..
                    353: *
                    354: *     Test the input parameters.
                    355: *
                    356:       INFO = 0
                    357:       N = NL + NR + 1
                    358:       M = N + SQRE
                    359: *
                    360:       IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN
                    361:          INFO = -1
                    362:       ELSE IF( NL.LT.1 ) THEN
                    363:          INFO = -2
                    364:       ELSE IF( NR.LT.1 ) THEN
                    365:          INFO = -3
                    366:       ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN
                    367:          INFO = -4
                    368:       ELSE IF( LDGCOL.LT.N ) THEN
                    369:          INFO = -14
                    370:       ELSE IF( LDGNUM.LT.N ) THEN
                    371:          INFO = -16
                    372:       END IF
                    373:       IF( INFO.NE.0 ) THEN
                    374:          CALL XERBLA( 'DLASD6', -INFO )
                    375:          RETURN
                    376:       END IF
                    377: *
                    378: *     The following values are for bookkeeping purposes only.  They are
                    379: *     integer pointers which indicate the portion of the workspace
                    380: *     used by a particular array in DLASD7 and DLASD8.
                    381: *
                    382:       ISIGMA = 1
                    383:       IW = ISIGMA + N
                    384:       IVFW = IW + M
                    385:       IVLW = IVFW + M
                    386: *
                    387:       IDX = 1
                    388:       IDXC = IDX + N
                    389:       IDXP = IDXC + N
                    390: *
                    391: *     Scale.
                    392: *
                    393:       ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) )
                    394:       D( NL+1 ) = ZERO
                    395:       DO 10 I = 1, N
                    396:          IF( ABS( D( I ) ).GT.ORGNRM ) THEN
                    397:             ORGNRM = ABS( D( I ) )
                    398:          END IF
                    399:    10 CONTINUE
                    400:       CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
                    401:       ALPHA = ALPHA / ORGNRM
                    402:       BETA = BETA / ORGNRM
                    403: *
                    404: *     Sort and Deflate singular values.
                    405: *
                    406:       CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF,
                    407:      $             WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA,
                    408:      $             WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ,
                    409:      $             PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S,
                    410:      $             INFO )
                    411: *
                    412: *     Solve Secular Equation, compute DIFL, DIFR, and update VF, VL.
                    413: *
                    414:       CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM,
                    415:      $             WORK( ISIGMA ), WORK( IW ), INFO )
                    416: *
1.17    ! bertrand  417: *     Report the possible convergence failure.
1.8       bertrand  418: *
                    419:       IF( INFO.NE.0 ) THEN
                    420:          RETURN
                    421:       END IF
                    422: *
1.1       bertrand  423: *     Save the poles if ICOMPQ = 1.
                    424: *
                    425:       IF( ICOMPQ.EQ.1 ) THEN
                    426:          CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 )
                    427:          CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 )
                    428:       END IF
                    429: *
                    430: *     Unscale.
                    431: *
                    432:       CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
                    433: *
                    434: *     Prepare the IDXQ sorting permutation.
                    435: *
                    436:       N1 = K
                    437:       N2 = N - K
                    438:       CALL DLAMRG( N1, N2, D, 1, -1, IDXQ )
                    439: *
                    440:       RETURN
                    441: *
                    442: *     End of DLASD6
                    443: *
                    444:       END

CVSweb interface <joel.bertrand@systella.fr>