1 SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 2 $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 3 $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 4 $ IWORK, INFO ) 5* 6* -- LAPACK auxiliary routine (version 3.0) -- 7* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 8* Courant Institute, Argonne National Lab, and Rice University 9* June 30, 1999 10* 11* .. Scalar Arguments .. 12 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 13 $ NR, SQRE 14 DOUBLE PRECISION ALPHA, BETA, C, S 15* .. 16* .. Array Arguments .. 17 INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 18 $ PERM( * ) 19 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ), 20 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 21 $ VF( * ), VL( * ), WORK( * ), Z( * ) 22* .. 23* 24* Purpose 25* ======= 26* 27* DLASD6 computes the SVD of an updated upper bidiagonal matrix B 28* obtained by merging two smaller ones by appending a row. This 29* routine is used only for the problem which requires all singular 30* values and optionally singular vector matrices in factored form. 31* B is an N-by-M matrix with N = NL + NR + 1 and M = N + SQRE. 32* A related subroutine, DLASD1, handles the case in which all singular 33* values and singular vectors of the bidiagonal matrix are desired. 34* 35* DLASD6 computes the SVD as follows: 36* 37* ( D1(in) 0 0 0 ) 38* B = U(in) * ( Z1' a Z2' b ) * VT(in) 39* ( 0 0 D2(in) 0 ) 40* 41* = U(out) * ( D(out) 0) * VT(out) 42* 43* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M 44* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 45* elsewhere; and the entry b is empty if SQRE = 0. 46* 47* The singular values of B can be computed using D1, D2, the first 48* components of all the right singular vectors of the lower block, and 49* the last components of all the right singular vectors of the upper 50* block. These components are stored and updated in VF and VL, 51* respectively, in DLASD6. Hence U and VT are not explicitly 52* referenced. 53* 54* The singular values are stored in D. The algorithm consists of two 55* stages: 56* 57* The first stage consists of deflating the size of the problem 58* when there are multiple singular values or if there is a zero 59* in the Z vector. For each such occurence the dimension of the 60* secular equation problem is reduced by one. This stage is 61* performed by the routine DLASD7. 62* 63* The second stage consists of calculating the updated 64* singular values. This is done by finding the roots of the 65* secular equation via the routine DLASD4 (as called by DLASD8). 66* This routine also updates VF and VL and computes the distances 67* between the updated singular values and the old singular 68* values. 69* 70* DLASD6 is called from DLASDA. 71* 72* Arguments 73* ========= 74* 75* ICOMPQ (input) INTEGER 76* Specifies whether singular vectors are to be computed in 77* factored form: 78* = 0: Compute singular values only. 79* = 1: Compute singular vectors in factored form as well. 80* 81* NL (input) INTEGER 82* The row dimension of the upper block. NL >= 1. 83* 84* NR (input) INTEGER 85* The row dimension of the lower block. NR >= 1. 86* 87* SQRE (input) INTEGER 88* = 0: the lower block is an NR-by-NR square matrix. 89* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 90* 91* The bidiagonal matrix has row dimension N = NL + NR + 1, 92* and column dimension M = N + SQRE. 93* 94* D (input/output) DOUBLE PRECISION array, dimension ( NL+NR+1 ). 95* On entry D(1:NL,1:NL) contains the singular values of the 96* upper block, and D(NL+2:N) contains the singular values 97* of the lower block. On exit D(1:N) contains the singular 98* values of the modified matrix. 99* 100* VF (input/output) DOUBLE PRECISION array, dimension ( M ) 101* On entry, VF(1:NL+1) contains the first components of all 102* right singular vectors of the upper block; and VF(NL+2:M) 103* contains the first components of all right singular vectors 104* of the lower block. On exit, VF contains the first components 105* of all right singular vectors of the bidiagonal matrix. 106* 107* VL (input/output) DOUBLE PRECISION array, dimension ( M ) 108* On entry, VL(1:NL+1) contains the last components of all 109* right singular vectors of the upper block; and VL(NL+2:M) 110* contains the last components of all right singular vectors of 111* the lower block. On exit, VL contains the last components of 112* all right singular vectors of the bidiagonal matrix. 113* 114* ALPHA (input) DOUBLE PRECISION 115* Contains the diagonal element associated with the added row. 116* 117* BETA (input) DOUBLE PRECISION 118* Contains the off-diagonal element associated with the added 119* row. 120* 121* IDXQ (output) INTEGER array, dimension ( N ) 122* This contains the permutation which will reintegrate the 123* subproblem just solved back into sorted order, i.e. 124* D( IDXQ( I = 1, N ) ) will be in ascending order. 125* 126* PERM (output) INTEGER array, dimension ( N ) 127* The permutations (from deflation and sorting) to be applied 128* to each block. Not referenced if ICOMPQ = 0. 129* 130* GIVPTR (output) INTEGER 131* The number of Givens rotations which took place in this 132* subproblem. Not referenced if ICOMPQ = 0. 133* 134* GIVCOL (output) INTEGER array, dimension ( LDGCOL, 2 ) 135* Each pair of numbers indicates a pair of columns to take place 136* in a Givens rotation. Not referenced if ICOMPQ = 0. 137* 138* LDGCOL (input) INTEGER 139* leading dimension of GIVCOL, must be at least N. 140* 141* GIVNUM (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 142* Each number indicates the C or S value to be used in the 143* corresponding Givens rotation. Not referenced if ICOMPQ = 0. 144* 145* LDGNUM (input) INTEGER 146* The leading dimension of GIVNUM and POLES, must be at least N. 147* 148* POLES (output) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 149* On exit, POLES(1,*) is an array containing the new singular 150* values obtained from solving the secular equation, and 151* POLES(2,*) is an array containing the poles in the secular 152* equation. Not referenced if ICOMPQ = 0. 153* 154* DIFL (output) DOUBLE PRECISION array, dimension ( N ) 155* On exit, DIFL(I) is the distance between I-th updated 156* (undeflated) singular value and the I-th (undeflated) old 157* singular value. 158* 159* DIFR (output) DOUBLE PRECISION array, 160* dimension ( LDGNUM, 2 ) if ICOMPQ = 1 and 161* dimension ( N ) if ICOMPQ = 0. 162* On exit, DIFR(I, 1) is the distance between I-th updated 163* (undeflated) singular value and the I+1-th (undeflated) old 164* singular value. 165* 166* If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 167* normalizing factors for the right singular vector matrix. 168* 169* See DLASD8 for details on DIFL and DIFR. 170* 171* Z (output) DOUBLE PRECISION array, dimension ( M ) 172* The first elements of this array contain the components 173* of the deflation-adjusted updating row vector. 174* 175* K (output) INTEGER 176* Contains the dimension of the non-deflated matrix, 177* This is the order of the related secular equation. 1 <= K <=N. 178* 179* C (output) DOUBLE PRECISION 180* C contains garbage if SQRE =0 and the C-value of a Givens 181* rotation related to the right null space if SQRE = 1. 182* 183* S (output) DOUBLE PRECISION 184* S contains garbage if SQRE =0 and the S-value of a Givens 185* rotation related to the right null space if SQRE = 1. 186* 187* WORK (workspace) DOUBLE PRECISION array, dimension ( 4 * M ) 188* 189* IWORK (workspace) INTEGER array, dimension ( 3 * N ) 190* 191* INFO (output) INTEGER 192* = 0: successful exit. 193* < 0: if INFO = -i, the i-th argument had an illegal value. 194* > 0: if INFO = 1, an singular value did not converge 195* 196* Further Details 197* =============== 198* 199* Based on contributions by 200* Ming Gu and Huan Ren, Computer Science Division, University of 201* California at Berkeley, USA 202* 203* ===================================================================== 204* 205* .. Parameters .. 206 DOUBLE PRECISION ONE, ZERO 207 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 208* .. 209* .. Local Scalars .. 210 INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M, 211 $ N, N1, N2 212 DOUBLE PRECISION ORGNRM 213* .. 214* .. External Subroutines .. 215 EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA 216* .. 217* .. Intrinsic Functions .. 218 INTRINSIC ABS, MAX 219* .. 220* .. Executable Statements .. 221* 222* Test the input parameters. 223* 224 INFO = 0 225 N = NL + NR + 1 226 M = N + SQRE 227* 228 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 229 INFO = -1 230 ELSE IF( NL.LT.1 ) THEN 231 INFO = -2 232 ELSE IF( NR.LT.1 ) THEN 233 INFO = -3 234 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 235 INFO = -4 236 ELSE IF( LDGCOL.LT.N ) THEN 237 INFO = -14 238 ELSE IF( LDGNUM.LT.N ) THEN 239 INFO = -16 240 END IF 241 IF( INFO.NE.0 ) THEN 242 CALL XERBLA( 'DLASD6', -INFO ) 243 RETURN 244 END IF 245* 246* The following values are for bookkeeping purposes only. They are 247* integer pointers which indicate the portion of the workspace 248* used by a particular array in DLASD7 and DLASD8. 249* 250 ISIGMA = 1 251 IW = ISIGMA + N 252 IVFW = IW + M 253 IVLW = IVFW + M 254* 255 IDX = 1 256 IDXC = IDX + N 257 IDXP = IDXC + N 258* 259* Scale. 260* 261 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 262 D( NL+1 ) = ZERO 263 DO 10 I = 1, N 264 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 265 ORGNRM = ABS( D( I ) ) 266 END IF 267 10 CONTINUE 268 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 269 ALPHA = ALPHA / ORGNRM 270 BETA = BETA / ORGNRM 271* 272* Sort and Deflate singular values. 273* 274 CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF, 275 $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA, 276 $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ, 277 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 278 $ INFO ) 279* 280* Solve Secular Equation, compute DIFL, DIFR, and update VF, VL. 281* 282 CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, 283 $ WORK( ISIGMA ), WORK( IW ), INFO ) 284* 285* Save the poles if ICOMPQ = 1. 286* 287 IF( ICOMPQ.EQ.1 ) THEN 288 CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 ) 289 CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 ) 290 END IF 291* 292* Unscale. 293* 294 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 295* 296* Prepare the IDXQ sorting permutation. 297* 298 N1 = K 299 N2 = N - K 300 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 301* 302 RETURN 303* 304* End of DLASD6 305* 306 END 307