|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, 00002 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 00003 $ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO ) 00004 * 00005 * -- LAPACK routine (version 3.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 * November 2006 00009 * 00010 * .. Scalar Arguments .. 00011 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL, 00012 $ LDGNUM, NL, NR, NRHS, SQRE 00013 DOUBLE PRECISION C, S 00014 * .. 00015 * .. Array Arguments .. 00016 INTEGER GIVCOL( LDGCOL, * ), PERM( * ) 00017 DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ), 00018 $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ), 00019 $ POLES( LDGNUM, * ), WORK( * ), Z( * ) 00020 * .. 00021 * 00022 * Purpose 00023 * ======= 00024 * 00025 * DLALS0 applies back the multiplying factors of either the left or the 00026 * right singular vector matrix of a diagonal matrix appended by a row 00027 * to the right hand side matrix B in solving the least squares problem 00028 * using the divide-and-conquer SVD approach. 00029 * 00030 * For the left singular vector matrix, three types of orthogonal 00031 * matrices are involved: 00032 * 00033 * (1L) Givens rotations: the number of such rotations is GIVPTR; the 00034 * pairs of columns/rows they were applied to are stored in GIVCOL; 00035 * and the C- and S-values of these rotations are stored in GIVNUM. 00036 * 00037 * (2L) Permutation. The (NL+1)-st row of B is to be moved to the first 00038 * row, and for J=2:N, PERM(J)-th row of B is to be moved to the 00039 * J-th row. 00040 * 00041 * (3L) The left singular vector matrix of the remaining matrix. 00042 * 00043 * For the right singular vector matrix, four types of orthogonal 00044 * matrices are involved: 00045 * 00046 * (1R) The right singular vector matrix of the remaining matrix. 00047 * 00048 * (2R) If SQRE = 1, one extra Givens rotation to generate the right 00049 * null space. 00050 * 00051 * (3R) The inverse transformation of (2L). 00052 * 00053 * (4R) The inverse transformation of (1L). 00054 * 00055 * Arguments 00056 * ========= 00057 * 00058 * ICOMPQ (input) INTEGER 00059 * Specifies whether singular vectors are to be computed in 00060 * factored form: 00061 * = 0: Left singular vector matrix. 00062 * = 1: Right singular vector matrix. 00063 * 00064 * NL (input) INTEGER 00065 * The row dimension of the upper block. NL >= 1. 00066 * 00067 * NR (input) INTEGER 00068 * The row dimension of the lower block. NR >= 1. 00069 * 00070 * SQRE (input) INTEGER 00071 * = 0: the lower block is an NR-by-NR square matrix. 00072 * = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00073 * 00074 * The bidiagonal matrix has row dimension N = NL + NR + 1, 00075 * and column dimension M = N + SQRE. 00076 * 00077 * NRHS (input) INTEGER 00078 * The number of columns of B and BX. NRHS must be at least 1. 00079 * 00080 * B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS ) 00081 * On input, B contains the right hand sides of the least 00082 * squares problem in rows 1 through M. On output, B contains 00083 * the solution X in rows 1 through N. 00084 * 00085 * LDB (input) INTEGER 00086 * The leading dimension of B. LDB must be at least 00087 * max(1,MAX( M, N ) ). 00088 * 00089 * BX (workspace) DOUBLE PRECISION array, dimension ( LDBX, NRHS ) 00090 * 00091 * LDBX (input) INTEGER 00092 * The leading dimension of BX. 00093 * 00094 * PERM (input) INTEGER array, dimension ( N ) 00095 * The permutations (from deflation and sorting) applied 00096 * to the two blocks. 00097 * 00098 * GIVPTR (input) INTEGER 00099 * The number of Givens rotations which took place in this 00100 * subproblem. 00101 * 00102 * GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 ) 00103 * Each pair of numbers indicates a pair of rows/columns 00104 * involved in a Givens rotation. 00105 * 00106 * LDGCOL (input) INTEGER 00107 * The leading dimension of GIVCOL, must be at least N. 00108 * 00109 * GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00110 * Each number indicates the C or S value used in the 00111 * corresponding Givens rotation. 00112 * 00113 * LDGNUM (input) INTEGER 00114 * The leading dimension of arrays DIFR, POLES and 00115 * GIVNUM, must be at least K. 00116 * 00117 * POLES (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00118 * On entry, POLES(1:K, 1) contains the new singular 00119 * values obtained from solving the secular equation, and 00120 * POLES(1:K, 2) is an array containing the poles in the secular 00121 * equation. 00122 * 00123 * DIFL (input) DOUBLE PRECISION array, dimension ( K ). 00124 * On entry, DIFL(I) is the distance between I-th updated 00125 * (undeflated) singular value and the I-th (undeflated) old 00126 * singular value. 00127 * 00128 * DIFR (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ). 00129 * On entry, DIFR(I, 1) contains the distances between I-th 00130 * updated (undeflated) singular value and the I+1-th 00131 * (undeflated) old singular value. And DIFR(I, 2) is the 00132 * normalizing factor for the I-th right singular vector. 00133 * 00134 * Z (input) DOUBLE PRECISION array, dimension ( K ) 00135 * Contain the components of the deflation-adjusted updating row 00136 * vector. 00137 * 00138 * K (input) INTEGER 00139 * Contains the dimension of the non-deflated matrix, 00140 * This is the order of the related secular equation. 1 <= K <=N. 00141 * 00142 * C (input) DOUBLE PRECISION 00143 * C contains garbage if SQRE =0 and the C-value of a Givens 00144 * rotation related to the right null space if SQRE = 1. 00145 * 00146 * S (input) DOUBLE PRECISION 00147 * S contains garbage if SQRE =0 and the S-value of a Givens 00148 * rotation related to the right null space if SQRE = 1. 00149 * 00150 * WORK (workspace) DOUBLE PRECISION array, dimension ( K ) 00151 * 00152 * INFO (output) INTEGER 00153 * = 0: successful exit. 00154 * < 0: if INFO = -i, the i-th argument had an illegal value. 00155 * 00156 * Further Details 00157 * =============== 00158 * 00159 * Based on contributions by 00160 * Ming Gu and Ren-Cang Li, Computer Science Division, University of 00161 * California at Berkeley, USA 00162 * Osni Marques, LBNL/NERSC, USA 00163 * 00164 * ===================================================================== 00165 * 00166 * .. Parameters .. 00167 DOUBLE PRECISION ONE, ZERO, NEGONE 00168 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 ) 00169 * .. 00170 * .. Local Scalars .. 00171 INTEGER I, J, M, N, NLP1 00172 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP 00173 * .. 00174 * .. External Subroutines .. 00175 EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DROT, DSCAL, 00176 $ XERBLA 00177 * .. 00178 * .. External Functions .. 00179 DOUBLE PRECISION DLAMC3, DNRM2 00180 EXTERNAL DLAMC3, DNRM2 00181 * .. 00182 * .. Intrinsic Functions .. 00183 INTRINSIC MAX 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( NL.LT.1 ) THEN 00194 INFO = -2 00195 ELSE IF( NR.LT.1 ) THEN 00196 INFO = -3 00197 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00198 INFO = -4 00199 END IF 00200 * 00201 N = NL + NR + 1 00202 * 00203 IF( NRHS.LT.1 ) THEN 00204 INFO = -5 00205 ELSE IF( LDB.LT.N ) THEN 00206 INFO = -7 00207 ELSE IF( LDBX.LT.N ) THEN 00208 INFO = -9 00209 ELSE IF( GIVPTR.LT.0 ) THEN 00210 INFO = -11 00211 ELSE IF( LDGCOL.LT.N ) THEN 00212 INFO = -13 00213 ELSE IF( LDGNUM.LT.N ) THEN 00214 INFO = -15 00215 ELSE IF( K.LT.1 ) THEN 00216 INFO = -20 00217 END IF 00218 IF( INFO.NE.0 ) THEN 00219 CALL XERBLA( 'DLALS0', -INFO ) 00220 RETURN 00221 END IF 00222 * 00223 M = N + SQRE 00224 NLP1 = NL + 1 00225 * 00226 IF( ICOMPQ.EQ.0 ) THEN 00227 * 00228 * Apply back orthogonal transformations from the left. 00229 * 00230 * Step (1L): apply back the Givens rotations performed. 00231 * 00232 DO 10 I = 1, GIVPTR 00233 CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, 00234 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), 00235 $ GIVNUM( I, 1 ) ) 00236 10 CONTINUE 00237 * 00238 * Step (2L): permute rows of B. 00239 * 00240 CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX ) 00241 DO 20 I = 2, N 00242 CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX ) 00243 20 CONTINUE 00244 * 00245 * Step (3L): apply the inverse of the left singular vector 00246 * matrix to BX. 00247 * 00248 IF( K.EQ.1 ) THEN 00249 CALL DCOPY( NRHS, BX, LDBX, B, LDB ) 00250 IF( Z( 1 ).LT.ZERO ) THEN 00251 CALL DSCAL( NRHS, NEGONE, B, LDB ) 00252 END IF 00253 ELSE 00254 DO 50 J = 1, K 00255 DIFLJ = DIFL( J ) 00256 DJ = POLES( J, 1 ) 00257 DSIGJ = -POLES( J, 2 ) 00258 IF( J.LT.K ) THEN 00259 DIFRJ = -DIFR( J, 1 ) 00260 DSIGJP = -POLES( J+1, 2 ) 00261 END IF 00262 IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) ) 00263 $ THEN 00264 WORK( J ) = ZERO 00265 ELSE 00266 WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ / 00267 $ ( POLES( J, 2 )+DJ ) 00268 END IF 00269 DO 30 I = 1, J - 1 00270 IF( ( Z( I ).EQ.ZERO ) .OR. 00271 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN 00272 WORK( I ) = ZERO 00273 ELSE 00274 WORK( I ) = POLES( I, 2 )*Z( I ) / 00275 $ ( DLAMC3( POLES( I, 2 ), DSIGJ )- 00276 $ DIFLJ ) / ( POLES( I, 2 )+DJ ) 00277 END IF 00278 30 CONTINUE 00279 DO 40 I = J + 1, K 00280 IF( ( Z( I ).EQ.ZERO ) .OR. 00281 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN 00282 WORK( I ) = ZERO 00283 ELSE 00284 WORK( I ) = POLES( I, 2 )*Z( I ) / 00285 $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+ 00286 $ DIFRJ ) / ( POLES( I, 2 )+DJ ) 00287 END IF 00288 40 CONTINUE 00289 WORK( 1 ) = NEGONE 00290 TEMP = DNRM2( K, WORK, 1 ) 00291 CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO, 00292 $ B( J, 1 ), LDB ) 00293 CALL DLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ), 00294 $ LDB, INFO ) 00295 50 CONTINUE 00296 END IF 00297 * 00298 * Move the deflated rows of BX to B also. 00299 * 00300 IF( K.LT.MAX( M, N ) ) 00301 $ CALL DLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX, 00302 $ B( K+1, 1 ), LDB ) 00303 ELSE 00304 * 00305 * Apply back the right orthogonal transformations. 00306 * 00307 * Step (1R): apply back the new right singular vector matrix 00308 * to B. 00309 * 00310 IF( K.EQ.1 ) THEN 00311 CALL DCOPY( NRHS, B, LDB, BX, LDBX ) 00312 ELSE 00313 DO 80 J = 1, K 00314 DSIGJ = POLES( J, 2 ) 00315 IF( Z( J ).EQ.ZERO ) THEN 00316 WORK( J ) = ZERO 00317 ELSE 00318 WORK( J ) = -Z( J ) / DIFL( J ) / 00319 $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 ) 00320 END IF 00321 DO 60 I = 1, J - 1 00322 IF( Z( J ).EQ.ZERO ) THEN 00323 WORK( I ) = ZERO 00324 ELSE 00325 WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1, 00326 $ 2 ) )-DIFR( I, 1 ) ) / 00327 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) 00328 END IF 00329 60 CONTINUE 00330 DO 70 I = J + 1, K 00331 IF( Z( J ).EQ.ZERO ) THEN 00332 WORK( I ) = ZERO 00333 ELSE 00334 WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I, 00335 $ 2 ) )-DIFL( I ) ) / 00336 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) 00337 END IF 00338 70 CONTINUE 00339 CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO, 00340 $ BX( J, 1 ), LDBX ) 00341 80 CONTINUE 00342 END IF 00343 * 00344 * Step (2R): if SQRE = 1, apply back the rotation that is 00345 * related to the right null space of the subproblem. 00346 * 00347 IF( SQRE.EQ.1 ) THEN 00348 CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX ) 00349 CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S ) 00350 END IF 00351 IF( K.LT.MAX( M, N ) ) 00352 $ CALL DLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ), 00353 $ LDBX ) 00354 * 00355 * Step (3R): permute rows of B. 00356 * 00357 CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB ) 00358 IF( SQRE.EQ.1 ) THEN 00359 CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB ) 00360 END IF 00361 DO 90 I = 2, N 00362 CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB ) 00363 90 CONTINUE 00364 * 00365 * Step (4R): apply back the Givens rotations performed. 00366 * 00367 DO 100 I = GIVPTR, 1, -1 00368 CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, 00369 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), 00370 $ -GIVNUM( I, 1 ) ) 00371 100 CONTINUE 00372 END IF 00373 * 00374 RETURN 00375 * 00376 * End of DLALS0 00377 * 00378 END
1.7.3 
