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