 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DLASDA( ICOMPQ, SMLSIZ, N, SQRE, D, E, U, LDU, VT, K, 00002 $ DIFL, DIFR, Z, POLES, GIVPTR, GIVCOL, LDGCOL, 00003 $ PERM, GIVNUM, C, S, WORK, IWORK, INFO ) 00004 * 00005 * -- LAPACK auxiliary routine (version 3.2.2) -- 00006 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00007 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00008 * June 2010 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER ICOMPQ, INFO, LDGCOL, LDU, N, SMLSIZ, SQRE 00012 * .. 00013 * .. Array Arguments .. 00014 INTEGER GIVCOL( LDGCOL, * ), GIVPTR( * ), IWORK( * ), 00015 $ K( * ), PERM( LDGCOL, * ) 00016 DOUBLE PRECISION C( * ), D( * ), DIFL( LDU, * ), DIFR( LDU, * ), 00017 $ E( * ), GIVNUM( LDU, * ), POLES( LDU, * ), 00018 $ S( * ), U( LDU, * ), VT( LDU, * ), WORK( * ), 00019 $ Z( LDU, * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * Using a divide and conquer approach, DLASDA computes the singular 00026 * value decomposition (SVD) of a real upper bidiagonal N-by-M matrix 00027 * B with diagonal D and offdiagonal E, where M = N + SQRE. The 00028 * algorithm computes the singular values in the SVD B = U * S * VT. 00029 * The orthogonal matrices U and VT are optionally computed in 00030 * compact form. 00031 * 00032 * A related subroutine, DLASD0, computes the singular values and 00033 * the singular vectors in explicit form. 00034 * 00035 * Arguments 00036 * ========= 00037 * 00038 * ICOMPQ (input) INTEGER 00039 * Specifies whether singular vectors are to be computed 00040 * in compact form, as follows 00041 * = 0: Compute singular values only. 00042 * = 1: Compute singular vectors of upper bidiagonal 00043 * matrix in compact form. 00044 * 00045 * SMLSIZ (input) INTEGER 00046 * The maximum size of the subproblems at the bottom of the 00047 * computation tree. 00048 * 00049 * N (input) INTEGER 00050 * The row dimension of the upper bidiagonal matrix. This is 00051 * also the dimension of the main diagonal array D. 00052 * 00053 * SQRE (input) INTEGER 00054 * Specifies the column dimension of the bidiagonal matrix. 00055 * = 0: The bidiagonal matrix has column dimension M = N; 00056 * = 1: The bidiagonal matrix has column dimension M = N + 1. 00057 * 00058 * D (input/output) DOUBLE PRECISION array, dimension ( N ) 00059 * On entry D contains the main diagonal of the bidiagonal 00060 * matrix. On exit D, if INFO = 0, contains its singular values. 00061 * 00062 * E (input) DOUBLE PRECISION array, dimension ( M-1 ) 00063 * Contains the subdiagonal entries of the bidiagonal matrix. 00064 * On exit, E has been destroyed. 00065 * 00066 * U (output) DOUBLE PRECISION array, 00067 * dimension ( LDU, SMLSIZ ) if ICOMPQ = 1, and not referenced 00068 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, U contains the left 00069 * singular vector matrices of all subproblems at the bottom 00070 * level. 00071 * 00072 * LDU (input) INTEGER, LDU = > N. 00073 * The leading dimension of arrays U, VT, DIFL, DIFR, POLES, 00074 * GIVNUM, and Z. 00075 * 00076 * VT (output) DOUBLE PRECISION array, 00077 * dimension ( LDU, SMLSIZ+1 ) if ICOMPQ = 1, and not referenced 00078 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, VT**T contains the right 00079 * singular vector matrices of all subproblems at the bottom 00080 * level. 00081 * 00082 * K (output) INTEGER array, 00083 * dimension ( N ) if ICOMPQ = 1 and dimension 1 if ICOMPQ = 0. 00084 * If ICOMPQ = 1, on exit, K(I) is the dimension of the I-th 00085 * secular equation on the computation tree. 00086 * 00087 * DIFL (output) DOUBLE PRECISION array, dimension ( LDU, NLVL ), 00088 * where NLVL = floor(log_2 (N/SMLSIZ))). 00089 * 00090 * DIFR (output) DOUBLE PRECISION array, 00091 * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1 and 00092 * dimension ( N ) if ICOMPQ = 0. 00093 * If ICOMPQ = 1, on exit, DIFL(1:N, I) and DIFR(1:N, 2 * I - 1) 00094 * record distances between singular values on the I-th 00095 * level and singular values on the (I -1)-th level, and 00096 * DIFR(1:N, 2 * I ) contains the normalizing factors for 00097 * the right singular vector matrix. See DLASD8 for details. 00098 * 00099 * Z (output) DOUBLE PRECISION array, 00100 * dimension ( LDU, NLVL ) if ICOMPQ = 1 and 00101 * dimension ( N ) if ICOMPQ = 0. 00102 * The first K elements of Z(1, I) contain the components of 00103 * the deflation-adjusted updating row vector for subproblems 00104 * on the I-th level. 00105 * 00106 * POLES (output) DOUBLE PRECISION array, 00107 * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not referenced 00108 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, POLES(1, 2*I - 1) and 00109 * POLES(1, 2*I) contain the new and old singular values 00110 * involved in the secular equations on the I-th level. 00111 * 00112 * GIVPTR (output) INTEGER array, 00113 * dimension ( N ) if ICOMPQ = 1, and not referenced if 00114 * ICOMPQ = 0. If ICOMPQ = 1, on exit, GIVPTR( I ) records 00115 * the number of Givens rotations performed on the I-th 00116 * problem on the computation tree. 00117 * 00118 * GIVCOL (output) INTEGER array, 00119 * dimension ( LDGCOL, 2 * NLVL ) if ICOMPQ = 1, and not 00120 * referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I, 00121 * GIVCOL(1, 2 *I - 1) and GIVCOL(1, 2 *I) record the locations 00122 * of Givens rotations performed on the I-th level on the 00123 * computation tree. 00124 * 00125 * LDGCOL (input) INTEGER, LDGCOL = > N. 00126 * The leading dimension of arrays GIVCOL and PERM. 00127 * 00128 * PERM (output) INTEGER array, 00129 * dimension ( LDGCOL, NLVL ) if ICOMPQ = 1, and not referenced 00130 * if ICOMPQ = 0. If ICOMPQ = 1, on exit, PERM(1, I) records 00131 * permutations done on the I-th level of the computation tree. 00132 * 00133 * GIVNUM (output) DOUBLE PRECISION array, 00134 * dimension ( LDU, 2 * NLVL ) if ICOMPQ = 1, and not 00135 * referenced if ICOMPQ = 0. If ICOMPQ = 1, on exit, for each I, 00136 * GIVNUM(1, 2 *I - 1) and GIVNUM(1, 2 *I) record the C- and S- 00137 * values of Givens rotations performed on the I-th level on 00138 * the computation tree. 00139 * 00140 * C (output) DOUBLE PRECISION array, 00141 * dimension ( N ) if ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. 00142 * If ICOMPQ = 1 and the I-th subproblem is not square, on exit, 00143 * C( I ) contains the C-value of a Givens rotation related to 00144 * the right null space of the I-th subproblem. 00145 * 00146 * S (output) DOUBLE PRECISION array, dimension ( N ) if 00147 * ICOMPQ = 1, and dimension 1 if ICOMPQ = 0. If ICOMPQ = 1 00148 * and the I-th subproblem is not square, on exit, S( I ) 00149 * contains the S-value of a Givens rotation related to 00150 * the right null space of the I-th subproblem. 00151 * 00152 * WORK (workspace) DOUBLE PRECISION array, dimension 00153 * (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)). 00154 * 00155 * IWORK (workspace) INTEGER array. 00156 * Dimension must be at least (7 * N). 00157 * 00158 * INFO (output) INTEGER 00159 * = 0: successful exit. 00160 * < 0: if INFO = -i, the i-th argument had an illegal value. 00161 * > 0: if INFO = 1, a singular value did not converge 00162 * 00163 * Further Details 00164 * =============== 00165 * 00166 * Based on contributions by 00167 * Ming Gu and Huan Ren, Computer Science Division, University of 00168 * California at Berkeley, USA 00169 * 00170 * ===================================================================== 00171 * 00172 * .. Parameters .. 00173 DOUBLE PRECISION ZERO, ONE 00174 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) 00175 * .. 00176 * .. Local Scalars .. 00177 INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK, 00178 $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML, 00179 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU, 00180 $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI 00181 DOUBLE PRECISION ALPHA, BETA 00182 * .. 00183 * .. External Subroutines .. 00184 EXTERNAL DCOPY, DLASD6, DLASDQ, DLASDT, DLASET, XERBLA 00185 * .. 00186 * .. Executable Statements .. 00187 * 00188 * Test the input parameters. 00189 * 00190 INFO = 0 00191 * 00192 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00193 INFO = -1 00194 ELSE IF( SMLSIZ.LT.3 ) THEN 00195 INFO = -2 00196 ELSE IF( N.LT.0 ) THEN 00197 INFO = -3 00198 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00199 INFO = -4 00200 ELSE IF( LDU.LT.( N+SQRE ) ) THEN 00201 INFO = -8 00202 ELSE IF( LDGCOL.LT.N ) THEN 00203 INFO = -17 00204 END IF 00205 IF( INFO.NE.0 ) THEN 00206 CALL XERBLA( 'DLASDA', -INFO ) 00207 RETURN 00208 END IF 00209 * 00210 M = N + SQRE 00211 * 00212 * If the input matrix is too small, call DLASDQ to find the SVD. 00213 * 00214 IF( N.LE.SMLSIZ ) THEN 00215 IF( ICOMPQ.EQ.0 ) THEN 00216 CALL DLASDQ( 'U', SQRE, N, 0, 0, 0, D, E, VT, LDU, U, LDU, 00217 $ U, LDU, WORK, INFO ) 00218 ELSE 00219 CALL DLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDU, U, LDU, 00220 $ U, LDU, WORK, INFO ) 00221 END IF 00222 RETURN 00223 END IF 00224 * 00225 * Book-keeping and set up the computation tree. 00226 * 00227 INODE = 1 00228 NDIML = INODE + N 00229 NDIMR = NDIML + N 00230 IDXQ = NDIMR + N 00231 IWK = IDXQ + N 00232 * 00233 NCC = 0 00234 NRU = 0 00235 * 00236 SMLSZP = SMLSIZ + 1 00237 VF = 1 00238 VL = VF + M 00239 NWORK1 = VL + M 00240 NWORK2 = NWORK1 + SMLSZP*SMLSZP 00241 * 00242 CALL DLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00243 $ IWORK( NDIMR ), SMLSIZ ) 00244 * 00245 * for the nodes on bottom level of the tree, solve 00246 * their subproblems by DLASDQ. 00247 * 00248 NDB1 = ( ND+1 ) / 2 00249 DO 30 I = NDB1, ND 00250 * 00251 * IC : center row of each node 00252 * NL : number of rows of left subproblem 00253 * NR : number of rows of right subproblem 00254 * NLF: starting row of the left subproblem 00255 * NRF: starting row of the right subproblem 00256 * 00257 I1 = I - 1 00258 IC = IWORK( INODE+I1 ) 00259 NL = IWORK( NDIML+I1 ) 00260 NLP1 = NL + 1 00261 NR = IWORK( NDIMR+I1 ) 00262 NLF = IC - NL 00263 NRF = IC + 1 00264 IDXQI = IDXQ + NLF - 2 00265 VFI = VF + NLF - 1 00266 VLI = VL + NLF - 1 00267 SQREI = 1 00268 IF( ICOMPQ.EQ.0 ) THEN 00269 CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, WORK( NWORK1 ), 00270 $ SMLSZP ) 00271 CALL DLASDQ( 'U', SQREI, NL, NLP1, NRU, NCC, D( NLF ), 00272 $ E( NLF ), WORK( NWORK1 ), SMLSZP, 00273 $ WORK( NWORK2 ), NL, WORK( NWORK2 ), NL, 00274 $ WORK( NWORK2 ), INFO ) 00275 ITEMP = NWORK1 + NL*SMLSZP 00276 CALL DCOPY( NLP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 ) 00277 CALL DCOPY( NLP1, WORK( ITEMP ), 1, WORK( VLI ), 1 ) 00278 ELSE 00279 CALL DLASET( 'A', NL, NL, ZERO, ONE, U( NLF, 1 ), LDU ) 00280 CALL DLASET( 'A', NLP1, NLP1, ZERO, ONE, VT( NLF, 1 ), LDU ) 00281 CALL DLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), 00282 $ E( NLF ), VT( NLF, 1 ), LDU, U( NLF, 1 ), LDU, 00283 $ U( NLF, 1 ), LDU, WORK( NWORK1 ), INFO ) 00284 CALL DCOPY( NLP1, VT( NLF, 1 ), 1, WORK( VFI ), 1 ) 00285 CALL DCOPY( NLP1, VT( NLF, NLP1 ), 1, WORK( VLI ), 1 ) 00286 END IF 00287 IF( INFO.NE.0 ) THEN 00288 RETURN 00289 END IF 00290 DO 10 J = 1, NL 00291 IWORK( IDXQI+J ) = J 00292 10 CONTINUE 00293 IF( ( I.EQ.ND ) .AND. ( SQRE.EQ.0 ) ) THEN 00294 SQREI = 0 00295 ELSE 00296 SQREI = 1 00297 END IF 00298 IDXQI = IDXQI + NLP1 00299 VFI = VFI + NLP1 00300 VLI = VLI + NLP1 00301 NRP1 = NR + SQREI 00302 IF( ICOMPQ.EQ.0 ) THEN 00303 CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, WORK( NWORK1 ), 00304 $ SMLSZP ) 00305 CALL DLASDQ( 'U', SQREI, NR, NRP1, NRU, NCC, D( NRF ), 00306 $ E( NRF ), WORK( NWORK1 ), SMLSZP, 00307 $ WORK( NWORK2 ), NR, WORK( NWORK2 ), NR, 00308 $ WORK( NWORK2 ), INFO ) 00309 ITEMP = NWORK1 + ( NRP1-1 )*SMLSZP 00310 CALL DCOPY( NRP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 ) 00311 CALL DCOPY( NRP1, WORK( ITEMP ), 1, WORK( VLI ), 1 ) 00312 ELSE 00313 CALL DLASET( 'A', NR, NR, ZERO, ONE, U( NRF, 1 ), LDU ) 00314 CALL DLASET( 'A', NRP1, NRP1, ZERO, ONE, VT( NRF, 1 ), LDU ) 00315 CALL DLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), 00316 $ E( NRF ), VT( NRF, 1 ), LDU, U( NRF, 1 ), LDU, 00317 $ U( NRF, 1 ), LDU, WORK( NWORK1 ), INFO ) 00318 CALL DCOPY( NRP1, VT( NRF, 1 ), 1, WORK( VFI ), 1 ) 00319 CALL DCOPY( NRP1, VT( NRF, NRP1 ), 1, WORK( VLI ), 1 ) 00320 END IF 00321 IF( INFO.NE.0 ) THEN 00322 RETURN 00323 END IF 00324 DO 20 J = 1, NR 00325 IWORK( IDXQI+J ) = J 00326 20 CONTINUE 00327 30 CONTINUE 00328 * 00329 * Now conquer each subproblem bottom-up. 00330 * 00331 J = 2**NLVL 00332 DO 50 LVL = NLVL, 1, -1 00333 LVL2 = LVL*2 - 1 00334 * 00335 * Find the first node LF and last node LL on 00336 * the current level LVL. 00337 * 00338 IF( LVL.EQ.1 ) THEN 00339 LF = 1 00340 LL = 1 00341 ELSE 00342 LF = 2**( LVL-1 ) 00343 LL = 2*LF - 1 00344 END IF 00345 DO 40 I = LF, LL 00346 IM1 = I - 1 00347 IC = IWORK( INODE+IM1 ) 00348 NL = IWORK( NDIML+IM1 ) 00349 NR = IWORK( NDIMR+IM1 ) 00350 NLF = IC - NL 00351 NRF = IC + 1 00352 IF( I.EQ.LL ) THEN 00353 SQREI = SQRE 00354 ELSE 00355 SQREI = 1 00356 END IF 00357 VFI = VF + NLF - 1 00358 VLI = VL + NLF - 1 00359 IDXQI = IDXQ + NLF - 1 00360 ALPHA = D( IC ) 00361 BETA = E( IC ) 00362 IF( ICOMPQ.EQ.0 ) THEN 00363 CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ), 00364 $ WORK( VFI ), WORK( VLI ), ALPHA, BETA, 00365 $ IWORK( IDXQI ), PERM, GIVPTR( 1 ), GIVCOL, 00366 $ LDGCOL, GIVNUM, LDU, POLES, DIFL, DIFR, Z, 00367 $ K( 1 ), C( 1 ), S( 1 ), WORK( NWORK1 ), 00368 $ IWORK( IWK ), INFO ) 00369 ELSE 00370 J = J - 1 00371 CALL DLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ), 00372 $ WORK( VFI ), WORK( VLI ), ALPHA, BETA, 00373 $ IWORK( IDXQI ), PERM( NLF, LVL ), 00374 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00375 $ GIVNUM( NLF, LVL2 ), LDU, 00376 $ POLES( NLF, LVL2 ), DIFL( NLF, LVL ), 00377 $ DIFR( NLF, LVL2 ), Z( NLF, LVL ), K( J ), 00378 $ C( J ), S( J ), WORK( NWORK1 ), 00379 $ IWORK( IWK ), INFO ) 00380 END IF 00381 IF( INFO.NE.0 ) THEN 00382 RETURN 00383 END IF 00384 40 CONTINUE 00385 50 CONTINUE 00386 * 00387 RETURN 00388 * 00389 * End of DLASDA 00390 * 00391 END
1.7.3 
