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*> \date June 2016 301* 302*> \ingroup OTHERauxiliary 303* 304*> \par Contributors: 305* ================== 306*> 307*> Ming Gu and Huan Ren, Computer Science Division, University of 308*> California at Berkeley, USA 309*> 310* ===================================================================== 311 SUBROUTINE DLASD6( ICOMPQ, NL, NR, SQRE, D, VF, VL, ALPHA, BETA, 312 $ IDXQ, PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, 313 $ LDGNUM, POLES, DIFL, DIFR, Z, K, C, S, WORK, 314 $ IWORK, INFO ) 315* 316* -- LAPACK auxiliary routine (version 3.7.0) -- 317* -- LAPACK is a software package provided by Univ. of Tennessee, -- 318* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 319* June 2016 320* 321* .. Scalar Arguments .. 322 INTEGER GIVPTR, ICOMPQ, INFO, K, LDGCOL, LDGNUM, NL, 323 $ NR, SQRE 324 DOUBLE PRECISION ALPHA, BETA, C, S 325* .. 326* .. Array Arguments .. 327 INTEGER GIVCOL( LDGCOL, * ), IDXQ( * ), IWORK( * ), 328 $ PERM( * ) 329 DOUBLE PRECISION D( * ), DIFL( * ), DIFR( * ), 330 $ GIVNUM( LDGNUM, * ), POLES( LDGNUM, * ), 331 $ VF( * ), VL( * ), WORK( * ), Z( * ) 332* .. 333* 334* ===================================================================== 335* 336* .. Parameters .. 337 DOUBLE PRECISION ONE, ZERO 338 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 339* .. 340* .. Local Scalars .. 341 INTEGER I, IDX, IDXC, IDXP, ISIGMA, IVFW, IVLW, IW, M, 342 $ N, N1, N2 343 DOUBLE PRECISION ORGNRM 344* .. 345* .. External Subroutines .. 346 EXTERNAL DCOPY, DLAMRG, DLASCL, DLASD7, DLASD8, XERBLA 347* .. 348* .. Intrinsic Functions .. 349 INTRINSIC ABS, MAX 350* .. 351* .. Executable Statements .. 352* 353* Test the input parameters. 354* 355 INFO = 0 356 N = NL + NR + 1 357 M = N + SQRE 358* 359 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 360 INFO = -1 361 ELSE IF( NL.LT.1 ) THEN 362 INFO = -2 363 ELSE IF( NR.LT.1 ) THEN 364 INFO = -3 365 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 366 INFO = -4 367 ELSE IF( LDGCOL.LT.N ) THEN 368 INFO = -14 369 ELSE IF( LDGNUM.LT.N ) THEN 370 INFO = -16 371 END IF 372 IF( INFO.NE.0 ) THEN 373 CALL XERBLA( 'DLASD6', -INFO ) 374 RETURN 375 END IF 376* 377* The following values are for bookkeeping purposes only. They are 378* integer pointers which indicate the portion of the workspace 379* used by a particular array in DLASD7 and DLASD8. 380* 381 ISIGMA = 1 382 IW = ISIGMA + N 383 IVFW = IW + M 384 IVLW = IVFW + M 385* 386 IDX = 1 387 IDXC = IDX + N 388 IDXP = IDXC + N 389* 390* Scale. 391* 392 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 393 D( NL+1 ) = ZERO 394 DO 10 I = 1, N 395 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 396 ORGNRM = ABS( D( I ) ) 397 END IF 398 10 CONTINUE 399 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 400 ALPHA = ALPHA / ORGNRM 401 BETA = BETA / ORGNRM 402* 403* Sort and Deflate singular values. 404* 405 CALL DLASD7( ICOMPQ, NL, NR, SQRE, K, D, Z, WORK( IW ), VF, 406 $ WORK( IVFW ), VL, WORK( IVLW ), ALPHA, BETA, 407 $ WORK( ISIGMA ), IWORK( IDX ), IWORK( IDXP ), IDXQ, 408 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, C, S, 409 $ INFO ) 410* 411* Solve Secular Equation, compute DIFL, DIFR, and update VF, VL. 412* 413 CALL DLASD8( ICOMPQ, K, D, Z, VF, VL, DIFL, DIFR, LDGNUM, 414 $ WORK( ISIGMA ), WORK( IW ), INFO ) 415* 416* Report the possible convergence failure. 417* 418 IF( INFO.NE.0 ) THEN 419 RETURN 420 END IF 421* 422* Save the poles if ICOMPQ = 1. 423* 424 IF( ICOMPQ.EQ.1 ) THEN 425 CALL DCOPY( K, D, 1, POLES( 1, 1 ), 1 ) 426 CALL DCOPY( K, WORK( ISIGMA ), 1, POLES( 1, 2 ), 1 ) 427 END IF 428* 429* Unscale. 430* 431 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 432* 433* Prepare the IDXQ sorting permutation. 434* 435 N1 = K 436 N2 = N - K 437 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 438* 439 RETURN 440* 441* End of DLASD6 442* 443 END 444