1 SUBROUTINE DLASD1( NL, NR, SQRE, D, ALPHA, BETA, U, LDU, VT, LDVT, 2 $ IDXQ, IWORK, WORK, INFO ) 3* 4* -- LAPACK auxiliary routine (instrumented to count ops, version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 INTEGER INFO, LDU, LDVT, NL, NR, SQRE 11 DOUBLE PRECISION ALPHA, BETA 12* .. 13* .. Array Arguments .. 14 INTEGER IDXQ( * ), IWORK( * ) 15 DOUBLE PRECISION D( * ), U( LDU, * ), VT( LDVT, * ), WORK( * ) 16* .. 17* .. Common block to return operation count .. 18 COMMON / LATIME / OPS, ITCNT 19* .. 20* .. Scalars in Common .. 21 DOUBLE PRECISION ITCNT, OPS 22* .. 23* 24* Purpose 25* ======= 26* 27* DLASD1 computes the SVD of an upper bidiagonal N-by-M matrix B, 28* where N = NL + NR + 1 and M = N + SQRE. DLASD1 is called from DLASD0. 29* 30* A related subroutine DLASD7 handles the case in which the singular 31* values (and the singular vectors in factored form) are desired. 32* 33* DLASD1 computes the SVD as follows: 34* 35* ( D1(in) 0 0 0 ) 36* B = U(in) * ( Z1' a Z2' b ) * VT(in) 37* ( 0 0 D2(in) 0 ) 38* 39* = U(out) * ( D(out) 0) * VT(out) 40* 41* where Z' = (Z1' a Z2' b) = u' VT', and u is a vector of dimension M 42* with ALPHA and BETA in the NL+1 and NL+2 th entries and zeros 43* elsewhere; and the entry b is empty if SQRE = 0. 44* 45* The left singular vectors of the original matrix are stored in U, and 46* the transpose of the right singular vectors are stored in VT, and the 47* singular values are in D. The algorithm consists of three stages: 48* 49* The first stage consists of deflating the size of the problem 50* when there are multiple singular values or when there are zeros in 51* the Z vector. For each such occurence the dimension of the 52* secular equation problem is reduced by one. This stage is 53* performed by the routine DLASD2. 54* 55* The second stage consists of calculating the updated 56* singular values. This is done by finding the square roots of the 57* roots of the secular equation via the routine DLASD4 (as called 58* by DLASD3). This routine also calculates the singular vectors of 59* the current problem. 60* 61* The final stage consists of computing the updated singular vectors 62* directly using the updated singular values. The singular vectors 63* for the current problem are multiplied with the singular vectors 64* from the overall problem. 65* 66* Arguments 67* ========= 68* 69* NL (input) INTEGER 70* The row dimension of the upper block. NL >= 1. 71* 72* NR (input) INTEGER 73* The row dimension of the lower block. NR >= 1. 74* 75* SQRE (input) INTEGER 76* = 0: the lower block is an NR-by-NR square matrix. 77* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 78* 79* The bidiagonal matrix has row dimension N = NL + NR + 1, 80* and column dimension M = N + SQRE. 81* 82* D (input/output) DOUBLE PRECISION array, 83* dimension (N = NL+NR+1). 84* On entry D(1:NL,1:NL) contains the singular values of the 85* upper block; and D(NL+2:N) contains the singular values of 86* the lower block. On exit D(1:N) contains the singular values 87* of the modified matrix. 88* 89* ALPHA (input) DOUBLE PRECISION 90* Contains the diagonal element associated with the added row. 91* 92* BETA (input) DOUBLE PRECISION 93* Contains the off-diagonal element associated with the added 94* row. 95* 96* U (input/output) DOUBLE PRECISION array, dimension(LDU,N) 97* On entry U(1:NL, 1:NL) contains the left singular vectors of 98* the upper block; U(NL+2:N, NL+2:N) contains the left singular 99* vectors of the lower block. On exit U contains the left 100* singular vectors of the bidiagonal matrix. 101* 102* LDU (input) INTEGER 103* The leading dimension of the array U. LDU >= max( 1, N ). 104* 105* VT (input/output) DOUBLE PRECISION array, dimension(LDVT,M) 106* where M = N + SQRE. 107* On entry VT(1:NL+1, 1:NL+1)' contains the right singular 108* vectors of the upper block; VT(NL+2:M, NL+2:M)' contains 109* the right singular vectors of the lower block. On exit 110* VT' contains the right singular vectors of the 111* bidiagonal matrix. 112* 113* LDVT (input) INTEGER 114* The leading dimension of the array VT. LDVT >= max( 1, M ). 115* 116* IDXQ (output) INTEGER array, dimension(N) 117* This contains the permutation which will reintegrate the 118* subproblem just solved back into sorted order, i.e. 119* D( IDXQ( I = 1, N ) ) will be in ascending order. 120* 121* IWORK (workspace) INTEGER array, dimension( 4 * N ) 122* 123* WORK (workspace) DOUBLE PRECISION array, dimension( 3*M**2 + 2*M ) 124* 125* INFO (output) INTEGER 126* = 0: successful exit. 127* < 0: if INFO = -i, the i-th argument had an illegal value. 128* > 0: if INFO = 1, an singular value did not converge 129* 130* Further Details 131* =============== 132* 133* Based on contributions by 134* Ming Gu and Huan Ren, Computer Science Division, University of 135* California at Berkeley, USA 136* 137* ===================================================================== 138* 139* .. Parameters .. 140* 141 DOUBLE PRECISION ONE, ZERO 142 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0 ) 143* .. 144* .. Local Scalars .. 145 INTEGER COLTYP, I, IDX, IDXC, IDXP, IQ, ISIGMA, IU2, 146 $ IVT2, IZ, K, LDQ, LDU2, LDVT2, M, N, N1, N2 147 DOUBLE PRECISION ORGNRM 148* .. 149* .. External Subroutines .. 150 EXTERNAL DLAMRG, DLASCL, DLASD2, DLASD3, XERBLA 151* .. 152* .. Intrinsic Functions .. 153 INTRINSIC DBLE, ABS, MAX 154* .. 155* .. Executable Statements .. 156* 157* Test the input parameters. 158* 159 INFO = 0 160* 161 IF( NL.LT.1 ) THEN 162 INFO = -1 163 ELSE IF( NR.LT.1 ) THEN 164 INFO = -2 165 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 166 INFO = -3 167 END IF 168 IF( INFO.NE.0 ) THEN 169 CALL XERBLA( 'DLASD1', -INFO ) 170 RETURN 171 END IF 172* 173 N = NL + NR + 1 174 M = N + SQRE 175* 176* The following values are for bookkeeping purposes only. They are 177* integer pointers which indicate the portion of the workspace 178* used by a particular array in DLASD2 and DLASD3. 179* 180 LDU2 = N 181 LDVT2 = M 182* 183 IZ = 1 184 ISIGMA = IZ + M 185 IU2 = ISIGMA + N 186 IVT2 = IU2 + LDU2*N 187 IQ = IVT2 + LDVT2*M 188* 189 IDX = 1 190 IDXC = IDX + N 191 COLTYP = IDXC + N 192 IDXP = COLTYP + N 193* 194* Scale. 195* 196 ORGNRM = MAX( ABS( ALPHA ), ABS( BETA ) ) 197 D( NL+1 ) = ZERO 198 DO 10 I = 1, N 199 IF( ABS( D( I ) ).GT.ORGNRM ) THEN 200 ORGNRM = ABS( D( I ) ) 201 END IF 202 10 CONTINUE 203 OPS = OPS + DBLE( N + 2 ) 204 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO ) 205 ALPHA = ALPHA / ORGNRM 206 BETA = BETA / ORGNRM 207* 208* Deflate singular values. 209* 210 CALL DLASD2( NL, NR, SQRE, K, D, WORK( IZ ), ALPHA, BETA, U, LDU, 211 $ VT, LDVT, WORK( ISIGMA ), WORK( IU2 ), LDU2, 212 $ WORK( IVT2 ), LDVT2, IWORK( IDXP ), IWORK( IDX ), 213 $ IWORK( IDXC ), IDXQ, IWORK( COLTYP ), INFO ) 214* 215* Solve Secular Equation and update singular vectors. 216* 217 LDQ = K 218 CALL DLASD3( NL, NR, SQRE, K, D, WORK( IQ ), LDQ, WORK( ISIGMA ), 219 $ U, LDU, WORK( IU2 ), LDU2, VT, LDVT, WORK( IVT2 ), 220 $ LDVT2, IWORK( IDXC ), IWORK( COLTYP ), WORK( IZ ), 221 $ INFO ) 222 IF( INFO.NE.0 ) THEN 223 RETURN 224 END IF 225* 226* Unscale. 227* 228 OPS = OPS + DBLE( N ) 229 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO ) 230* 231* Prepare the IDXQ sorting permutation. 232* 233 N1 = K 234 N2 = N - K 235 CALL DLAMRG( N1, N2, D, 1, -1, IDXQ ) 236* 237 RETURN 238* 239* End of DLASD1 240* 241 END 242