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