version 1.1.1.1, 2010/01/26 15:22:45
|
version 1.20, 2017/06/17 10:53:57
|
Line 1
|
Line 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. |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
*> \htmlonly |
|
*> Download DLASD6 + dependencies |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlasd6.f"> |
|
*> [TGZ]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlasd6.f"> |
|
*> [ZIP]</a> |
|
*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlasd6.f"> |
|
*> [TXT]</a> |
|
*> \endhtmlonly |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, |
|
* IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, |
|
* LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, |
|
* IWORK, INFO ) |
|
* |
|
* .. Scalar Arguments .. |
|
* INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, |
|
* $ NR, SQRE |
|
* DOUBLE PRECISION ALPHA, BETA, C, S |
|
* .. |
|
* .. Array Arguments .. |
|
* INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), |
|
* $ PERM( * ) |
|
* DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ), |
|
* $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), |
|
* $ VF( * ), VL( * ), WORK( * ), Z( * ) |
|
* .. |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DLASD6 computes the SVD of an updated upper bidiagonal matrix B |
|
*> obtained by merging two smaller ones by appending a row. This |
|
*> routine is used only for the problem which requires all singular |
|
*> values and optionally singular vector matrices in factored form. |
|
*> B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. |
|
*> A related subroutine, DLASD1, handles the case in which all singular |
|
*> values and singular vectors of the bidiagonal matrix are desired. |
|
*> |
|
*> DLASD6 computes the SVD as follows: |
|
*> |
|
*> ( D1(in) 0 0 0 ) |
|
*> B = U(in) * ( Z1**T a Z2**T b ) * VT(in) |
|
*> ( 0 0 D2(in) 0 ) |
|
*> |
|
*> = U(out) * ( D(out) 0) * VT(out) |
|
*> |
|
*> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M |
|
*> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros |
|
*> elsewhere; and the entry b is empty if SQRE = 0. |
|
*> |
|
*> The singular values of B can be computed using D1, D2, the first |
|
*> components of all the right singular vectors of the lower block, and |
|
*> the last components of all the right singular vectors of the upper |
|
*> block. These components are stored and updated in VF and VL, |
|
*> respectively, in DLASD6. Hence U and VT are not explicitly |
|
*> referenced. |
|
*> |
|
*> The singular values are stored in D. The algorithm consists of two |
|
*> stages: |
|
*> |
|
*> The first stage consists of deflating the size of the problem |
|
*> when there are multiple singular values or if there is a zero |
|
*> in the Z vector. For each such occurrence the dimension of the |
|
*> secular equation problem is reduced by one. This stage is |
|
*> performed by the routine DLASD7. |
|
*> |
|
*> The second stage consists of calculating the updated |
|
*> singular values. This is done by finding the roots of the |
|
*> secular equation via the routine DLASD4 (as called by DLASD8). |
|
*> This routine also updates VF and VL and computes the distances |
|
*> between the updated singular values and the old singular |
|
*> values. |
|
*> |
|
*> DLASD6 is called from DLASDA. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] ICOMPQ |
|
*> \verbatim |
|
*> ICOMPQ is INTEGER |
|
*> Specifies whether singular vectors are to be computed in |
|
*> factored form: |
|
*> = 0: Compute singular values only. |
|
*> = 1: Compute singular vectors in factored form as well. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] NL |
|
*> \verbatim |
|
*> NL is INTEGER |
|
*> The row dimension of the upper block. NL >= 1. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] NR |
|
*> \verbatim |
|
*> NR is INTEGER |
|
*> The row dimension of the lower block. NR >= 1. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] SQRE |
|
*> \verbatim |
|
*> SQRE is INTEGER |
|
*> = 0: the lower block is an NR-by-NR square matrix. |
|
*> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. |
|
*> |
|
*> The bidiagonal matrix has row dimension N = NL + NR + 1, |
|
*> and column dimension M = N + SQRE. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] D |
|
*> \verbatim |
|
*> D is DOUBLE PRECISION array, dimension ( NL+NR+1 ). |
|
*> On entry D(1:NL,1:NL) contains the singular values of the |
|
*> upper block, and D(NL+2:N) contains the singular values |
|
*> of the lower block. On exit D(1:N) contains the singular |
|
*> values of the modified matrix. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] VF |
|
*> \verbatim |
|
*> VF is DOUBLE PRECISION array, dimension ( M ) |
|
*> On entry, VF(1:NL+1) contains the first components of all |
|
*> right singular vectors of the upper block; and VF(NL+2:M) |
|
*> contains the first components of all right singular vectors |
|
*> of the lower block. On exit, VF contains the first components |
|
*> of all right singular vectors of the bidiagonal matrix. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] VL |
|
*> \verbatim |
|
*> VL is DOUBLE PRECISION array, dimension ( M ) |
|
*> On entry, VL(1:NL+1) contains the last components of all |
|
*> right singular vectors of the upper block; and VL(NL+2:M) |
|
*> contains the last components of all right singular vectors of |
|
*> the lower block. On exit, VL contains the last components of |
|
*> all right singular vectors of the bidiagonal matrix. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] ALPHA |
|
*> \verbatim |
|
*> ALPHA is DOUBLE PRECISION |
|
*> Contains the diagonal element associated with the added row. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] BETA |
|
*> \verbatim |
|
*> BETA is DOUBLE PRECISION |
|
*> Contains the off-diagonal element associated with the added |
|
*> row. |
|
*> \endverbatim |
|
*> |
|
*> \param[in,out] IDXQ |
|
*> \verbatim |
|
*> IDXQ is INTEGER array, dimension ( N ) |
|
*> This contains the permutation which will reintegrate the |
|
*> subproblem just solved back into sorted order, i.e. |
|
*> D( IDXQ( I = 1, N ) ) will be in ascending order. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] PERM |
|
*> \verbatim |
|
*> PERM is INTEGER array, dimension ( N ) |
|
*> The permutations (from deflation and sorting) to be applied |
|
*> to each block. Not referenced if ICOMPQ = 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] GIVPTR |
|
*> \verbatim |
|
*> GIVPTR is INTEGER |
|
*> The number of Givens rotations which took place in this |
|
*> subproblem. Not referenced if ICOMPQ = 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] GIVCOL |
|
*> \verbatim |
|
*> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 ) |
|
*> Each pair of numbers indicates a pair of columns to take place |
|
*> in a Givens rotation. Not referenced if ICOMPQ = 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDGCOL |
|
*> \verbatim |
|
*> LDGCOL is INTEGER |
|
*> leading dimension of GIVCOL, must be at least N. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] GIVNUM |
|
*> \verbatim |
|
*> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) |
|
*> Each number indicates the C or S value to be used in the |
|
*> corresponding Givens rotation. Not referenced if ICOMPQ = 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[in] LDGNUM |
|
*> \verbatim |
|
*> LDGNUM is INTEGER |
|
*> The leading dimension of GIVNUM and POLES, must be at least N. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] POLES |
|
*> \verbatim |
|
*> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) |
|
*> On exit, POLES(1,*) is an array containing the new singular |
|
*> values obtained from solving the secular equation, and |
|
*> POLES(2,*) is an array containing the poles in the secular |
|
*> equation. Not referenced if ICOMPQ = 0. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] DIFL |
|
*> \verbatim |
|
*> DIFL is DOUBLE PRECISION array, dimension ( N ) |
|
*> On exit, DIFL(I) is the distance between I-th updated |
|
*> (undeflated) singular value and the I-th (undeflated) old |
|
*> singular value. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] DIFR |
|
*> \verbatim |
|
*> DIFR is DOUBLE PRECISION array, |
|
*> dimension ( LDDIFR, 2 ) if ICOMPQ = 1 and |
|
*> dimension ( K ) if ICOMPQ = 0. |
|
*> On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not |
|
*> defined and will not be referenced. |
|
*> |
|
*> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the |
|
*> normalizing factors for the right singular vector matrix. |
|
*> |
|
*> See DLASD8 for details on DIFL and DIFR. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] Z |
|
*> \verbatim |
|
*> Z is DOUBLE PRECISION array, dimension ( M ) |
|
*> The first elements of this array contain the components |
|
*> of the deflation-adjusted updating row vector. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] K |
|
*> \verbatim |
|
*> K is INTEGER |
|
*> Contains the dimension of the non-deflated matrix, |
|
*> This is the order of the related secular equation. 1 <= K <=N. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] C |
|
*> \verbatim |
|
*> C is DOUBLE PRECISION |
|
*> C contains garbage if SQRE =0 and the C-value of a Givens |
|
*> rotation related to the right null space if SQRE = 1. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] S |
|
*> \verbatim |
|
*> S is DOUBLE PRECISION |
|
*> S contains garbage if SQRE =0 and the S-value of a Givens |
|
*> rotation related to the right null space if SQRE = 1. |
|
*> \endverbatim |
|
*> |
|
*> \param[out] WORK |
|
*> \verbatim |
|
*> WORK is DOUBLE PRECISION array, dimension ( 4 * M ) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] IWORK |
|
*> \verbatim |
|
*> IWORK is INTEGER array, dimension ( 3 * N ) |
|
*> \endverbatim |
|
*> |
|
*> \param[out] INFO |
|
*> \verbatim |
|
*> INFO is INTEGER |
|
*> = 0: successful exit. |
|
*> < 0: if INFO = -i, the i-th argument had an illegal value. |
|
*> > 0: if INFO = 1, a singular value did not converge |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date June 2016 |
|
* |
|
*> \ingroup OTHERauxiliary |
|
* |
|
*> \par Contributors: |
|
* ================== |
|
*> |
|
*> Ming Gu and Huan Ren, Computer Science Division, University of |
|
*> California at Berkeley, USA |
|
*> |
|
* ===================================================================== |
SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, |
SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, |
$ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, |
$ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, |
$ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, |
$ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, |
$ IWORK, INFO ) |
$ IWORK, INFO ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.2) -- |
* -- LAPACK auxiliary routine (version 3.7.0) -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
* November 2006 |
* June 2016 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, |
INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, |
Line 21
|
Line 331
|
$ VF( * ), VL( * ), WORK( * ), Z( * ) |
$ VF( * ), VL( * ), WORK( * ), Z( * ) |
* .. |
* .. |
* |
* |
* Purpose |
|
* ======= |
|
* |
|
* DLASD6 computes the SVD of an updated upper bidiagonal matrix B |
|
* obtained by merging two smaller ones by appending a row. This |
|
* routine is used only for the problem which requires all singular |
|
* values and optionally singular vector matrices in factored form. |
|
* B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. |
|
* A related subroutine, DLASD1, handles the case in which all singular |
|
* values and singular vectors of the bidiagonal matrix are desired. |
|
* |
|
* DLASD6 computes the SVD as follows: |
|
* |
|
* ( D1(in) 0 0 0 ) |
|
* B = U(in) * ( Z1' a Z2' b ) * VT(in) |
|
* ( 0 0 D2(in) 0 ) |
|
* |
|
* = U(out) * ( D(out) 0) * VT(out) |
|
* |
|
* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M |
|
* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros |
|
* elsewhere; and the entry b is empty if SQRE = 0. |
|
* |
|
* The singular values of B can be computed using D1, D2, the first |
|
* components of all the right singular vectors of the lower block, and |
|
* the last components of all the right singular vectors of the upper |
|
* block. These components are stored and updated in VF and VL, |
|
* respectively, in DLASD6. Hence U and VT are not explicitly |
|
* referenced. |
|
* |
|
* The singular values are stored in D. The algorithm consists of two |
|
* stages: |
|
* |
|
* The first stage consists of deflating the size of the problem |
|
* when there are multiple singular values or if there is a zero |
|
* in the Z vector. For each such occurence the dimension of the |
|
* secular equation problem is reduced by one. This stage is |
|
* performed by the routine DLASD7. |
|
* |
|
* The second stage consists of calculating the updated |
|
* singular values. This is done by finding the roots of the |
|
* secular equation via the routine DLASD4 (as called by DLASD8). |
|
* This routine also updates VF and VL and computes the distances |
|
* between the updated singular values and the old singular |
|
* values. |
|
* |
|
* DLASD6 is called from DLASDA. |
|
* |
|
* Arguments |
|
* ========= |
|
* |
|
* ICOMPQ (input) INTEGER |
|
* Specifies whether singular vectors are to be computed in |
|
* factored form: |
|
* = 0: Compute singular values only. |
|
* = 1: Compute singular vectors in factored form as well. |
|
* |
|
* NL (input) INTEGER |
|
* The row dimension of the upper block. NL >= 1. |
|
* |
|
* NR (input) INTEGER |
|
* The row dimension of the lower block. NR >= 1. |
|
* |
|
* SQRE (input) INTEGER |
|
* = 0: the lower block is an NR-by-NR square matrix. |
|
* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. |
|
* |
|
* The bidiagonal matrix has row dimension N = NL + NR + 1, |
|
* and column dimension M = N + SQRE. |
|
* |
|
* D (input/output) DOUBLE PRECISION array, dimension ( NL+NR+1 ). |
|
* On entry D(1:NL,1:NL) contains the singular values of the |
|
* upper block, and D(NL+2:N) contains the singular values |
|
* of the lower block. On exit D(1:N) contains the singular |
|
* values of the modified matrix. |
|
* |
|
* VF (input/output) DOUBLE PRECISION array, dimension ( M ) |
|
* On entry, VF(1:NL+1) contains the first components of all |
|
* right singular vectors of the upper block; and VF(NL+2:M) |
|
* contains the first components of all right singular vectors |
|
* of the lower block. On exit, VF contains the first components |
|
* of all right singular vectors of the bidiagonal matrix. |
|
* |
|
* VL (input/output) DOUBLE PRECISION array, dimension ( M ) |
|
* On entry, VL(1:NL+1) contains the last components of all |
|
* right singular vectors of the upper block; and VL(NL+2:M) |
|
* contains the last components of all right singular vectors of |
|
* the lower block. On exit, VL contains the last components of |
|
* all right singular vectors of the bidiagonal matrix. |
|
* |
|
* ALPHA (input/output) DOUBLE PRECISION |
|
* Contains the diagonal element associated with the added row. |
|
* |
|
* BETA (input/output) DOUBLE PRECISION |
|
* Contains the off-diagonal element associated with the added |
|
* row. |
|
* |
|
* IDXQ (output) INTEGER array, dimension ( N ) |
|
* This contains the permutation which will reintegrate the |
|
* subproblem just solved back into sorted order, i.e. |
|
* D( IDXQ( I = 1, N ) ) will be in ascending order. |
|
* |
|
* PERM (output) INTEGER array, dimension ( N ) |
|
* The permutations (from deflation and sorting) to be applied |
|
* to each block. Not referenced if ICOMPQ = 0. |
|
* |
|
* GIVPTR (output) INTEGER |
|
* The number of Givens rotations which took place in this |
|
* subproblem. Not referenced if ICOMPQ = 0. |
|
* |
|
* GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) |
|
* Each pair of numbers indicates a pair of columns to take place |
|
* in a Givens rotation. Not referenced if ICOMPQ = 0. |
|
* |
|
* LDGCOL (input) INTEGER |
|
* leading dimension of GIVCOL, must be at least N. |
|
* |
|
* GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) |
|
* Each number indicates the C or S value to be used in the |
|
* corresponding Givens rotation. Not referenced if ICOMPQ = 0. |
|
* |
|
* LDGNUM (input) INTEGER |
|
* The leading dimension of GIVNUM and POLES, must be at least N. |
|
* |
|
* POLES (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) |
|
* On exit, POLES(1,*) is an array containing the new singular |
|
* values obtained from solving the secular equation, and |
|
* POLES(2,*) is an array containing the poles in the secular |
|
* equation. Not referenced if ICOMPQ = 0. |
|
* |
|
* DIFL (output) DOUBLE PRECISION array, dimension ( N ) |
|
* On exit, DIFL(I) is the distance between I-th updated |
|
* (undeflated) singular value and the I-th (undeflated) old |
|
* singular value. |
|
* |
|
* DIFR (output) DOUBLE PRECISION array, |
|
* dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and |
|
* dimension ( N ) if ICOMPQ = 0. |
|
* On exit, DIFR(I, 1) is the distance between I-th updated |
|
* (undeflated) singular value and the I+1-th (undeflated) old |
|
* singular value. |
|
* |
|
* If ICOMPQ = 1, DIFR(1:K,2) is an array containing the |
|
* normalizing factors for the right singular vector matrix. |
|
* |
|
* See DLASD8 for details on DIFL and DIFR. |
|
* |
|
* Z (output) DOUBLE PRECISION array, dimension ( M ) |
|
* The first elements of this array contain the components |
|
* of the deflation-adjusted updating row vector. |
|
* |
|
* K (output) INTEGER |
|
* Contains the dimension of the non-deflated matrix, |
|
* This is the order of the related secular equation. 1 <= K <=N. |
|
* |
|
* C (output) DOUBLE PRECISION |
|
* C contains garbage if SQRE =0 and the C-value of a Givens |
|
* rotation related to the right null space if SQRE = 1. |
|
* |
|
* S (output) DOUBLE PRECISION |
|
* S contains garbage if SQRE =0 and the S-value of a Givens |
|
* rotation related to the right null space if SQRE = 1. |
|
* |
|
* WORK (workspace) DOUBLE PRECISION array, dimension ( 4 * M ) |
|
* |
|
* IWORK (workspace) INTEGER array, dimension ( 3 * N ) |
|
* |
|
* INFO (output) INTEGER |
|
* = 0: successful exit. |
|
* < 0: if INFO = -i, the i-th argument had an illegal value. |
|
* > 0: if INFO = 1, an singular value did not converge |
|
* |
|
* Further Details |
|
* =============== |
|
* |
|
* Based on contributions by |
|
* Ming Gu and Huan Ren, Computer Science Division, University of |
|
* California at Berkeley, USA |
|
* |
|
* ===================================================================== |
* ===================================================================== |
* |
* |
* .. Parameters .. |
* .. Parameters .. |
Line 282
|
Line 413
|
CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, |
CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, |
$ WORK( ISIGMA ), WORK( IW ), INFO ) |
$ WORK( ISIGMA ), WORK( IW ), INFO ) |
* |
* |
|
* Report the possible convergence failure. |
|
* |
|
IF( INFO.NE.0 ) THEN |
|
RETURN |
|
END IF |
|
* |
* Save the poles if ICOMPQ = 1. |
* Save the poles if ICOMPQ = 1. |
* |
* |
IF( ICOMPQ.EQ.1 ) THEN |
IF( ICOMPQ.EQ.1 ) THEN |