![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLALS0 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLALS0 + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlals0.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlals0.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlals0.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, 00022 * PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 00023 * POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL, 00027 * $ LDGNUM, NL, NR, NRHS, SQRE 00028 * DOUBLE PRECISION C, S 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER GIVCOL( LDGCOL, * ), PERM( * ) 00032 * DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ), 00033 * $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ), 00034 * $ POLES( LDGNUM, * ), WORK( * ), Z( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> DLALS0 applies back the multiplying factors of either the left or the 00044 *> right singular vector matrix of a diagonal matrix appended by a row 00045 *> to the right hand side matrix B in solving the least squares problem 00046 *> using the divide-and-conquer SVD approach. 00047 *> 00048 *> For the left singular vector matrix, three types of orthogonal 00049 *> matrices are involved: 00050 *> 00051 *> (1L) Givens rotations: the number of such rotations is GIVPTR; the 00052 *> pairs of columns/rows they were applied to are stored in GIVCOL; 00053 *> and the C- and S-values of these rotations are stored in GIVNUM. 00054 *> 00055 *> (2L) Permutation. The (NL+1)-st row of B is to be moved to the first 00056 *> row, and for J=2:N, PERM(J)-th row of B is to be moved to the 00057 *> J-th row. 00058 *> 00059 *> (3L) The left singular vector matrix of the remaining matrix. 00060 *> 00061 *> For the right singular vector matrix, four types of orthogonal 00062 *> matrices are involved: 00063 *> 00064 *> (1R) The right singular vector matrix of the remaining matrix. 00065 *> 00066 *> (2R) If SQRE = 1, one extra Givens rotation to generate the right 00067 *> null space. 00068 *> 00069 *> (3R) The inverse transformation of (2L). 00070 *> 00071 *> (4R) The inverse transformation of (1L). 00072 *> \endverbatim 00073 * 00074 * Arguments: 00075 * ========== 00076 * 00077 *> \param[in] ICOMPQ 00078 *> \verbatim 00079 *> ICOMPQ is INTEGER 00080 *> Specifies whether singular vectors are to be computed in 00081 *> factored form: 00082 *> = 0: Left singular vector matrix. 00083 *> = 1: Right singular vector matrix. 00084 *> \endverbatim 00085 *> 00086 *> \param[in] NL 00087 *> \verbatim 00088 *> NL is INTEGER 00089 *> The row dimension of the upper block. NL >= 1. 00090 *> \endverbatim 00091 *> 00092 *> \param[in] NR 00093 *> \verbatim 00094 *> NR is INTEGER 00095 *> The row dimension of the lower block. NR >= 1. 00096 *> \endverbatim 00097 *> 00098 *> \param[in] SQRE 00099 *> \verbatim 00100 *> SQRE is INTEGER 00101 *> = 0: the lower block is an NR-by-NR square matrix. 00102 *> = 1: the lower block is an NR-by-(NR+1) rectangular matrix. 00103 *> 00104 *> The bidiagonal matrix has row dimension N = NL + NR + 1, 00105 *> and column dimension M = N + SQRE. 00106 *> \endverbatim 00107 *> 00108 *> \param[in] NRHS 00109 *> \verbatim 00110 *> NRHS is INTEGER 00111 *> The number of columns of B and BX. NRHS must be at least 1. 00112 *> \endverbatim 00113 *> 00114 *> \param[in,out] B 00115 *> \verbatim 00116 *> B is DOUBLE PRECISION array, dimension ( LDB, NRHS ) 00117 *> On input, B contains the right hand sides of the least 00118 *> squares problem in rows 1 through M. On output, B contains 00119 *> the solution X in rows 1 through N. 00120 *> \endverbatim 00121 *> 00122 *> \param[in] LDB 00123 *> \verbatim 00124 *> LDB is INTEGER 00125 *> The leading dimension of B. LDB must be at least 00126 *> max(1,MAX( M, N ) ). 00127 *> \endverbatim 00128 *> 00129 *> \param[out] BX 00130 *> \verbatim 00131 *> BX is DOUBLE PRECISION array, dimension ( LDBX, NRHS ) 00132 *> \endverbatim 00133 *> 00134 *> \param[in] LDBX 00135 *> \verbatim 00136 *> LDBX is INTEGER 00137 *> The leading dimension of BX. 00138 *> \endverbatim 00139 *> 00140 *> \param[in] PERM 00141 *> \verbatim 00142 *> PERM is INTEGER array, dimension ( N ) 00143 *> The permutations (from deflation and sorting) applied 00144 *> to the two blocks. 00145 *> \endverbatim 00146 *> 00147 *> \param[in] GIVPTR 00148 *> \verbatim 00149 *> GIVPTR is INTEGER 00150 *> The number of Givens rotations which took place in this 00151 *> subproblem. 00152 *> \endverbatim 00153 *> 00154 *> \param[in] GIVCOL 00155 *> \verbatim 00156 *> GIVCOL is INTEGER array, dimension ( LDGCOL, 2 ) 00157 *> Each pair of numbers indicates a pair of rows/columns 00158 *> involved in a Givens rotation. 00159 *> \endverbatim 00160 *> 00161 *> \param[in] LDGCOL 00162 *> \verbatim 00163 *> LDGCOL is INTEGER 00164 *> The leading dimension of GIVCOL, must be at least N. 00165 *> \endverbatim 00166 *> 00167 *> \param[in] GIVNUM 00168 *> \verbatim 00169 *> GIVNUM is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00170 *> Each number indicates the C or S value used in the 00171 *> corresponding Givens rotation. 00172 *> \endverbatim 00173 *> 00174 *> \param[in] LDGNUM 00175 *> \verbatim 00176 *> LDGNUM is INTEGER 00177 *> The leading dimension of arrays DIFR, POLES and 00178 *> GIVNUM, must be at least K. 00179 *> \endverbatim 00180 *> 00181 *> \param[in] POLES 00182 *> \verbatim 00183 *> POLES is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) 00184 *> On entry, POLES(1:K, 1) contains the new singular 00185 *> values obtained from solving the secular equation, and 00186 *> POLES(1:K, 2) is an array containing the poles in the secular 00187 *> equation. 00188 *> \endverbatim 00189 *> 00190 *> \param[in] DIFL 00191 *> \verbatim 00192 *> DIFL is DOUBLE PRECISION array, dimension ( K ). 00193 *> On entry, DIFL(I) is the distance between I-th updated 00194 *> (undeflated) singular value and the I-th (undeflated) old 00195 *> singular value. 00196 *> \endverbatim 00197 *> 00198 *> \param[in] DIFR 00199 *> \verbatim 00200 *> DIFR is DOUBLE PRECISION array, dimension ( LDGNUM, 2 ). 00201 *> On entry, DIFR(I, 1) contains the distances between I-th 00202 *> updated (undeflated) singular value and the I+1-th 00203 *> (undeflated) old singular value. And DIFR(I, 2) is the 00204 *> normalizing factor for the I-th right singular vector. 00205 *> \endverbatim 00206 *> 00207 *> \param[in] Z 00208 *> \verbatim 00209 *> Z is DOUBLE PRECISION array, dimension ( K ) 00210 *> Contain the components of the deflation-adjusted updating row 00211 *> vector. 00212 *> \endverbatim 00213 *> 00214 *> \param[in] K 00215 *> \verbatim 00216 *> K is INTEGER 00217 *> Contains the dimension of the non-deflated matrix, 00218 *> This is the order of the related secular equation. 1 <= K <=N. 00219 *> \endverbatim 00220 *> 00221 *> \param[in] C 00222 *> \verbatim 00223 *> C is DOUBLE PRECISION 00224 *> C contains garbage if SQRE =0 and the C-value of a Givens 00225 *> rotation related to the right null space if SQRE = 1. 00226 *> \endverbatim 00227 *> 00228 *> \param[in] S 00229 *> \verbatim 00230 *> S is DOUBLE PRECISION 00231 *> S contains garbage if SQRE =0 and the S-value of a Givens 00232 *> rotation related to the right null space if SQRE = 1. 00233 *> \endverbatim 00234 *> 00235 *> \param[out] WORK 00236 *> \verbatim 00237 *> WORK is DOUBLE PRECISION array, dimension ( K ) 00238 *> \endverbatim 00239 *> 00240 *> \param[out] INFO 00241 *> \verbatim 00242 *> INFO is INTEGER 00243 *> = 0: successful exit. 00244 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00245 *> \endverbatim 00246 * 00247 * Authors: 00248 * ======== 00249 * 00250 *> \author Univ. of Tennessee 00251 *> \author Univ. of California Berkeley 00252 *> \author Univ. of Colorado Denver 00253 *> \author NAG Ltd. 00254 * 00255 *> \date November 2011 00256 * 00257 *> \ingroup doubleOTHERcomputational 00258 * 00259 *> \par Contributors: 00260 * ================== 00261 *> 00262 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of 00263 *> California at Berkeley, USA \n 00264 *> Osni Marques, LBNL/NERSC, USA \n 00265 * 00266 * ===================================================================== 00267 SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, 00268 $ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, 00269 $ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO ) 00270 * 00271 * -- LAPACK computational routine (version 3.4.0) -- 00272 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00273 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00274 * November 2011 00275 * 00276 * .. Scalar Arguments .. 00277 INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL, 00278 $ LDGNUM, NL, NR, NRHS, SQRE 00279 DOUBLE PRECISION C, S 00280 * .. 00281 * .. Array Arguments .. 00282 INTEGER GIVCOL( LDGCOL, * ), PERM( * ) 00283 DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ), 00284 $ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ), 00285 $ POLES( LDGNUM, * ), WORK( * ), Z( * ) 00286 * .. 00287 * 00288 * ===================================================================== 00289 * 00290 * .. Parameters .. 00291 DOUBLE PRECISION ONE, ZERO, NEGONE 00292 PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 ) 00293 * .. 00294 * .. Local Scalars .. 00295 INTEGER I, J, M, N, NLP1 00296 DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP 00297 * .. 00298 * .. External Subroutines .. 00299 EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DROT, DSCAL, 00300 $ XERBLA 00301 * .. 00302 * .. External Functions .. 00303 DOUBLE PRECISION DLAMC3, DNRM2 00304 EXTERNAL DLAMC3, DNRM2 00305 * .. 00306 * .. Intrinsic Functions .. 00307 INTRINSIC MAX 00308 * .. 00309 * .. Executable Statements .. 00310 * 00311 * Test the input parameters. 00312 * 00313 INFO = 0 00314 * 00315 IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN 00316 INFO = -1 00317 ELSE IF( NL.LT.1 ) THEN 00318 INFO = -2 00319 ELSE IF( NR.LT.1 ) THEN 00320 INFO = -3 00321 ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN 00322 INFO = -4 00323 END IF 00324 * 00325 N = NL + NR + 1 00326 * 00327 IF( NRHS.LT.1 ) THEN 00328 INFO = -5 00329 ELSE IF( LDB.LT.N ) THEN 00330 INFO = -7 00331 ELSE IF( LDBX.LT.N ) THEN 00332 INFO = -9 00333 ELSE IF( GIVPTR.LT.0 ) THEN 00334 INFO = -11 00335 ELSE IF( LDGCOL.LT.N ) THEN 00336 INFO = -13 00337 ELSE IF( LDGNUM.LT.N ) THEN 00338 INFO = -15 00339 ELSE IF( K.LT.1 ) THEN 00340 INFO = -20 00341 END IF 00342 IF( INFO.NE.0 ) THEN 00343 CALL XERBLA( 'DLALS0', -INFO ) 00344 RETURN 00345 END IF 00346 * 00347 M = N + SQRE 00348 NLP1 = NL + 1 00349 * 00350 IF( ICOMPQ.EQ.0 ) THEN 00351 * 00352 * Apply back orthogonal transformations from the left. 00353 * 00354 * Step (1L): apply back the Givens rotations performed. 00355 * 00356 DO 10 I = 1, GIVPTR 00357 CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, 00358 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), 00359 $ GIVNUM( I, 1 ) ) 00360 10 CONTINUE 00361 * 00362 * Step (2L): permute rows of B. 00363 * 00364 CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX ) 00365 DO 20 I = 2, N 00366 CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX ) 00367 20 CONTINUE 00368 * 00369 * Step (3L): apply the inverse of the left singular vector 00370 * matrix to BX. 00371 * 00372 IF( K.EQ.1 ) THEN 00373 CALL DCOPY( NRHS, BX, LDBX, B, LDB ) 00374 IF( Z( 1 ).LT.ZERO ) THEN 00375 CALL DSCAL( NRHS, NEGONE, B, LDB ) 00376 END IF 00377 ELSE 00378 DO 50 J = 1, K 00379 DIFLJ = DIFL( J ) 00380 DJ = POLES( J, 1 ) 00381 DSIGJ = -POLES( J, 2 ) 00382 IF( J.LT.K ) THEN 00383 DIFRJ = -DIFR( J, 1 ) 00384 DSIGJP = -POLES( J+1, 2 ) 00385 END IF 00386 IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) ) 00387 $ THEN 00388 WORK( J ) = ZERO 00389 ELSE 00390 WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ / 00391 $ ( POLES( J, 2 )+DJ ) 00392 END IF 00393 DO 30 I = 1, J - 1 00394 IF( ( Z( I ).EQ.ZERO ) .OR. 00395 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN 00396 WORK( I ) = ZERO 00397 ELSE 00398 WORK( I ) = POLES( I, 2 )*Z( I ) / 00399 $ ( DLAMC3( POLES( I, 2 ), DSIGJ )- 00400 $ DIFLJ ) / ( POLES( I, 2 )+DJ ) 00401 END IF 00402 30 CONTINUE 00403 DO 40 I = J + 1, K 00404 IF( ( Z( I ).EQ.ZERO ) .OR. 00405 $ ( POLES( I, 2 ).EQ.ZERO ) ) THEN 00406 WORK( I ) = ZERO 00407 ELSE 00408 WORK( I ) = POLES( I, 2 )*Z( I ) / 00409 $ ( DLAMC3( POLES( I, 2 ), DSIGJP )+ 00410 $ DIFRJ ) / ( POLES( I, 2 )+DJ ) 00411 END IF 00412 40 CONTINUE 00413 WORK( 1 ) = NEGONE 00414 TEMP = DNRM2( K, WORK, 1 ) 00415 CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO, 00416 $ B( J, 1 ), LDB ) 00417 CALL DLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ), 00418 $ LDB, INFO ) 00419 50 CONTINUE 00420 END IF 00421 * 00422 * Move the deflated rows of BX to B also. 00423 * 00424 IF( K.LT.MAX( M, N ) ) 00425 $ CALL DLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX, 00426 $ B( K+1, 1 ), LDB ) 00427 ELSE 00428 * 00429 * Apply back the right orthogonal transformations. 00430 * 00431 * Step (1R): apply back the new right singular vector matrix 00432 * to B. 00433 * 00434 IF( K.EQ.1 ) THEN 00435 CALL DCOPY( NRHS, B, LDB, BX, LDBX ) 00436 ELSE 00437 DO 80 J = 1, K 00438 DSIGJ = POLES( J, 2 ) 00439 IF( Z( J ).EQ.ZERO ) THEN 00440 WORK( J ) = ZERO 00441 ELSE 00442 WORK( J ) = -Z( J ) / DIFL( J ) / 00443 $ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 ) 00444 END IF 00445 DO 60 I = 1, J - 1 00446 IF( Z( J ).EQ.ZERO ) THEN 00447 WORK( I ) = ZERO 00448 ELSE 00449 WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1, 00450 $ 2 ) )-DIFR( I, 1 ) ) / 00451 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) 00452 END IF 00453 60 CONTINUE 00454 DO 70 I = J + 1, K 00455 IF( Z( J ).EQ.ZERO ) THEN 00456 WORK( I ) = ZERO 00457 ELSE 00458 WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I, 00459 $ 2 ) )-DIFL( I ) ) / 00460 $ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) 00461 END IF 00462 70 CONTINUE 00463 CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO, 00464 $ BX( J, 1 ), LDBX ) 00465 80 CONTINUE 00466 END IF 00467 * 00468 * Step (2R): if SQRE = 1, apply back the rotation that is 00469 * related to the right null space of the subproblem. 00470 * 00471 IF( SQRE.EQ.1 ) THEN 00472 CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX ) 00473 CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S ) 00474 END IF 00475 IF( K.LT.MAX( M, N ) ) 00476 $ CALL DLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ), 00477 $ LDBX ) 00478 * 00479 * Step (3R): permute rows of B. 00480 * 00481 CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB ) 00482 IF( SQRE.EQ.1 ) THEN 00483 CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB ) 00484 END IF 00485 DO 90 I = 2, N 00486 CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB ) 00487 90 CONTINUE 00488 * 00489 * Step (4R): apply back the Givens rotations performed. 00490 * 00491 DO 100 I = GIVPTR, 1, -1 00492 CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, 00493 $ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), 00494 $ -GIVNUM( I, 1 ) ) 00495 100 CONTINUE 00496 END IF 00497 * 00498 RETURN 00499 * 00500 * End of DLALS0 00501 * 00502 END