 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE SLASDA( 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 REAL 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, SLASDA 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, SLASD0, 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) REAL 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) REAL array, dimension ( M-1 ) 00063 * Contains the subdiagonal entries of the bidiagonal matrix. 00064 * On exit, E has been destroyed. 00065 * 00066 * U (output) REAL 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) REAL 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, dimension ( N ) 00083 * 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) REAL array, dimension ( LDU, NLVL ), 00088 * where NLVL = floor(log_2 (N/SMLSIZ))). 00089 * 00090 * DIFR (output) REAL 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 SLASD8 for details. 00098 * 00099 * Z (output) REAL 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) REAL 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, dimension ( LDGCOL, NLVL ) 00129 * 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) REAL 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) REAL 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) REAL 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) REAL array, dimension 00153 * (6 * N + (SMLSIZ + 1)*(SMLSIZ + 1)). 00154 * 00155 * IWORK (workspace) INTEGER array, dimension (7*N). 00156 * 00157 * INFO (output) INTEGER 00158 * = 0: successful exit. 00159 * < 0: if INFO = -i, the i-th argument had an illegal value. 00160 * > 0: if INFO = 1, a singular value did not converge 00161 * 00162 * Further Details 00163 * =============== 00164 * 00165 * Based on contributions by 00166 * Ming Gu and Huan Ren, Computer Science Division, University of 00167 * California at Berkeley, USA 00168 * 00169 * ===================================================================== 00170 * 00171 * .. Parameters .. 00172 REAL ZERO, ONE 00173 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0 ) 00174 * .. 00175 * .. Local Scalars .. 00176 INTEGER I, I1, IC, IDXQ, IDXQI, IM1, INODE, ITEMP, IWK, 00177 $ J, LF, LL, LVL, LVL2, M, NCC, ND, NDB1, NDIML, 00178 $ NDIMR, NL, NLF, NLP1, NLVL, NR, NRF, NRP1, NRU, 00179 $ NWORK1, NWORK2, SMLSZP, SQREI, VF, VFI, VL, VLI 00180 REAL ALPHA, BETA 00181 * .. 00182 * .. External Subroutines .. 00183 EXTERNAL SCOPY, SLASD6, SLASDQ, SLASDT, SLASET, XERBLA 00184 * .. 00185 * .. Executable Statements .. 00186 * 00187 * Test the input parameters. 00188 * 00189 INFO = 0 00190 * 00191 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00192 INFO = -1 00193 ELSE IF( SMLSIZ.LT.3 ) THEN 00194 INFO = -2 00195 ELSE IF( N.LT.0 ) THEN 00196 INFO = -3 00197 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00198 INFO = -4 00199 ELSE IF( LDU.LT.( N+SQRE ) ) THEN 00200 INFO = -8 00201 ELSE IF( LDGCOL.LT.N ) THEN 00202 INFO = -17 00203 END IF 00204 IF( INFO.NE.0 ) THEN 00205 CALL XERBLA( 'SLASDA', -INFO ) 00206 RETURN 00207 END IF 00208 * 00209 M = N + SQRE 00210 * 00211 * If the input matrix is too small, call SLASDQ to find the SVD. 00212 * 00213 IF( N.LE.SMLSIZ ) THEN 00214 IF( ICOMPQ.EQ.0 ) THEN 00215 CALL SLASDQ( 'U', SQRE, N, 0, 0, 0, D, E, VT, LDU, U, LDU, 00216 $ U, LDU, WORK, INFO ) 00217 ELSE 00218 CALL SLASDQ( 'U', SQRE, N, M, N, 0, D, E, VT, LDU, U, LDU, 00219 $ U, LDU, WORK, INFO ) 00220 END IF 00221 RETURN 00222 END IF 00223 * 00224 * Book-keeping and set up the computation tree. 00225 * 00226 INODE = 1 00227 NDIML = INODE + N 00228 NDIMR = NDIML + N 00229 IDXQ = NDIMR + N 00230 IWK = IDXQ + N 00231 * 00232 NCC = 0 00233 NRU = 0 00234 * 00235 SMLSZP = SMLSIZ + 1 00236 VF = 1 00237 VL = VF + M 00238 NWORK1 = VL + M 00239 NWORK2 = NWORK1 + SMLSZP*SMLSZP 00240 * 00241 CALL SLASDT( N, NLVL, ND, IWORK( INODE ), IWORK( NDIML ), 00242 $ IWORK( NDIMR ), SMLSIZ ) 00243 * 00244 * for the nodes on bottom level of the tree, solve 00245 * their subproblems by SLASDQ. 00246 * 00247 NDB1 = ( ND+1 ) / 2 00248 DO 30 I = NDB1, ND 00249 * 00250 * IC : center row of each node 00251 * NL : number of rows of left subproblem 00252 * NR : number of rows of right subproblem 00253 * NLF: starting row of the left subproblem 00254 * NRF: starting row of the right subproblem 00255 * 00256 I1 = I - 1 00257 IC = IWORK( INODE+I1 ) 00258 NL = IWORK( NDIML+I1 ) 00259 NLP1 = NL + 1 00260 NR = IWORK( NDIMR+I1 ) 00261 NLF = IC - NL 00262 NRF = IC + 1 00263 IDXQI = IDXQ + NLF - 2 00264 VFI = VF + NLF - 1 00265 VLI = VL + NLF - 1 00266 SQREI = 1 00267 IF( ICOMPQ.EQ.0 ) THEN 00268 CALL SLASET( 'A', NLP1, NLP1, ZERO, ONE, WORK( NWORK1 ), 00269 $ SMLSZP ) 00270 CALL SLASDQ( 'U', SQREI, NL, NLP1, NRU, NCC, D( NLF ), 00271 $ E( NLF ), WORK( NWORK1 ), SMLSZP, 00272 $ WORK( NWORK2 ), NL, WORK( NWORK2 ), NL, 00273 $ WORK( NWORK2 ), INFO ) 00274 ITEMP = NWORK1 + NL*SMLSZP 00275 CALL SCOPY( NLP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 ) 00276 CALL SCOPY( NLP1, WORK( ITEMP ), 1, WORK( VLI ), 1 ) 00277 ELSE 00278 CALL SLASET( 'A', NL, NL, ZERO, ONE, U( NLF, 1 ), LDU ) 00279 CALL SLASET( 'A', NLP1, NLP1, ZERO, ONE, VT( NLF, 1 ), LDU ) 00280 CALL SLASDQ( 'U', SQREI, NL, NLP1, NL, NCC, D( NLF ), 00281 $ E( NLF ), VT( NLF, 1 ), LDU, U( NLF, 1 ), LDU, 00282 $ U( NLF, 1 ), LDU, WORK( NWORK1 ), INFO ) 00283 CALL SCOPY( NLP1, VT( NLF, 1 ), 1, WORK( VFI ), 1 ) 00284 CALL SCOPY( NLP1, VT( NLF, NLP1 ), 1, WORK( VLI ), 1 ) 00285 END IF 00286 IF( INFO.NE.0 ) THEN 00287 RETURN 00288 END IF 00289 DO 10 J = 1, NL 00290 IWORK( IDXQI+J ) = J 00291 10 CONTINUE 00292 IF( ( I.EQ.ND ) .AND. ( SQRE.EQ.0 ) ) THEN 00293 SQREI = 0 00294 ELSE 00295 SQREI = 1 00296 END IF 00297 IDXQI = IDXQI + NLP1 00298 VFI = VFI + NLP1 00299 VLI = VLI + NLP1 00300 NRP1 = NR + SQREI 00301 IF( ICOMPQ.EQ.0 ) THEN 00302 CALL SLASET( 'A', NRP1, NRP1, ZERO, ONE, WORK( NWORK1 ), 00303 $ SMLSZP ) 00304 CALL SLASDQ( 'U', SQREI, NR, NRP1, NRU, NCC, D( NRF ), 00305 $ E( NRF ), WORK( NWORK1 ), SMLSZP, 00306 $ WORK( NWORK2 ), NR, WORK( NWORK2 ), NR, 00307 $ WORK( NWORK2 ), INFO ) 00308 ITEMP = NWORK1 + ( NRP1-1 )*SMLSZP 00309 CALL SCOPY( NRP1, WORK( NWORK1 ), 1, WORK( VFI ), 1 ) 00310 CALL SCOPY( NRP1, WORK( ITEMP ), 1, WORK( VLI ), 1 ) 00311 ELSE 00312 CALL SLASET( 'A', NR, NR, ZERO, ONE, U( NRF, 1 ), LDU ) 00313 CALL SLASET( 'A', NRP1, NRP1, ZERO, ONE, VT( NRF, 1 ), LDU ) 00314 CALL SLASDQ( 'U', SQREI, NR, NRP1, NR, NCC, D( NRF ), 00315 $ E( NRF ), VT( NRF, 1 ), LDU, U( NRF, 1 ), LDU, 00316 $ U( NRF, 1 ), LDU, WORK( NWORK1 ), INFO ) 00317 CALL SCOPY( NRP1, VT( NRF, 1 ), 1, WORK( VFI ), 1 ) 00318 CALL SCOPY( NRP1, VT( NRF, NRP1 ), 1, WORK( VLI ), 1 ) 00319 END IF 00320 IF( INFO.NE.0 ) THEN 00321 RETURN 00322 END IF 00323 DO 20 J = 1, NR 00324 IWORK( IDXQI+J ) = J 00325 20 CONTINUE 00326 30 CONTINUE 00327 * 00328 * Now conquer each subproblem bottom-up. 00329 * 00330 J = 2**NLVL 00331 DO 50 LVL = NLVL, 1, -1 00332 LVL2 = LVL*2 - 1 00333 * 00334 * Find the first node LF and last node LL on 00335 * the current level LVL. 00336 * 00337 IF( LVL.EQ.1 ) THEN 00338 LF = 1 00339 LL = 1 00340 ELSE 00341 LF = 2**( LVL-1 ) 00342 LL = 2*LF - 1 00343 END IF 00344 DO 40 I = LF, LL 00345 IM1 = I - 1 00346 IC = IWORK( INODE+IM1 ) 00347 NL = IWORK( NDIML+IM1 ) 00348 NR = IWORK( NDIMR+IM1 ) 00349 NLF = IC - NL 00350 NRF = IC + 1 00351 IF( I.EQ.LL ) THEN 00352 SQREI = SQRE 00353 ELSE 00354 SQREI = 1 00355 END IF 00356 VFI = VF + NLF - 1 00357 VLI = VL + NLF - 1 00358 IDXQI = IDXQ + NLF - 1 00359 ALPHA = D( IC ) 00360 BETA = E( IC ) 00361 IF( ICOMPQ.EQ.0 ) THEN 00362 CALL SLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ), 00363 $ WORK( VFI ), WORK( VLI ), ALPHA, BETA, 00364 $ IWORK( IDXQI ), PERM, GIVPTR( 1 ), GIVCOL, 00365 $ LDGCOL, GIVNUM, LDU, POLES, DIFL, DIFR, Z, 00366 $ K( 1 ), C( 1 ), S( 1 ), WORK( NWORK1 ), 00367 $ IWORK( IWK ), INFO ) 00368 ELSE 00369 J = J - 1 00370 CALL SLASD6( ICOMPQ, NL, NR, SQREI, D( NLF ), 00371 $ WORK( VFI ), WORK( VLI ), ALPHA, BETA, 00372 $ IWORK( IDXQI ), PERM( NLF, LVL ), 00373 $ GIVPTR( J ), GIVCOL( NLF, LVL2 ), LDGCOL, 00374 $ GIVNUM( NLF, LVL2 ), LDU, 00375 $ POLES( NLF, LVL2 ), DIFL( NLF, LVL ), 00376 $ DIFR( NLF, LVL2 ), Z( NLF, LVL ), K( J ), 00377 $ C( J ), S( J ), WORK( NWORK1 ), 00378 $ IWORK( IWK ), INFO ) 00379 END IF 00380 IF( INFO.NE.0 ) THEN 00381 RETURN 00382 END IF 00383 40 CONTINUE 00384 50 CONTINUE 00385 * 00386 RETURN 00387 * 00388 * End of SLASDA 00389 * 00390 END
1.7.3 
