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
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