root / src / lapack / double / dlals0.f @ 10
Historique | Voir | Annoter | Télécharger (12,78 ko)
| 1 |
SUBROUTINE DLALS0( ICOMPQ, NL, NR, SQRE, NRHS, B, LDB, BX, LDBX, |
|---|---|
| 2 |
$ PERM, GIVPTR, GIVCOL, LDGCOL, GIVNUM, LDGNUM, |
| 3 |
$ POLES, DIFL, DIFR, Z, K, C, S, WORK, INFO ) |
| 4 |
* |
| 5 |
* -- LAPACK routine (version 3.2) -- |
| 6 |
* -- LAPACK is a software package provided by Univ. of Tennessee, -- |
| 7 |
* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- |
| 8 |
* November 2006 |
| 9 |
* |
| 10 |
* .. Scalar Arguments .. |
| 11 |
INTEGER GIVPTR, ICOMPQ, INFO, K, LDB, LDBX, LDGCOL, |
| 12 |
$ LDGNUM, NL, NR, NRHS, SQRE |
| 13 |
DOUBLE PRECISION C, S |
| 14 |
* .. |
| 15 |
* .. Array Arguments .. |
| 16 |
INTEGER GIVCOL( LDGCOL, * ), PERM( * ) |
| 17 |
DOUBLE PRECISION B( LDB, * ), BX( LDBX, * ), DIFL( * ), |
| 18 |
$ DIFR( LDGNUM, * ), GIVNUM( LDGNUM, * ), |
| 19 |
$ POLES( LDGNUM, * ), WORK( * ), Z( * ) |
| 20 |
* .. |
| 21 |
* |
| 22 |
* Purpose |
| 23 |
* ======= |
| 24 |
* |
| 25 |
* DLALS0 applies back the multiplying factors of either the left or the |
| 26 |
* right singular vector matrix of a diagonal matrix appended by a row |
| 27 |
* to the right hand side matrix B in solving the least squares problem |
| 28 |
* using the divide-and-conquer SVD approach. |
| 29 |
* |
| 30 |
* For the left singular vector matrix, three types of orthogonal |
| 31 |
* matrices are involved: |
| 32 |
* |
| 33 |
* (1L) Givens rotations: the number of such rotations is GIVPTR; the |
| 34 |
* pairs of columns/rows they were applied to are stored in GIVCOL; |
| 35 |
* and the C- and S-values of these rotations are stored in GIVNUM. |
| 36 |
* |
| 37 |
* (2L) Permutation. The (NL+1)-st row of B is to be moved to the first |
| 38 |
* row, and for J=2:N, PERM(J)-th row of B is to be moved to the |
| 39 |
* J-th row. |
| 40 |
* |
| 41 |
* (3L) The left singular vector matrix of the remaining matrix. |
| 42 |
* |
| 43 |
* For the right singular vector matrix, four types of orthogonal |
| 44 |
* matrices are involved: |
| 45 |
* |
| 46 |
* (1R) The right singular vector matrix of the remaining matrix. |
| 47 |
* |
| 48 |
* (2R) If SQRE = 1, one extra Givens rotation to generate the right |
| 49 |
* null space. |
| 50 |
* |
| 51 |
* (3R) The inverse transformation of (2L). |
| 52 |
* |
| 53 |
* (4R) The inverse transformation of (1L). |
| 54 |
* |
| 55 |
* Arguments |
| 56 |
* ========= |
| 57 |
* |
| 58 |
* ICOMPQ (input) INTEGER |
| 59 |
* Specifies whether singular vectors are to be computed in |
| 60 |
* factored form: |
| 61 |
* = 0: Left singular vector matrix. |
| 62 |
* = 1: Right singular vector matrix. |
| 63 |
* |
| 64 |
* NL (input) INTEGER |
| 65 |
* The row dimension of the upper block. NL >= 1. |
| 66 |
* |
| 67 |
* NR (input) INTEGER |
| 68 |
* The row dimension of the lower block. NR >= 1. |
| 69 |
* |
| 70 |
* SQRE (input) INTEGER |
| 71 |
* = 0: the lower block is an NR-by-NR square matrix. |
| 72 |
* = 1: the lower block is an NR-by-(NR+1) rectangular matrix. |
| 73 |
* |
| 74 |
* The bidiagonal matrix has row dimension N = NL + NR + 1, |
| 75 |
* and column dimension M = N + SQRE. |
| 76 |
* |
| 77 |
* NRHS (input) INTEGER |
| 78 |
* The number of columns of B and BX. NRHS must be at least 1. |
| 79 |
* |
| 80 |
* B (input/output) DOUBLE PRECISION array, dimension ( LDB, NRHS ) |
| 81 |
* On input, B contains the right hand sides of the least |
| 82 |
* squares problem in rows 1 through M. On output, B contains |
| 83 |
* the solution X in rows 1 through N. |
| 84 |
* |
| 85 |
* LDB (input) INTEGER |
| 86 |
* The leading dimension of B. LDB must be at least |
| 87 |
* max(1,MAX( M, N ) ). |
| 88 |
* |
| 89 |
* BX (workspace) DOUBLE PRECISION array, dimension ( LDBX, NRHS ) |
| 90 |
* |
| 91 |
* LDBX (input) INTEGER |
| 92 |
* The leading dimension of BX. |
| 93 |
* |
| 94 |
* PERM (input) INTEGER array, dimension ( N ) |
| 95 |
* The permutations (from deflation and sorting) applied |
| 96 |
* to the two blocks. |
| 97 |
* |
| 98 |
* GIVPTR (input) INTEGER |
| 99 |
* The number of Givens rotations which took place in this |
| 100 |
* subproblem. |
| 101 |
* |
| 102 |
* GIVCOL (input) INTEGER array, dimension ( LDGCOL, 2 ) |
| 103 |
* Each pair of numbers indicates a pair of rows/columns |
| 104 |
* involved in a Givens rotation. |
| 105 |
* |
| 106 |
* LDGCOL (input) INTEGER |
| 107 |
* The leading dimension of GIVCOL, must be at least N. |
| 108 |
* |
| 109 |
* GIVNUM (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) |
| 110 |
* Each number indicates the C or S value used in the |
| 111 |
* corresponding Givens rotation. |
| 112 |
* |
| 113 |
* LDGNUM (input) INTEGER |
| 114 |
* The leading dimension of arrays DIFR, POLES and |
| 115 |
* GIVNUM, must be at least K. |
| 116 |
* |
| 117 |
* POLES (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ) |
| 118 |
* On entry, POLES(1:K, 1) contains the new singular |
| 119 |
* values obtained from solving the secular equation, and |
| 120 |
* POLES(1:K, 2) is an array containing the poles in the secular |
| 121 |
* equation. |
| 122 |
* |
| 123 |
* DIFL (input) DOUBLE PRECISION array, dimension ( K ). |
| 124 |
* On entry, DIFL(I) is the distance between I-th updated |
| 125 |
* (undeflated) singular value and the I-th (undeflated) old |
| 126 |
* singular value. |
| 127 |
* |
| 128 |
* DIFR (input) DOUBLE PRECISION array, dimension ( LDGNUM, 2 ). |
| 129 |
* On entry, DIFR(I, 1) contains the distances between I-th |
| 130 |
* updated (undeflated) singular value and the I+1-th |
| 131 |
* (undeflated) old singular value. And DIFR(I, 2) is the |
| 132 |
* normalizing factor for the I-th right singular vector. |
| 133 |
* |
| 134 |
* Z (input) DOUBLE PRECISION array, dimension ( K ) |
| 135 |
* Contain the components of the deflation-adjusted updating row |
| 136 |
* vector. |
| 137 |
* |
| 138 |
* K (input) INTEGER |
| 139 |
* Contains the dimension of the non-deflated matrix, |
| 140 |
* This is the order of the related secular equation. 1 <= K <=N. |
| 141 |
* |
| 142 |
* C (input) DOUBLE PRECISION |
| 143 |
* C contains garbage if SQRE =0 and the C-value of a Givens |
| 144 |
* rotation related to the right null space if SQRE = 1. |
| 145 |
* |
| 146 |
* S (input) DOUBLE PRECISION |
| 147 |
* S contains garbage if SQRE =0 and the S-value of a Givens |
| 148 |
* rotation related to the right null space if SQRE = 1. |
| 149 |
* |
| 150 |
* WORK (workspace) DOUBLE PRECISION array, dimension ( K ) |
| 151 |
* |
| 152 |
* INFO (output) INTEGER |
| 153 |
* = 0: successful exit. |
| 154 |
* < 0: if INFO = -i, the i-th argument had an illegal value. |
| 155 |
* |
| 156 |
* Further Details |
| 157 |
* =============== |
| 158 |
* |
| 159 |
* Based on contributions by |
| 160 |
* Ming Gu and Ren-Cang Li, Computer Science Division, University of |
| 161 |
* California at Berkeley, USA |
| 162 |
* Osni Marques, LBNL/NERSC, USA |
| 163 |
* |
| 164 |
* ===================================================================== |
| 165 |
* |
| 166 |
* .. Parameters .. |
| 167 |
DOUBLE PRECISION ONE, ZERO, NEGONE |
| 168 |
PARAMETER ( ONE = 1.0D0, ZERO = 0.0D0, NEGONE = -1.0D0 ) |
| 169 |
* .. |
| 170 |
* .. Local Scalars .. |
| 171 |
INTEGER I, J, M, N, NLP1 |
| 172 |
DOUBLE PRECISION DIFLJ, DIFRJ, DJ, DSIGJ, DSIGJP, TEMP |
| 173 |
* .. |
| 174 |
* .. External Subroutines .. |
| 175 |
EXTERNAL DCOPY, DGEMV, DLACPY, DLASCL, DROT, DSCAL, |
| 176 |
$ XERBLA |
| 177 |
* .. |
| 178 |
* .. External Functions .. |
| 179 |
DOUBLE PRECISION DLAMC3, DNRM2 |
| 180 |
EXTERNAL DLAMC3, DNRM2 |
| 181 |
* .. |
| 182 |
* .. Intrinsic Functions .. |
| 183 |
INTRINSIC MAX |
| 184 |
* .. |
| 185 |
* .. Executable Statements .. |
| 186 |
* |
| 187 |
* Test the input parameters. |
| 188 |
* |
| 189 |
INFO = 0 |
| 190 |
* |
| 191 |
IF( ( ICOMPQ.LT.0 ) .OR. ( ICOMPQ.GT.1 ) ) THEN |
| 192 |
INFO = -1 |
| 193 |
ELSE IF( NL.LT.1 ) THEN |
| 194 |
INFO = -2 |
| 195 |
ELSE IF( NR.LT.1 ) THEN |
| 196 |
INFO = -3 |
| 197 |
ELSE IF( ( SQRE.LT.0 ) .OR. ( SQRE.GT.1 ) ) THEN |
| 198 |
INFO = -4 |
| 199 |
END IF |
| 200 |
* |
| 201 |
N = NL + NR + 1 |
| 202 |
* |
| 203 |
IF( NRHS.LT.1 ) THEN |
| 204 |
INFO = -5 |
| 205 |
ELSE IF( LDB.LT.N ) THEN |
| 206 |
INFO = -7 |
| 207 |
ELSE IF( LDBX.LT.N ) THEN |
| 208 |
INFO = -9 |
| 209 |
ELSE IF( GIVPTR.LT.0 ) THEN |
| 210 |
INFO = -11 |
| 211 |
ELSE IF( LDGCOL.LT.N ) THEN |
| 212 |
INFO = -13 |
| 213 |
ELSE IF( LDGNUM.LT.N ) THEN |
| 214 |
INFO = -15 |
| 215 |
ELSE IF( K.LT.1 ) THEN |
| 216 |
INFO = -20 |
| 217 |
END IF |
| 218 |
IF( INFO.NE.0 ) THEN |
| 219 |
CALL XERBLA( 'DLALS0', -INFO ) |
| 220 |
RETURN |
| 221 |
END IF |
| 222 |
* |
| 223 |
M = N + SQRE |
| 224 |
NLP1 = NL + 1 |
| 225 |
* |
| 226 |
IF( ICOMPQ.EQ.0 ) THEN |
| 227 |
* |
| 228 |
* Apply back orthogonal transformations from the left. |
| 229 |
* |
| 230 |
* Step (1L): apply back the Givens rotations performed. |
| 231 |
* |
| 232 |
DO 10 I = 1, GIVPTR |
| 233 |
CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, |
| 234 |
$ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), |
| 235 |
$ GIVNUM( I, 1 ) ) |
| 236 |
10 CONTINUE |
| 237 |
* |
| 238 |
* Step (2L): permute rows of B. |
| 239 |
* |
| 240 |
CALL DCOPY( NRHS, B( NLP1, 1 ), LDB, BX( 1, 1 ), LDBX ) |
| 241 |
DO 20 I = 2, N |
| 242 |
CALL DCOPY( NRHS, B( PERM( I ), 1 ), LDB, BX( I, 1 ), LDBX ) |
| 243 |
20 CONTINUE |
| 244 |
* |
| 245 |
* Step (3L): apply the inverse of the left singular vector |
| 246 |
* matrix to BX. |
| 247 |
* |
| 248 |
IF( K.EQ.1 ) THEN |
| 249 |
CALL DCOPY( NRHS, BX, LDBX, B, LDB ) |
| 250 |
IF( Z( 1 ).LT.ZERO ) THEN |
| 251 |
CALL DSCAL( NRHS, NEGONE, B, LDB ) |
| 252 |
END IF |
| 253 |
ELSE |
| 254 |
DO 50 J = 1, K |
| 255 |
DIFLJ = DIFL( J ) |
| 256 |
DJ = POLES( J, 1 ) |
| 257 |
DSIGJ = -POLES( J, 2 ) |
| 258 |
IF( J.LT.K ) THEN |
| 259 |
DIFRJ = -DIFR( J, 1 ) |
| 260 |
DSIGJP = -POLES( J+1, 2 ) |
| 261 |
END IF |
| 262 |
IF( ( Z( J ).EQ.ZERO ) .OR. ( POLES( J, 2 ).EQ.ZERO ) ) |
| 263 |
$ THEN |
| 264 |
WORK( J ) = ZERO |
| 265 |
ELSE |
| 266 |
WORK( J ) = -POLES( J, 2 )*Z( J ) / DIFLJ / |
| 267 |
$ ( POLES( J, 2 )+DJ ) |
| 268 |
END IF |
| 269 |
DO 30 I = 1, J - 1 |
| 270 |
IF( ( Z( I ).EQ.ZERO ) .OR. |
| 271 |
$ ( POLES( I, 2 ).EQ.ZERO ) ) THEN |
| 272 |
WORK( I ) = ZERO |
| 273 |
ELSE |
| 274 |
WORK( I ) = POLES( I, 2 )*Z( I ) / |
| 275 |
$ ( DLAMC3( POLES( I, 2 ), DSIGJ )- |
| 276 |
$ DIFLJ ) / ( POLES( I, 2 )+DJ ) |
| 277 |
END IF |
| 278 |
30 CONTINUE |
| 279 |
DO 40 I = J + 1, K |
| 280 |
IF( ( Z( I ).EQ.ZERO ) .OR. |
| 281 |
$ ( POLES( I, 2 ).EQ.ZERO ) ) THEN |
| 282 |
WORK( I ) = ZERO |
| 283 |
ELSE |
| 284 |
WORK( I ) = POLES( I, 2 )*Z( I ) / |
| 285 |
$ ( DLAMC3( POLES( I, 2 ), DSIGJP )+ |
| 286 |
$ DIFRJ ) / ( POLES( I, 2 )+DJ ) |
| 287 |
END IF |
| 288 |
40 CONTINUE |
| 289 |
WORK( 1 ) = NEGONE |
| 290 |
TEMP = DNRM2( K, WORK, 1 ) |
| 291 |
CALL DGEMV( 'T', K, NRHS, ONE, BX, LDBX, WORK, 1, ZERO, |
| 292 |
$ B( J, 1 ), LDB ) |
| 293 |
CALL DLASCL( 'G', 0, 0, TEMP, ONE, 1, NRHS, B( J, 1 ), |
| 294 |
$ LDB, INFO ) |
| 295 |
50 CONTINUE |
| 296 |
END IF |
| 297 |
* |
| 298 |
* Move the deflated rows of BX to B also. |
| 299 |
* |
| 300 |
IF( K.LT.MAX( M, N ) ) |
| 301 |
$ CALL DLACPY( 'A', N-K, NRHS, BX( K+1, 1 ), LDBX, |
| 302 |
$ B( K+1, 1 ), LDB ) |
| 303 |
ELSE |
| 304 |
* |
| 305 |
* Apply back the right orthogonal transformations. |
| 306 |
* |
| 307 |
* Step (1R): apply back the new right singular vector matrix |
| 308 |
* to B. |
| 309 |
* |
| 310 |
IF( K.EQ.1 ) THEN |
| 311 |
CALL DCOPY( NRHS, B, LDB, BX, LDBX ) |
| 312 |
ELSE |
| 313 |
DO 80 J = 1, K |
| 314 |
DSIGJ = POLES( J, 2 ) |
| 315 |
IF( Z( J ).EQ.ZERO ) THEN |
| 316 |
WORK( J ) = ZERO |
| 317 |
ELSE |
| 318 |
WORK( J ) = -Z( J ) / DIFL( J ) / |
| 319 |
$ ( DSIGJ+POLES( J, 1 ) ) / DIFR( J, 2 ) |
| 320 |
END IF |
| 321 |
DO 60 I = 1, J - 1 |
| 322 |
IF( Z( J ).EQ.ZERO ) THEN |
| 323 |
WORK( I ) = ZERO |
| 324 |
ELSE |
| 325 |
WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I+1, |
| 326 |
$ 2 ) )-DIFR( I, 1 ) ) / |
| 327 |
$ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) |
| 328 |
END IF |
| 329 |
60 CONTINUE |
| 330 |
DO 70 I = J + 1, K |
| 331 |
IF( Z( J ).EQ.ZERO ) THEN |
| 332 |
WORK( I ) = ZERO |
| 333 |
ELSE |
| 334 |
WORK( I ) = Z( J ) / ( DLAMC3( DSIGJ, -POLES( I, |
| 335 |
$ 2 ) )-DIFL( I ) ) / |
| 336 |
$ ( DSIGJ+POLES( I, 1 ) ) / DIFR( I, 2 ) |
| 337 |
END IF |
| 338 |
70 CONTINUE |
| 339 |
CALL DGEMV( 'T', K, NRHS, ONE, B, LDB, WORK, 1, ZERO, |
| 340 |
$ BX( J, 1 ), LDBX ) |
| 341 |
80 CONTINUE |
| 342 |
END IF |
| 343 |
* |
| 344 |
* Step (2R): if SQRE = 1, apply back the rotation that is |
| 345 |
* related to the right null space of the subproblem. |
| 346 |
* |
| 347 |
IF( SQRE.EQ.1 ) THEN |
| 348 |
CALL DCOPY( NRHS, B( M, 1 ), LDB, BX( M, 1 ), LDBX ) |
| 349 |
CALL DROT( NRHS, BX( 1, 1 ), LDBX, BX( M, 1 ), LDBX, C, S ) |
| 350 |
END IF |
| 351 |
IF( K.LT.MAX( M, N ) ) |
| 352 |
$ CALL DLACPY( 'A', N-K, NRHS, B( K+1, 1 ), LDB, BX( K+1, 1 ), |
| 353 |
$ LDBX ) |
| 354 |
* |
| 355 |
* Step (3R): permute rows of B. |
| 356 |
* |
| 357 |
CALL DCOPY( NRHS, BX( 1, 1 ), LDBX, B( NLP1, 1 ), LDB ) |
| 358 |
IF( SQRE.EQ.1 ) THEN |
| 359 |
CALL DCOPY( NRHS, BX( M, 1 ), LDBX, B( M, 1 ), LDB ) |
| 360 |
END IF |
| 361 |
DO 90 I = 2, N |
| 362 |
CALL DCOPY( NRHS, BX( I, 1 ), LDBX, B( PERM( I ), 1 ), LDB ) |
| 363 |
90 CONTINUE |
| 364 |
* |
| 365 |
* Step (4R): apply back the Givens rotations performed. |
| 366 |
* |
| 367 |
DO 100 I = GIVPTR, 1, -1 |
| 368 |
CALL DROT( NRHS, B( GIVCOL( I, 2 ), 1 ), LDB, |
| 369 |
$ B( GIVCOL( I, 1 ), 1 ), LDB, GIVNUM( I, 2 ), |
| 370 |
$ -GIVNUM( I, 1 ) ) |
| 371 |
100 CONTINUE |
| 372 |
END IF |
| 373 |
* |
| 374 |
RETURN |
| 375 |
* |
| 376 |
* End of DLALS0 |
| 377 |
* |
| 378 |
END |