1 SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
2 $ RANK, WORK, RWORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER UPLO
11 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
12 DOUBLE PRECISION RCOND
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * )
16 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
17 COMPLEX*16 B( LDB, * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZLALSD uses the singular value decomposition of A to solve the least
24 * squares problem of finding X to minimize the Euclidean norm of each
25 * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
26 * are N-by-NRHS. The solution X overwrites B.
27 *
28 * The singular values of A smaller than RCOND times the largest
29 * singular value are treated as zero in solving the least squares
30 * problem; in this case a minimum norm solution is returned.
31 * The actual singular values are returned in D in ascending order.
32 *
33 * This code makes very mild assumptions about floating point
34 * arithmetic. It will work on machines with a guard digit in
35 * add/subtract, or on those binary machines without guard digits
36 * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
37 * It could conceivably fail on hexadecimal or decimal machines
38 * without guard digits, but we know of none.
39 *
40 * Arguments
41 * =========
42 *
43 * UPLO (input) CHARACTER*1
44 * = 'U': D and E define an upper bidiagonal matrix.
45 * = 'L': D and E define a lower bidiagonal matrix.
46 *
47 * SMLSIZ (input) INTEGER
48 * The maximum size of the subproblems at the bottom of the
49 * computation tree.
50 *
51 * N (input) INTEGER
52 * The dimension of the bidiagonal matrix. N >= 0.
53 *
54 * NRHS (input) INTEGER
55 * The number of columns of B. NRHS must be at least 1.
56 *
57 * D (input/output) DOUBLE PRECISION array, dimension (N)
58 * On entry D contains the main diagonal of the bidiagonal
59 * matrix. On exit, if INFO = 0, D contains its singular values.
60 *
61 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
62 * Contains the super-diagonal entries of the bidiagonal matrix.
63 * On exit, E has been destroyed.
64 *
65 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
66 * On input, B contains the right hand sides of the least
67 * squares problem. On output, B contains the solution X.
68 *
69 * LDB (input) INTEGER
70 * The leading dimension of B in the calling subprogram.
71 * LDB must be at least max(1,N).
72 *
73 * RCOND (input) DOUBLE PRECISION
74 * The singular values of A less than or equal to RCOND times
75 * the largest singular value are treated as zero in solving
76 * the least squares problem. If RCOND is negative,
77 * machine precision is used instead.
78 * For example, if diag(S)*X=B were the least squares problem,
79 * where diag(S) is a diagonal matrix of singular values, the
80 * solution would be X(i) = B(i) / S(i) if S(i) is greater than
81 * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
82 * RCOND*max(S).
83 *
84 * RANK (output) INTEGER
85 * The number of singular values of A greater than RCOND times
86 * the largest singular value.
87 *
88 * WORK (workspace) COMPLEX*16 array, dimension at least
89 * (N * NRHS).
90 *
91 * RWORK (workspace) DOUBLE PRECISION array, dimension at least
92 * (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
93 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
94 * where
95 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
96 *
97 * IWORK (workspace) INTEGER array, dimension at least
98 * (3*N*NLVL + 11*N).
99 *
100 * INFO (output) INTEGER
101 * = 0: successful exit.
102 * < 0: if INFO = -i, the i-th argument had an illegal value.
103 * > 0: The algorithm failed to compute a singular value while
104 * working on the submatrix lying in rows and columns
105 * INFO/(N+1) through MOD(INFO,N+1).
106 *
107 * Further Details
108 * ===============
109 *
110 * Based on contributions by
111 * Ming Gu and Ren-Cang Li, Computer Science Division, University of
112 * California at Berkeley, USA
113 * Osni Marques, LBNL/NERSC, USA
114 *
115 * =====================================================================
116 *
117 * .. Parameters ..
118 DOUBLE PRECISION ZERO, ONE, TWO
119 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
120 COMPLEX*16 CZERO
121 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) )
122 * ..
123 * .. Local Scalars ..
124 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
125 $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
126 $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
127 $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
128 $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
129 $ U, VT, Z
130 DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
131 * ..
132 * .. External Functions ..
133 INTEGER IDAMAX
134 DOUBLE PRECISION DLAMCH, DLANST
135 EXTERNAL IDAMAX, DLAMCH, DLANST
136 * ..
137 * .. External Subroutines ..
138 EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
139 $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
140 $ ZLASCL, ZLASET
141 * ..
142 * .. Intrinsic Functions ..
143 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
144 * ..
145 * .. Executable Statements ..
146 *
147 * Test the input parameters.
148 *
149 INFO = 0
150 *
151 IF( N.LT.0 ) THEN
152 INFO = -3
153 ELSE IF( NRHS.LT.1 ) THEN
154 INFO = -4
155 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
156 INFO = -8
157 END IF
158 IF( INFO.NE.0 ) THEN
159 CALL XERBLA( 'ZLALSD', -INFO )
160 RETURN
161 END IF
162 *
163 EPS = DLAMCH( 'Epsilon' )
164 *
165 * Set up the tolerance.
166 *
167 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
168 RCND = EPS
169 ELSE
170 RCND = RCOND
171 END IF
172 *
173 RANK = 0
174 *
175 * Quick return if possible.
176 *
177 IF( N.EQ.0 ) THEN
178 RETURN
179 ELSE IF( N.EQ.1 ) THEN
180 IF( D( 1 ).EQ.ZERO ) THEN
181 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
182 ELSE
183 RANK = 1
184 CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
185 D( 1 ) = ABS( D( 1 ) )
186 END IF
187 RETURN
188 END IF
189 *
190 * Rotate the matrix if it is lower bidiagonal.
191 *
192 IF( UPLO.EQ.'L' ) THEN
193 DO 10 I = 1, N - 1
194 CALL DLARTG( D( I ), E( I ), CS, SN, R )
195 D( I ) = R
196 E( I ) = SN*D( I+1 )
197 D( I+1 ) = CS*D( I+1 )
198 IF( NRHS.EQ.1 ) THEN
199 CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
200 ELSE
201 RWORK( I*2-1 ) = CS
202 RWORK( I*2 ) = SN
203 END IF
204 10 CONTINUE
205 IF( NRHS.GT.1 ) THEN
206 DO 30 I = 1, NRHS
207 DO 20 J = 1, N - 1
208 CS = RWORK( J*2-1 )
209 SN = RWORK( J*2 )
210 CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
211 20 CONTINUE
212 30 CONTINUE
213 END IF
214 END IF
215 *
216 * Scale.
217 *
218 NM1 = N - 1
219 ORGNRM = DLANST( 'M', N, D, E )
220 IF( ORGNRM.EQ.ZERO ) THEN
221 CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
222 RETURN
223 END IF
224 *
225 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
226 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
227 *
228 * If N is smaller than the minimum divide size SMLSIZ, then solve
229 * the problem with another solver.
230 *
231 IF( N.LE.SMLSIZ ) THEN
232 IRWU = 1
233 IRWVT = IRWU + N*N
234 IRWWRK = IRWVT + N*N
235 IRWRB = IRWWRK
236 IRWIB = IRWRB + N*NRHS
237 IRWB = IRWIB + N*NRHS
238 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
239 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
240 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
241 $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
242 $ RWORK( IRWWRK ), INFO )
243 IF( INFO.NE.0 ) THEN
244 RETURN
245 END IF
246 *
247 * In the real version, B is passed to DLASDQ and multiplied
248 * internally by Q**H. Here B is complex and that product is
249 * computed below in two steps (real and imaginary parts).
250 *
251 J = IRWB - 1
252 DO 50 JCOL = 1, NRHS
253 DO 40 JROW = 1, N
254 J = J + 1
255 RWORK( J ) = DBLE( B( JROW, JCOL ) )
256 40 CONTINUE
257 50 CONTINUE
258 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
259 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
260 J = IRWB - 1
261 DO 70 JCOL = 1, NRHS
262 DO 60 JROW = 1, N
263 J = J + 1
264 RWORK( J ) = DIMAG( B( JROW, JCOL ) )
265 60 CONTINUE
266 70 CONTINUE
267 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
268 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
269 JREAL = IRWRB - 1
270 JIMAG = IRWIB - 1
271 DO 90 JCOL = 1, NRHS
272 DO 80 JROW = 1, N
273 JREAL = JREAL + 1
274 JIMAG = JIMAG + 1
275 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
276 $ RWORK( JIMAG ) )
277 80 CONTINUE
278 90 CONTINUE
279 *
280 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
281 DO 100 I = 1, N
282 IF( D( I ).LE.TOL ) THEN
283 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
284 ELSE
285 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
286 $ LDB, INFO )
287 RANK = RANK + 1
288 END IF
289 100 CONTINUE
290 *
291 * Since B is complex, the following call to DGEMM is performed
292 * in two steps (real and imaginary parts). That is for V * B
293 * (in the real version of the code V**H is stored in WORK).
294 *
295 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
296 * $ WORK( NWORK ), N )
297 *
298 J = IRWB - 1
299 DO 120 JCOL = 1, NRHS
300 DO 110 JROW = 1, N
301 J = J + 1
302 RWORK( J ) = DBLE( B( JROW, JCOL ) )
303 110 CONTINUE
304 120 CONTINUE
305 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
306 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
307 J = IRWB - 1
308 DO 140 JCOL = 1, NRHS
309 DO 130 JROW = 1, N
310 J = J + 1
311 RWORK( J ) = DIMAG( B( JROW, JCOL ) )
312 130 CONTINUE
313 140 CONTINUE
314 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
315 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
316 JREAL = IRWRB - 1
317 JIMAG = IRWIB - 1
318 DO 160 JCOL = 1, NRHS
319 DO 150 JROW = 1, N
320 JREAL = JREAL + 1
321 JIMAG = JIMAG + 1
322 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
323 $ RWORK( JIMAG ) )
324 150 CONTINUE
325 160 CONTINUE
326 *
327 * Unscale.
328 *
329 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
330 CALL DLASRT( 'D', N, D, INFO )
331 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
332 *
333 RETURN
334 END IF
335 *
336 * Book-keeping and setting up some constants.
337 *
338 NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
339 *
340 SMLSZP = SMLSIZ + 1
341 *
342 U = 1
343 VT = 1 + SMLSIZ*N
344 DIFL = VT + SMLSZP*N
345 DIFR = DIFL + NLVL*N
346 Z = DIFR + NLVL*N*2
347 C = Z + NLVL*N
348 S = C + N
349 POLES = S + N
350 GIVNUM = POLES + 2*NLVL*N
351 NRWORK = GIVNUM + 2*NLVL*N
352 BX = 1
353 *
354 IRWRB = NRWORK
355 IRWIB = IRWRB + SMLSIZ*NRHS
356 IRWB = IRWIB + SMLSIZ*NRHS
357 *
358 SIZEI = 1 + N
359 K = SIZEI + N
360 GIVPTR = K + N
361 PERM = GIVPTR + N
362 GIVCOL = PERM + NLVL*N
363 IWK = GIVCOL + NLVL*N*2
364 *
365 ST = 1
366 SQRE = 0
367 ICMPQ1 = 1
368 ICMPQ2 = 0
369 NSUB = 0
370 *
371 DO 170 I = 1, N
372 IF( ABS( D( I ) ).LT.EPS ) THEN
373 D( I ) = SIGN( EPS, D( I ) )
374 END IF
375 170 CONTINUE
376 *
377 DO 240 I = 1, NM1
378 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
379 NSUB = NSUB + 1
380 IWORK( NSUB ) = ST
381 *
382 * Subproblem found. First determine its size and then
383 * apply divide and conquer on it.
384 *
385 IF( I.LT.NM1 ) THEN
386 *
387 * A subproblem with E(I) small for I < NM1.
388 *
389 NSIZE = I - ST + 1
390 IWORK( SIZEI+NSUB-1 ) = NSIZE
391 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
392 *
393 * A subproblem with E(NM1) not too small but I = NM1.
394 *
395 NSIZE = N - ST + 1
396 IWORK( SIZEI+NSUB-1 ) = NSIZE
397 ELSE
398 *
399 * A subproblem with E(NM1) small. This implies an
400 * 1-by-1 subproblem at D(N), which is not solved
401 * explicitly.
402 *
403 NSIZE = I - ST + 1
404 IWORK( SIZEI+NSUB-1 ) = NSIZE
405 NSUB = NSUB + 1
406 IWORK( NSUB ) = N
407 IWORK( SIZEI+NSUB-1 ) = 1
408 CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
409 END IF
410 ST1 = ST - 1
411 IF( NSIZE.EQ.1 ) THEN
412 *
413 * This is a 1-by-1 subproblem and is not solved
414 * explicitly.
415 *
416 CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
417 ELSE IF( NSIZE.LE.SMLSIZ ) THEN
418 *
419 * This is a small subproblem and is solved by DLASDQ.
420 *
421 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
422 $ RWORK( VT+ST1 ), N )
423 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
424 $ RWORK( U+ST1 ), N )
425 CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
426 $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
427 $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
428 $ INFO )
429 IF( INFO.NE.0 ) THEN
430 RETURN
431 END IF
432 *
433 * In the real version, B is passed to DLASDQ and multiplied
434 * internally by Q**H. Here B is complex and that product is
435 * computed below in two steps (real and imaginary parts).
436 *
437 J = IRWB - 1
438 DO 190 JCOL = 1, NRHS
439 DO 180 JROW = ST, ST + NSIZE - 1
440 J = J + 1
441 RWORK( J ) = DBLE( B( JROW, JCOL ) )
442 180 CONTINUE
443 190 CONTINUE
444 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
445 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
446 $ ZERO, RWORK( IRWRB ), NSIZE )
447 J = IRWB - 1
448 DO 210 JCOL = 1, NRHS
449 DO 200 JROW = ST, ST + NSIZE - 1
450 J = J + 1
451 RWORK( J ) = DIMAG( B( JROW, JCOL ) )
452 200 CONTINUE
453 210 CONTINUE
454 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
455 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
456 $ ZERO, RWORK( IRWIB ), NSIZE )
457 JREAL = IRWRB - 1
458 JIMAG = IRWIB - 1
459 DO 230 JCOL = 1, NRHS
460 DO 220 JROW = ST, ST + NSIZE - 1
461 JREAL = JREAL + 1
462 JIMAG = JIMAG + 1
463 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
464 $ RWORK( JIMAG ) )
465 220 CONTINUE
466 230 CONTINUE
467 *
468 CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
469 $ WORK( BX+ST1 ), N )
470 ELSE
471 *
472 * A large problem. Solve it using divide and conquer.
473 *
474 CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
475 $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
476 $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
477 $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
478 $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
479 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
480 $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
481 $ RWORK( S+ST1 ), RWORK( NRWORK ),
482 $ IWORK( IWK ), INFO )
483 IF( INFO.NE.0 ) THEN
484 RETURN
485 END IF
486 BXST = BX + ST1
487 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
488 $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
489 $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
490 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
491 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
492 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
493 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
494 $ RWORK( C+ST1 ), RWORK( S+ST1 ),
495 $ RWORK( NRWORK ), IWORK( IWK ), INFO )
496 IF( INFO.NE.0 ) THEN
497 RETURN
498 END IF
499 END IF
500 ST = I + 1
501 END IF
502 240 CONTINUE
503 *
504 * Apply the singular values and treat the tiny ones as zero.
505 *
506 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
507 *
508 DO 250 I = 1, N
509 *
510 * Some of the elements in D can be negative because 1-by-1
511 * subproblems were not solved explicitly.
512 *
513 IF( ABS( D( I ) ).LE.TOL ) THEN
514 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
515 ELSE
516 RANK = RANK + 1
517 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
518 $ WORK( BX+I-1 ), N, INFO )
519 END IF
520 D( I ) = ABS( D( I ) )
521 250 CONTINUE
522 *
523 * Now apply back the right singular vectors.
524 *
525 ICMPQ2 = 1
526 DO 320 I = 1, NSUB
527 ST = IWORK( I )
528 ST1 = ST - 1
529 NSIZE = IWORK( SIZEI+I-1 )
530 BXST = BX + ST1
531 IF( NSIZE.EQ.1 ) THEN
532 CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
533 ELSE IF( NSIZE.LE.SMLSIZ ) THEN
534 *
535 * Since B and BX are complex, the following call to DGEMM
536 * is performed in two steps (real and imaginary parts).
537 *
538 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
539 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
540 * $ B( ST, 1 ), LDB )
541 *
542 J = BXST - N - 1
543 JREAL = IRWB - 1
544 DO 270 JCOL = 1, NRHS
545 J = J + N
546 DO 260 JROW = 1, NSIZE
547 JREAL = JREAL + 1
548 RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
549 260 CONTINUE
550 270 CONTINUE
551 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
552 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
553 $ RWORK( IRWRB ), NSIZE )
554 J = BXST - N - 1
555 JIMAG = IRWB - 1
556 DO 290 JCOL = 1, NRHS
557 J = J + N
558 DO 280 JROW = 1, NSIZE
559 JIMAG = JIMAG + 1
560 RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
561 280 CONTINUE
562 290 CONTINUE
563 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
564 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
565 $ RWORK( IRWIB ), NSIZE )
566 JREAL = IRWRB - 1
567 JIMAG = IRWIB - 1
568 DO 310 JCOL = 1, NRHS
569 DO 300 JROW = ST, ST + NSIZE - 1
570 JREAL = JREAL + 1
571 JIMAG = JIMAG + 1
572 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
573 $ RWORK( JIMAG ) )
574 300 CONTINUE
575 310 CONTINUE
576 ELSE
577 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
578 $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
579 $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
580 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
581 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
582 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
583 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
584 $ RWORK( C+ST1 ), RWORK( S+ST1 ),
585 $ RWORK( NRWORK ), IWORK( IWK ), INFO )
586 IF( INFO.NE.0 ) THEN
587 RETURN
588 END IF
589 END IF
590 320 CONTINUE
591 *
592 * Unscale and sort the singular values.
593 *
594 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
595 CALL DLASRT( 'D', N, D, INFO )
596 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
597 *
598 RETURN
599 *
600 * End of ZLALSD
601 *
602 END
2 $ RANK, WORK, RWORK, IWORK, INFO )
3 *
4 * -- LAPACK routine (version 3.3.1) --
5 * -- LAPACK is a software package provided by Univ. of Tennessee, --
6 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7 * -- April 2011 --
8 *
9 * .. Scalar Arguments ..
10 CHARACTER UPLO
11 INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
12 DOUBLE PRECISION RCOND
13 * ..
14 * .. Array Arguments ..
15 INTEGER IWORK( * )
16 DOUBLE PRECISION D( * ), E( * ), RWORK( * )
17 COMPLEX*16 B( LDB, * ), WORK( * )
18 * ..
19 *
20 * Purpose
21 * =======
22 *
23 * ZLALSD uses the singular value decomposition of A to solve the least
24 * squares problem of finding X to minimize the Euclidean norm of each
25 * column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
26 * are N-by-NRHS. The solution X overwrites B.
27 *
28 * The singular values of A smaller than RCOND times the largest
29 * singular value are treated as zero in solving the least squares
30 * problem; in this case a minimum norm solution is returned.
31 * The actual singular values are returned in D in ascending order.
32 *
33 * This code makes very mild assumptions about floating point
34 * arithmetic. It will work on machines with a guard digit in
35 * add/subtract, or on those binary machines without guard digits
36 * which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
37 * It could conceivably fail on hexadecimal or decimal machines
38 * without guard digits, but we know of none.
39 *
40 * Arguments
41 * =========
42 *
43 * UPLO (input) CHARACTER*1
44 * = 'U': D and E define an upper bidiagonal matrix.
45 * = 'L': D and E define a lower bidiagonal matrix.
46 *
47 * SMLSIZ (input) INTEGER
48 * The maximum size of the subproblems at the bottom of the
49 * computation tree.
50 *
51 * N (input) INTEGER
52 * The dimension of the bidiagonal matrix. N >= 0.
53 *
54 * NRHS (input) INTEGER
55 * The number of columns of B. NRHS must be at least 1.
56 *
57 * D (input/output) DOUBLE PRECISION array, dimension (N)
58 * On entry D contains the main diagonal of the bidiagonal
59 * matrix. On exit, if INFO = 0, D contains its singular values.
60 *
61 * E (input/output) DOUBLE PRECISION array, dimension (N-1)
62 * Contains the super-diagonal entries of the bidiagonal matrix.
63 * On exit, E has been destroyed.
64 *
65 * B (input/output) COMPLEX*16 array, dimension (LDB,NRHS)
66 * On input, B contains the right hand sides of the least
67 * squares problem. On output, B contains the solution X.
68 *
69 * LDB (input) INTEGER
70 * The leading dimension of B in the calling subprogram.
71 * LDB must be at least max(1,N).
72 *
73 * RCOND (input) DOUBLE PRECISION
74 * The singular values of A less than or equal to RCOND times
75 * the largest singular value are treated as zero in solving
76 * the least squares problem. If RCOND is negative,
77 * machine precision is used instead.
78 * For example, if diag(S)*X=B were the least squares problem,
79 * where diag(S) is a diagonal matrix of singular values, the
80 * solution would be X(i) = B(i) / S(i) if S(i) is greater than
81 * RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
82 * RCOND*max(S).
83 *
84 * RANK (output) INTEGER
85 * The number of singular values of A greater than RCOND times
86 * the largest singular value.
87 *
88 * WORK (workspace) COMPLEX*16 array, dimension at least
89 * (N * NRHS).
90 *
91 * RWORK (workspace) DOUBLE PRECISION array, dimension at least
92 * (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
93 * MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
94 * where
95 * NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
96 *
97 * IWORK (workspace) INTEGER array, dimension at least
98 * (3*N*NLVL + 11*N).
99 *
100 * INFO (output) INTEGER
101 * = 0: successful exit.
102 * < 0: if INFO = -i, the i-th argument had an illegal value.
103 * > 0: The algorithm failed to compute a singular value while
104 * working on the submatrix lying in rows and columns
105 * INFO/(N+1) through MOD(INFO,N+1).
106 *
107 * Further Details
108 * ===============
109 *
110 * Based on contributions by
111 * Ming Gu and Ren-Cang Li, Computer Science Division, University of
112 * California at Berkeley, USA
113 * Osni Marques, LBNL/NERSC, USA
114 *
115 * =====================================================================
116 *
117 * .. Parameters ..
118 DOUBLE PRECISION ZERO, ONE, TWO
119 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
120 COMPLEX*16 CZERO
121 PARAMETER ( CZERO = ( 0.0D0, 0.0D0 ) )
122 * ..
123 * .. Local Scalars ..
124 INTEGER BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
125 $ GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
126 $ IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
127 $ JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
128 $ PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
129 $ U, VT, Z
130 DOUBLE PRECISION CS, EPS, ORGNRM, RCND, R, SN, TOL
131 * ..
132 * .. External Functions ..
133 INTEGER IDAMAX
134 DOUBLE PRECISION DLAMCH, DLANST
135 EXTERNAL IDAMAX, DLAMCH, DLANST
136 * ..
137 * .. External Subroutines ..
138 EXTERNAL DGEMM, DLARTG, DLASCL, DLASDA, DLASDQ, DLASET,
139 $ DLASRT, XERBLA, ZCOPY, ZDROT, ZLACPY, ZLALSA,
140 $ ZLASCL, ZLASET
141 * ..
142 * .. Intrinsic Functions ..
143 INTRINSIC ABS, DBLE, DCMPLX, DIMAG, INT, LOG, SIGN
144 * ..
145 * .. Executable Statements ..
146 *
147 * Test the input parameters.
148 *
149 INFO = 0
150 *
151 IF( N.LT.0 ) THEN
152 INFO = -3
153 ELSE IF( NRHS.LT.1 ) THEN
154 INFO = -4
155 ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
156 INFO = -8
157 END IF
158 IF( INFO.NE.0 ) THEN
159 CALL XERBLA( 'ZLALSD', -INFO )
160 RETURN
161 END IF
162 *
163 EPS = DLAMCH( 'Epsilon' )
164 *
165 * Set up the tolerance.
166 *
167 IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
168 RCND = EPS
169 ELSE
170 RCND = RCOND
171 END IF
172 *
173 RANK = 0
174 *
175 * Quick return if possible.
176 *
177 IF( N.EQ.0 ) THEN
178 RETURN
179 ELSE IF( N.EQ.1 ) THEN
180 IF( D( 1 ).EQ.ZERO ) THEN
181 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
182 ELSE
183 RANK = 1
184 CALL ZLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
185 D( 1 ) = ABS( D( 1 ) )
186 END IF
187 RETURN
188 END IF
189 *
190 * Rotate the matrix if it is lower bidiagonal.
191 *
192 IF( UPLO.EQ.'L' ) THEN
193 DO 10 I = 1, N - 1
194 CALL DLARTG( D( I ), E( I ), CS, SN, R )
195 D( I ) = R
196 E( I ) = SN*D( I+1 )
197 D( I+1 ) = CS*D( I+1 )
198 IF( NRHS.EQ.1 ) THEN
199 CALL ZDROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
200 ELSE
201 RWORK( I*2-1 ) = CS
202 RWORK( I*2 ) = SN
203 END IF
204 10 CONTINUE
205 IF( NRHS.GT.1 ) THEN
206 DO 30 I = 1, NRHS
207 DO 20 J = 1, N - 1
208 CS = RWORK( J*2-1 )
209 SN = RWORK( J*2 )
210 CALL ZDROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
211 20 CONTINUE
212 30 CONTINUE
213 END IF
214 END IF
215 *
216 * Scale.
217 *
218 NM1 = N - 1
219 ORGNRM = DLANST( 'M', N, D, E )
220 IF( ORGNRM.EQ.ZERO ) THEN
221 CALL ZLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
222 RETURN
223 END IF
224 *
225 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
226 CALL DLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
227 *
228 * If N is smaller than the minimum divide size SMLSIZ, then solve
229 * the problem with another solver.
230 *
231 IF( N.LE.SMLSIZ ) THEN
232 IRWU = 1
233 IRWVT = IRWU + N*N
234 IRWWRK = IRWVT + N*N
235 IRWRB = IRWWRK
236 IRWIB = IRWRB + N*NRHS
237 IRWB = IRWIB + N*NRHS
238 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
239 CALL DLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
240 CALL DLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
241 $ RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
242 $ RWORK( IRWWRK ), INFO )
243 IF( INFO.NE.0 ) THEN
244 RETURN
245 END IF
246 *
247 * In the real version, B is passed to DLASDQ and multiplied
248 * internally by Q**H. Here B is complex and that product is
249 * computed below in two steps (real and imaginary parts).
250 *
251 J = IRWB - 1
252 DO 50 JCOL = 1, NRHS
253 DO 40 JROW = 1, N
254 J = J + 1
255 RWORK( J ) = DBLE( B( JROW, JCOL ) )
256 40 CONTINUE
257 50 CONTINUE
258 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
259 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
260 J = IRWB - 1
261 DO 70 JCOL = 1, NRHS
262 DO 60 JROW = 1, N
263 J = J + 1
264 RWORK( J ) = DIMAG( B( JROW, JCOL ) )
265 60 CONTINUE
266 70 CONTINUE
267 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
268 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
269 JREAL = IRWRB - 1
270 JIMAG = IRWIB - 1
271 DO 90 JCOL = 1, NRHS
272 DO 80 JROW = 1, N
273 JREAL = JREAL + 1
274 JIMAG = JIMAG + 1
275 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
276 $ RWORK( JIMAG ) )
277 80 CONTINUE
278 90 CONTINUE
279 *
280 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
281 DO 100 I = 1, N
282 IF( D( I ).LE.TOL ) THEN
283 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
284 ELSE
285 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
286 $ LDB, INFO )
287 RANK = RANK + 1
288 END IF
289 100 CONTINUE
290 *
291 * Since B is complex, the following call to DGEMM is performed
292 * in two steps (real and imaginary parts). That is for V * B
293 * (in the real version of the code V**H is stored in WORK).
294 *
295 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
296 * $ WORK( NWORK ), N )
297 *
298 J = IRWB - 1
299 DO 120 JCOL = 1, NRHS
300 DO 110 JROW = 1, N
301 J = J + 1
302 RWORK( J ) = DBLE( B( JROW, JCOL ) )
303 110 CONTINUE
304 120 CONTINUE
305 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
306 $ RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
307 J = IRWB - 1
308 DO 140 JCOL = 1, NRHS
309 DO 130 JROW = 1, N
310 J = J + 1
311 RWORK( J ) = DIMAG( B( JROW, JCOL ) )
312 130 CONTINUE
313 140 CONTINUE
314 CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
315 $ RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
316 JREAL = IRWRB - 1
317 JIMAG = IRWIB - 1
318 DO 160 JCOL = 1, NRHS
319 DO 150 JROW = 1, N
320 JREAL = JREAL + 1
321 JIMAG = JIMAG + 1
322 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
323 $ RWORK( JIMAG ) )
324 150 CONTINUE
325 160 CONTINUE
326 *
327 * Unscale.
328 *
329 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
330 CALL DLASRT( 'D', N, D, INFO )
331 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
332 *
333 RETURN
334 END IF
335 *
336 * Book-keeping and setting up some constants.
337 *
338 NLVL = INT( LOG( DBLE( N ) / DBLE( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
339 *
340 SMLSZP = SMLSIZ + 1
341 *
342 U = 1
343 VT = 1 + SMLSIZ*N
344 DIFL = VT + SMLSZP*N
345 DIFR = DIFL + NLVL*N
346 Z = DIFR + NLVL*N*2
347 C = Z + NLVL*N
348 S = C + N
349 POLES = S + N
350 GIVNUM = POLES + 2*NLVL*N
351 NRWORK = GIVNUM + 2*NLVL*N
352 BX = 1
353 *
354 IRWRB = NRWORK
355 IRWIB = IRWRB + SMLSIZ*NRHS
356 IRWB = IRWIB + SMLSIZ*NRHS
357 *
358 SIZEI = 1 + N
359 K = SIZEI + N
360 GIVPTR = K + N
361 PERM = GIVPTR + N
362 GIVCOL = PERM + NLVL*N
363 IWK = GIVCOL + NLVL*N*2
364 *
365 ST = 1
366 SQRE = 0
367 ICMPQ1 = 1
368 ICMPQ2 = 0
369 NSUB = 0
370 *
371 DO 170 I = 1, N
372 IF( ABS( D( I ) ).LT.EPS ) THEN
373 D( I ) = SIGN( EPS, D( I ) )
374 END IF
375 170 CONTINUE
376 *
377 DO 240 I = 1, NM1
378 IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
379 NSUB = NSUB + 1
380 IWORK( NSUB ) = ST
381 *
382 * Subproblem found. First determine its size and then
383 * apply divide and conquer on it.
384 *
385 IF( I.LT.NM1 ) THEN
386 *
387 * A subproblem with E(I) small for I < NM1.
388 *
389 NSIZE = I - ST + 1
390 IWORK( SIZEI+NSUB-1 ) = NSIZE
391 ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
392 *
393 * A subproblem with E(NM1) not too small but I = NM1.
394 *
395 NSIZE = N - ST + 1
396 IWORK( SIZEI+NSUB-1 ) = NSIZE
397 ELSE
398 *
399 * A subproblem with E(NM1) small. This implies an
400 * 1-by-1 subproblem at D(N), which is not solved
401 * explicitly.
402 *
403 NSIZE = I - ST + 1
404 IWORK( SIZEI+NSUB-1 ) = NSIZE
405 NSUB = NSUB + 1
406 IWORK( NSUB ) = N
407 IWORK( SIZEI+NSUB-1 ) = 1
408 CALL ZCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
409 END IF
410 ST1 = ST - 1
411 IF( NSIZE.EQ.1 ) THEN
412 *
413 * This is a 1-by-1 subproblem and is not solved
414 * explicitly.
415 *
416 CALL ZCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
417 ELSE IF( NSIZE.LE.SMLSIZ ) THEN
418 *
419 * This is a small subproblem and is solved by DLASDQ.
420 *
421 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
422 $ RWORK( VT+ST1 ), N )
423 CALL DLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
424 $ RWORK( U+ST1 ), N )
425 CALL DLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
426 $ E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
427 $ N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
428 $ INFO )
429 IF( INFO.NE.0 ) THEN
430 RETURN
431 END IF
432 *
433 * In the real version, B is passed to DLASDQ and multiplied
434 * internally by Q**H. Here B is complex and that product is
435 * computed below in two steps (real and imaginary parts).
436 *
437 J = IRWB - 1
438 DO 190 JCOL = 1, NRHS
439 DO 180 JROW = ST, ST + NSIZE - 1
440 J = J + 1
441 RWORK( J ) = DBLE( B( JROW, JCOL ) )
442 180 CONTINUE
443 190 CONTINUE
444 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
445 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
446 $ ZERO, RWORK( IRWRB ), NSIZE )
447 J = IRWB - 1
448 DO 210 JCOL = 1, NRHS
449 DO 200 JROW = ST, ST + NSIZE - 1
450 J = J + 1
451 RWORK( J ) = DIMAG( B( JROW, JCOL ) )
452 200 CONTINUE
453 210 CONTINUE
454 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
455 $ RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
456 $ ZERO, RWORK( IRWIB ), NSIZE )
457 JREAL = IRWRB - 1
458 JIMAG = IRWIB - 1
459 DO 230 JCOL = 1, NRHS
460 DO 220 JROW = ST, ST + NSIZE - 1
461 JREAL = JREAL + 1
462 JIMAG = JIMAG + 1
463 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
464 $ RWORK( JIMAG ) )
465 220 CONTINUE
466 230 CONTINUE
467 *
468 CALL ZLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
469 $ WORK( BX+ST1 ), N )
470 ELSE
471 *
472 * A large problem. Solve it using divide and conquer.
473 *
474 CALL DLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
475 $ E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
476 $ IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
477 $ RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
478 $ RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
479 $ IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
480 $ RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
481 $ RWORK( S+ST1 ), RWORK( NRWORK ),
482 $ IWORK( IWK ), INFO )
483 IF( INFO.NE.0 ) THEN
484 RETURN
485 END IF
486 BXST = BX + ST1
487 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
488 $ LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
489 $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
490 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
491 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
492 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
493 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
494 $ RWORK( C+ST1 ), RWORK( S+ST1 ),
495 $ RWORK( NRWORK ), IWORK( IWK ), INFO )
496 IF( INFO.NE.0 ) THEN
497 RETURN
498 END IF
499 END IF
500 ST = I + 1
501 END IF
502 240 CONTINUE
503 *
504 * Apply the singular values and treat the tiny ones as zero.
505 *
506 TOL = RCND*ABS( D( IDAMAX( N, D, 1 ) ) )
507 *
508 DO 250 I = 1, N
509 *
510 * Some of the elements in D can be negative because 1-by-1
511 * subproblems were not solved explicitly.
512 *
513 IF( ABS( D( I ) ).LE.TOL ) THEN
514 CALL ZLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
515 ELSE
516 RANK = RANK + 1
517 CALL ZLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
518 $ WORK( BX+I-1 ), N, INFO )
519 END IF
520 D( I ) = ABS( D( I ) )
521 250 CONTINUE
522 *
523 * Now apply back the right singular vectors.
524 *
525 ICMPQ2 = 1
526 DO 320 I = 1, NSUB
527 ST = IWORK( I )
528 ST1 = ST - 1
529 NSIZE = IWORK( SIZEI+I-1 )
530 BXST = BX + ST1
531 IF( NSIZE.EQ.1 ) THEN
532 CALL ZCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
533 ELSE IF( NSIZE.LE.SMLSIZ ) THEN
534 *
535 * Since B and BX are complex, the following call to DGEMM
536 * is performed in two steps (real and imaginary parts).
537 *
538 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
539 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
540 * $ B( ST, 1 ), LDB )
541 *
542 J = BXST - N - 1
543 JREAL = IRWB - 1
544 DO 270 JCOL = 1, NRHS
545 J = J + N
546 DO 260 JROW = 1, NSIZE
547 JREAL = JREAL + 1
548 RWORK( JREAL ) = DBLE( WORK( J+JROW ) )
549 260 CONTINUE
550 270 CONTINUE
551 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
552 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
553 $ RWORK( IRWRB ), NSIZE )
554 J = BXST - N - 1
555 JIMAG = IRWB - 1
556 DO 290 JCOL = 1, NRHS
557 J = J + N
558 DO 280 JROW = 1, NSIZE
559 JIMAG = JIMAG + 1
560 RWORK( JIMAG ) = DIMAG( WORK( J+JROW ) )
561 280 CONTINUE
562 290 CONTINUE
563 CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
564 $ RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
565 $ RWORK( IRWIB ), NSIZE )
566 JREAL = IRWRB - 1
567 JIMAG = IRWIB - 1
568 DO 310 JCOL = 1, NRHS
569 DO 300 JROW = ST, ST + NSIZE - 1
570 JREAL = JREAL + 1
571 JIMAG = JIMAG + 1
572 B( JROW, JCOL ) = DCMPLX( RWORK( JREAL ),
573 $ RWORK( JIMAG ) )
574 300 CONTINUE
575 310 CONTINUE
576 ELSE
577 CALL ZLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
578 $ B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
579 $ RWORK( VT+ST1 ), IWORK( K+ST1 ),
580 $ RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
581 $ RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
582 $ IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
583 $ IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
584 $ RWORK( C+ST1 ), RWORK( S+ST1 ),
585 $ RWORK( NRWORK ), IWORK( IWK ), INFO )
586 IF( INFO.NE.0 ) THEN
587 RETURN
588 END IF
589 END IF
590 320 CONTINUE
591 *
592 * Unscale and sort the singular values.
593 *
594 CALL DLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
595 CALL DLASRT( 'D', N, D, INFO )
596 CALL ZLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
597 *
598 RETURN
599 *
600 * End of ZLALSD
601 *
602 END