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. 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 occurrence 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*> 167*> \param[in,out] IDXQ 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 ( LDDIFR, 2 ) if ICOMPQ = 1 and 236*> dimension ( K ) if ICOMPQ = 0. 237*> On exit, DIFR(I,1) = D(I) - DSIGMA(I+1), DIFR(K,1) is not 238*> defined and will not be referenced. 239*> 240*> If ICOMPQ = 1, DIFR(1:K,2) is an array containing the 241*> normalizing factors for the right singular vector matrix. 242*> 243*> See DLASD8 for details on DIFL and DIFR. 244*> \endverbatim 245*> 246*> \param[out] Z 247*> \verbatim 248*> Z is DOUBLE PRECISION array, dimension ( M ) 249*> The first elements of this array contain the components 250*> of the deflation-adjusted updating row vector. 251*> \endverbatim 252*> 253*> \param[out] K 254*> \verbatim 255*> K is INTEGER 256*> Contains the dimension of the non-deflated matrix, 257*> This is the order of the related secular equation. 1 <= K <=N. 258*> \endverbatim 259*> 260*> \param[out] C 261*> \verbatim 262*> C is DOUBLE PRECISION 263*> C contains garbage if SQRE =0 and the C-value of a Givens 264*> rotation related to the right null space if SQRE = 1. 265*> \endverbatim 266*> 267*> \param[out] S 268*> \verbatim 269*> S is DOUBLE PRECISION 270*> S contains garbage if SQRE =0 and the S-value of a Givens 271*> rotation related to the right null space if SQRE = 1. 272*> \endverbatim 273*> 274*> \param[out] WORK 275*> \verbatim 276*> WORK is DOUBLE PRECISION array, dimension ( 4 * M ) 277*> \endverbatim 278*> 279*> \param[out] IWORK 280*> \verbatim 281*> IWORK is INTEGER array, dimension ( 3 * N ) 282*> \endverbatim 283*> 284*> \param[out] INFO 285*> \verbatim 286*> INFO is INTEGER 287*> = 0: successful exit. 288*> < 0: if INFO = -i, the i-th argument had an illegal value. 289*> > 0: if INFO = 1, a singular value did not converge 290*> \endverbatim 291* 292* Authors: 293* ======== 294* 295*> \author Univ. of Tennessee 296*> \author Univ. of California Berkeley 297*> \author Univ. of Colorado Denver 298*> \author NAG Ltd. 299* 300*> \ingroup OTHERauxiliary 301* 302*> \par Contributors: 303* ================== 304*> 305*> Ming Gu and Huan Ren, Computer Science Division, University of 306*> California at Berkeley, USA 307*> 308* ===================================================================== 309 SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 310 $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 311 $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 312 $ IWORK, INFO ) 313* 314* -- LAPACK auxiliary routine -- 315* -- LAPACK is a software package provided by Univ. of Tennessee, -- 316* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 317* 318* .. Scalar Arguments .. 319 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 320 $ NR, SQRE 321 DOUBLE PRECISION ALPHA, BETA, C, S 322* .. 323* .. Array Arguments .. 324 INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 325 $ PERM( * ) 326 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ), 327 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 328 $ VF( * ), VL( * ), WORK( * ), Z( * ) 329* .. 330* 331* ===================================================================== 332* 333* .. Parameters .. 334 DOUBLE PRECISION ONE, ZERO 335 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 336* .. 337* .. Local Scalars .. 338 INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M, 339 $ N, N1, N2 340 DOUBLE PRECISION ORGNRM 341* .. 342* .. External Subroutines .. 343 EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA 344* .. 345* .. Intrinsic Functions .. 346 INTRINSIC ABS, MAX 347* .. 348* .. Executable Statements .. 349* 350* Test the input parameters. 351* 352 INFO = 0 353 N = NL + NR + 1 354 M = N + SQRE 355* 356 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 357 INFO = -1 358 ELSE IF( NL.LT.1 ) THEN 359 INFO = -2 360 ELSE IF( NR.LT.1 ) THEN 361 INFO = -3 362 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 363 INFO = -4 364 ELSE IF( LDGCOL.LT.N ) THEN 365 INFO = -14 366 ELSE IF( LDGNUM.LT.N ) THEN 367 INFO = -16 368 END IF 369 IF( INFO.NE.0 ) THEN 370 CALL XERBLA( 'DLASD6', -INFO ) 371 RETURN 372 END IF 373* 374* The following values are for bookkeeping purposes only. They are 375* integer pointers which indicate the portion of the workspace 376* used by a particular array in DLASD7 and DLASD8. 377* 378 ISIGMA = 1 379 IW = ISIGMA + N 380 IVFW = IW + M 381 IVLW = IVFW + M 382* 383 IDX = 1 384 IDXC = IDX + N 385 IDXP = IDXC + N 386* 387* Scale. 388* 389 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 390 D( NL+1 ) = ZERO 391 DO 10 I = 1, N 392 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 393 ORGNRM = ABS( D( I ) ) 394 END IF 395 10 CONTINUE 396 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 397 ALPHA = ALPHA / ORGNRM 398 BETA = BETA / ORGNRM 399* 400* Sort and Deflate singular values. 401* 402 CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF, 403 $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA, 404 $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ, 405 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 406 $ INFO ) 407* 408* Solve Secular Equation, compute DIFL, DIFR, and update VF, VL. 409* 410 CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, 411 $ WORK( ISIGMA ), WORK( IW ), INFO ) 412* 413* Report the possible convergence failure. 414* 415 IF( INFO.NE.0 ) THEN 416 RETURN 417 END IF 418* 419* Save the poles if ICOMPQ = 1. 420* 421 IF( ICOMPQ.EQ.1 ) THEN 422 CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 ) 423 CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 ) 424 END IF 425* 426* Unscale. 427* 428 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 429* 430* Prepare the IDXQ sorting permutation. 431* 432 N1 = K 433 N2 = N - K 434 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 435* 436 RETURN 437* 438* End of DLASD6 439* 440 END 441