1*> \brief \b SLASD1 computes the SVD of an upper bidiagonal matrix B of the specified size. 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 SLASD1 + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slasd1.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slasd1.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slasd1.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 22* IDXQ, IWORK, WORK, INFO ) 23* 24* .. Scalar Arguments .. 25* INTEGER INFO, LDU, LDVT, NL, NR, SQRE 26* REAL ALPHA, BETA 27* .. 28* .. Array Arguments .. 29* INTEGER IDXQ( * ), IWORK( * ) 30* REAL D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, 40*> where N = NL + NR + 1 and M = N + SQRE. SLASD1 is called from SLASD0. 41*> 42*> A related subroutine SLASD7 handles the case in which the singular 43*> values (and the singular vectors in factored form) are desired. 44*> 45*> SLASD1 computes the SVD as follows: 46*> 47*> ( D1(in) 0 0 0 ) 48*> B = U(in) * ( Z1**T a Z2**T b ) * VT(in) 49*> ( 0 0 D2(in) 0 ) 50*> 51*> = U(out) * ( D(out) 0) * VT(out) 52*> 53*> where Z**T = (Z1**T a Z2**T b) = u**T VT**T, and u is a vector of dimension M 54*> with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 55*> elsewhere; and the entry b is empty if SQRE = 0. 56*> 57*> The left singular vectors of the original matrix are stored in U, and 58*> the transpose of the right singular vectors are stored in VT, and the 59*> singular values are in D. The algorithm consists of three stages: 60*> 61*> The first stage consists of deflating the size of the problem 62*> when there are multiple singular values or when there are zeros in 63*> the Z vector. For each such occurrence the dimension of the 64*> secular equation problem is reduced by one. This stage is 65*> performed by the routine SLASD2. 66*> 67*> The second stage consists of calculating the updated 68*> singular values. This is done by finding the square roots of the 69*> roots of the secular equation via the routine SLASD4 (as called 70*> by SLASD3). This routine also calculates the singular vectors of 71*> the current problem. 72*> 73*> The final stage consists of computing the updated singular vectors 74*> directly using the updated singular values. The singular vectors 75*> for the current problem are multiplied with the singular vectors 76*> from the overall problem. 77*> \endverbatim 78* 79* Arguments: 80* ========== 81* 82*> \param[in] NL 83*> \verbatim 84*> NL is INTEGER 85*> The row dimension of the upper block. NL >= 1. 86*> \endverbatim 87*> 88*> \param[in] NR 89*> \verbatim 90*> NR is INTEGER 91*> The row dimension of the lower block. NR >= 1. 92*> \endverbatim 93*> 94*> \param[in] SQRE 95*> \verbatim 96*> SQRE is INTEGER 97*> = 0: the lower block is an NR-by-NR square matrix. 98*> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 99*> 100*> The bidiagonal matrix has row dimension N = NL + NR + 1, 101*> and column dimension M = N + SQRE. 102*> \endverbatim 103*> 104*> \param[in,out] D 105*> \verbatim 106*> D is REAL array, dimension (NL+NR+1). 107*> N = NL+NR+1 108*> On entry D(1:NL,1:NL) contains the singular values of the 109*> upper block; and D(NL+2:N) contains the singular values of 110*> the lower block. On exit D(1:N) contains the singular values 111*> of the modified matrix. 112*> \endverbatim 113*> 114*> \param[in,out] ALPHA 115*> \verbatim 116*> ALPHA is REAL 117*> Contains the diagonal element associated with the added row. 118*> \endverbatim 119*> 120*> \param[in,out] BETA 121*> \verbatim 122*> BETA is REAL 123*> Contains the off-diagonal element associated with the added 124*> row. 125*> \endverbatim 126*> 127*> \param[in,out] U 128*> \verbatim 129*> U is REAL array, dimension (LDU,N) 130*> On entry U(1:NL, 1:NL) contains the left singular vectors of 131*> the upper block; U(NL+2:N, NL+2:N) contains the left singular 132*> vectors of the lower block. On exit U contains the left 133*> singular vectors of the bidiagonal matrix. 134*> \endverbatim 135*> 136*> \param[in] LDU 137*> \verbatim 138*> LDU is INTEGER 139*> The leading dimension of the array U. LDU >= max( 1, N ). 140*> \endverbatim 141*> 142*> \param[in,out] VT 143*> \verbatim 144*> VT is REAL array, dimension (LDVT,M) 145*> where M = N + SQRE. 146*> On entry VT(1:NL+1, 1:NL+1)**T contains the right singular 147*> vectors of the upper block; VT(NL+2:M, NL+2:M)**T contains 148*> the right singular vectors of the lower block. On exit 149*> VT**T contains the right singular vectors of the 150*> bidiagonal matrix. 151*> \endverbatim 152*> 153*> \param[in] LDVT 154*> \verbatim 155*> LDVT is INTEGER 156*> The leading dimension of the array VT. LDVT >= max( 1, M ). 157*> \endverbatim 158*> 159*> \param[in,out] IDXQ 160*> \verbatim 161*> IDXQ is INTEGER array, dimension (N) 162*> This contains the permutation which will reintegrate the 163*> subproblem just solved back into sorted order, i.e. 164*> D( IDXQ( I = 1, N ) ) will be in ascending order. 165*> \endverbatim 166*> 167*> \param[out] IWORK 168*> \verbatim 169*> IWORK is INTEGER array, dimension (4*N) 170*> \endverbatim 171*> 172*> \param[out] WORK 173*> \verbatim 174*> WORK is REAL array, dimension (3*M**2+2*M) 175*> \endverbatim 176*> 177*> \param[out] INFO 178*> \verbatim 179*> INFO is INTEGER 180*> = 0: successful exit. 181*> < 0: if INFO = -i, the i-th argument had an illegal value. 182*> > 0: if INFO = 1, a singular value did not converge 183*> \endverbatim 184* 185* Authors: 186* ======== 187* 188*> \author Univ. of Tennessee 189*> \author Univ. of California Berkeley 190*> \author Univ. of Colorado Denver 191*> \author NAG Ltd. 192* 193*> \ingroup OTHERauxiliary 194* 195*> \par Contributors: 196* ================== 197*> 198*> Ming Gu and Huan Ren, Computer Science Division, University of 199*> California at Berkeley, USA 200*> 201* ===================================================================== 202 SUBROUTINE SLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 203 $ IDXQ, IWORK, WORK, INFO ) 204* 205* -- LAPACK auxiliary routine -- 206* -- LAPACK is a software package provided by Univ. of Tennessee, -- 207* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 208* 209* .. Scalar Arguments .. 210 INTEGER INFO, LDU, LDVT, NL, NR, SQRE 211 REAL ALPHA, BETA 212* .. 213* .. Array Arguments .. 214 INTEGER IDXQ( * ), IWORK( * ) 215 REAL D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 216* .. 217* 218* ===================================================================== 219* 220* .. Parameters .. 221* 222 REAL ONE, ZERO 223 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 224* .. 225* .. Local Scalars .. 226 INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2, 227 $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2 228 REAL ORGNRM 229* .. 230* .. External Subroutines .. 231 EXTERNAL SLAMRG, SLASCL, SLASD2, SLASD3, XERBLA 232* .. 233* .. Intrinsic Functions .. 234 INTRINSIC ABS, MAX 235* .. 236* .. Executable Statements .. 237* 238* Test the input parameters. 239* 240 INFO = 0 241* 242 IF( NL.LT.1 ) THEN 243 INFO = -1 244 ELSE IF( NR.LT.1 ) THEN 245 INFO = -2 246 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 247 INFO = -3 248 END IF 249 IF( INFO.NE.0 ) THEN 250 CALL XERBLA( 'SLASD1', -INFO ) 251 RETURN 252 END IF 253* 254 N = NL + NR + 1 255 M = N + SQRE 256* 257* The following values are for bookkeeping purposes only. They are 258* integer pointers which indicate the portion of the workspace 259* used by a particular array in SLASD2 and SLASD3. 260* 261 LDU2 = N 262 LDVT2 = M 263* 264 IZ = 1 265 ISIGMA = IZ + M 266 IU2 = ISIGMA + N 267 IVT2 = IU2 + LDU2*N 268 IQ = IVT2 + LDVT2*M 269* 270 IDX = 1 271 IDXC = IDX + N 272 COLTYP = IDXC + N 273 IDXP = COLTYP + N 274* 275* Scale. 276* 277 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 278 D( NL+1 ) = ZERO 279 DO 10 I = 1, N 280 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 281 ORGNRM = ABS( D( I ) ) 282 END IF 283 10 CONTINUE 284 CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 285 ALPHA = ALPHA / ORGNRM 286 BETA = BETA / ORGNRM 287* 288* Deflate singular values. 289* 290 CALL SLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU, 291 $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2, 292 $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ), 293 $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO ) 294* 295* Solve Secular Equation and update singular vectors. 296* 297 LDQ = K 298 CALL SLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ), 299 $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ), 300 $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ), 301 $ INFO ) 302* 303* Report the possible convergence failure. 304* 305 IF( INFO.NE.0 ) THEN 306 RETURN 307 END IF 308* 309* Unscale. 310* 311 CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 312* 313* Prepare the IDXQ sorting permutation. 314* 315 N1 = K 316 N2 = N - K 317 CALL SLAMRG( N1, N2, D, 1, -1, IDXQ ) 318* 319 RETURN 320* 321* End of SLASD1 322* 323 END 324