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