LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zlalsd.f
Go to the documentation of this file.
1 *> \brief \b ZLALSD uses the singular value decomposition of A to solve the least squares problem.
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download ZLALSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zlalsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zlalsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zlalsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
22 * RANK, WORK, RWORK, IWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER UPLO
26 * INTEGER INFO, LDB, N, NRHS, RANK, SMLSIZ
27 * DOUBLE PRECISION RCOND
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * DOUBLE PRECISION D( * ), E( * ), RWORK( * )
32 * COMPLEX*16 B( LDB, * ), WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> ZLALSD uses the singular value decomposition of A to solve the least
42 *> squares problem of finding X to minimize the Euclidean norm of each
43 *> column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
44 *> are N-by-NRHS. The solution X overwrites B.
45 *>
46 *> The singular values of A smaller than RCOND times the largest
47 *> singular value are treated as zero in solving the least squares
48 *> problem; in this case a minimum norm solution is returned.
49 *> The actual singular values are returned in D in ascending order.
50 *>
51 *> This code makes very mild assumptions about floating point
52 *> arithmetic. It will work on machines with a guard digit in
53 *> add/subtract, or on those binary machines without guard digits
54 *> which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
55 *> It could conceivably fail on hexadecimal or decimal machines
56 *> without guard digits, but we know of none.
57 *> \endverbatim
58 *
59 * Arguments:
60 * ==========
61 *
62 *> \param[in] UPLO
63 *> \verbatim
64 *> UPLO is CHARACTER*1
65 *> = 'U': D and E define an upper bidiagonal matrix.
66 *> = 'L': D and E define a lower bidiagonal matrix.
67 *> \endverbatim
68 *>
69 *> \param[in] SMLSIZ
70 *> \verbatim
71 *> SMLSIZ is INTEGER
72 *> The maximum size of the subproblems at the bottom of the
73 *> computation tree.
74 *> \endverbatim
75 *>
76 *> \param[in] N
77 *> \verbatim
78 *> N is INTEGER
79 *> The dimension of the bidiagonal matrix. N >= 0.
80 *> \endverbatim
81 *>
82 *> \param[in] NRHS
83 *> \verbatim
84 *> NRHS is INTEGER
85 *> The number of columns of B. NRHS must be at least 1.
86 *> \endverbatim
87 *>
88 *> \param[in,out] D
89 *> \verbatim
90 *> D is DOUBLE PRECISION array, dimension (N)
91 *> On entry D contains the main diagonal of the bidiagonal
92 *> matrix. On exit, if INFO = 0, D contains its singular values.
93 *> \endverbatim
94 *>
95 *> \param[in,out] E
96 *> \verbatim
97 *> E is DOUBLE PRECISION array, dimension (N-1)
98 *> Contains the super-diagonal entries of the bidiagonal matrix.
99 *> On exit, E has been destroyed.
100 *> \endverbatim
101 *>
102 *> \param[in,out] B
103 *> \verbatim
104 *> B is COMPLEX*16 array, dimension (LDB,NRHS)
105 *> On input, B contains the right hand sides of the least
106 *> squares problem. On output, B contains the solution X.
107 *> \endverbatim
108 *>
109 *> \param[in] LDB
110 *> \verbatim
111 *> LDB is INTEGER
112 *> The leading dimension of B in the calling subprogram.
113 *> LDB must be at least max(1,N).
114 *> \endverbatim
115 *>
116 *> \param[in] RCOND
117 *> \verbatim
118 *> RCOND is DOUBLE PRECISION
119 *> The singular values of A less than or equal to RCOND times
120 *> the largest singular value are treated as zero in solving
121 *> the least squares problem. If RCOND is negative,
122 *> machine precision is used instead.
123 *> For example, if diag(S)*X=B were the least squares problem,
124 *> where diag(S) is a diagonal matrix of singular values, the
125 *> solution would be X(i) = B(i) / S(i) if S(i) is greater than
126 *> RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
127 *> RCOND*max(S).
128 *> \endverbatim
129 *>
130 *> \param[out] RANK
131 *> \verbatim
132 *> RANK is INTEGER
133 *> The number of singular values of A greater than RCOND times
134 *> the largest singular value.
135 *> \endverbatim
136 *>
137 *> \param[out] WORK
138 *> \verbatim
139 *> WORK is COMPLEX*16 array, dimension at least
140 *> (N * NRHS).
141 *> \endverbatim
142 *>
143 *> \param[out] RWORK
144 *> \verbatim
145 *> RWORK is DOUBLE PRECISION array, dimension at least
146 *> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
147 *> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
148 *> where
149 *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
150 *> \endverbatim
151 *>
152 *> \param[out] IWORK
153 *> \verbatim
154 *> IWORK is INTEGER array, dimension at least
155 *> (3*N*NLVL + 11*N).
156 *> \endverbatim
157 *>
158 *> \param[out] INFO
159 *> \verbatim
160 *> INFO is INTEGER
161 *> = 0: successful exit.
162 *> < 0: if INFO = -i, the i-th argument had an illegal value.
163 *> > 0: The algorithm failed to compute a singular value while
164 *> working on the submatrix lying in rows and columns
165 *> INFO/(N+1) through MOD(INFO,N+1).
166 *> \endverbatim
167 *
168 * Authors:
169 * ========
170 *
171 *> \author Univ. of Tennessee
172 *> \author Univ. of California Berkeley
173 *> \author Univ. of Colorado Denver
174 *> \author NAG Ltd.
175 *
176 *> \date September 2012
177 *
178 *> \ingroup complex16OTHERcomputational
179 *
180 *> \par Contributors:
181 * ==================
182 *>
183 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
184 *> California at Berkeley, USA \n
185 *> Osni Marques, LBNL/NERSC, USA \n
186 *
187 * =====================================================================
188  SUBROUTINE zlalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
189  $ rank, work, rwork, iwork, info )
190 *
191 * -- LAPACK computational routine (version 3.4.2) --
192 * -- LAPACK is a software package provided by Univ. of Tennessee, --
193 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
194 * September 2012
195 *
196 * .. Scalar Arguments ..
197  CHARACTER uplo
198  INTEGER info, ldb, n, nrhs, rank, smlsiz
199  DOUBLE PRECISION rcond
200 * ..
201 * .. Array Arguments ..
202  INTEGER iwork( * )
203  DOUBLE PRECISION d( * ), e( * ), rwork( * )
204  COMPLEX*16 b( ldb, * ), work( * )
205 * ..
206 *
207 * =====================================================================
208 *
209 * .. Parameters ..
210  DOUBLE PRECISION zero, one, two
211  parameter( zero = 0.0d0, one = 1.0d0, two = 2.0d0 )
212  COMPLEX*16 czero
213  parameter( czero = ( 0.0d0, 0.0d0 ) )
214 * ..
215 * .. Local Scalars ..
216  INTEGER bx, bxst, c, difl, difr, givcol, givnum,
217  $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
218  $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
219  $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
220  $ perm, poles, s, sizei, smlszp, sqre, st, st1,
221  $ u, vt, z
222  DOUBLE PRECISION cs, eps, orgnrm, rcnd, r, sn, tol
223 * ..
224 * .. External Functions ..
225  INTEGER idamax
226  DOUBLE PRECISION dlamch, dlanst
227  EXTERNAL idamax, dlamch, dlanst
228 * ..
229 * .. External Subroutines ..
230  EXTERNAL dgemm, dlartg, dlascl, dlasda, dlasdq, dlaset,
231  $ dlasrt, xerbla, zcopy, zdrot, zlacpy, zlalsa,
232  $ zlascl, zlaset
233 * ..
234 * .. Intrinsic Functions ..
235  INTRINSIC abs, dble, dcmplx, dimag, int, log, sign
236 * ..
237 * .. Executable Statements ..
238 *
239 * Test the input parameters.
240 *
241  info = 0
242 *
243  IF( n.LT.0 ) THEN
244  info = -3
245  ELSE IF( nrhs.LT.1 ) THEN
246  info = -4
247  ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
248  info = -8
249  END IF
250  IF( info.NE.0 ) THEN
251  CALL xerbla( 'ZLALSD', -info )
252  return
253  END IF
254 *
255  eps = dlamch( 'Epsilon' )
256 *
257 * Set up the tolerance.
258 *
259  IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
260  rcnd = eps
261  ELSE
262  rcnd = rcond
263  END IF
264 *
265  rank = 0
266 *
267 * Quick return if possible.
268 *
269  IF( n.EQ.0 ) THEN
270  return
271  ELSE IF( n.EQ.1 ) THEN
272  IF( d( 1 ).EQ.zero ) THEN
273  CALL zlaset( 'A', 1, nrhs, czero, czero, b, ldb )
274  ELSE
275  rank = 1
276  CALL zlascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
277  d( 1 ) = abs( d( 1 ) )
278  END IF
279  return
280  END IF
281 *
282 * Rotate the matrix if it is lower bidiagonal.
283 *
284  IF( uplo.EQ.'L' ) THEN
285  DO 10 i = 1, n - 1
286  CALL dlartg( d( i ), e( i ), cs, sn, r )
287  d( i ) = r
288  e( i ) = sn*d( i+1 )
289  d( i+1 ) = cs*d( i+1 )
290  IF( nrhs.EQ.1 ) THEN
291  CALL zdrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
292  ELSE
293  rwork( i*2-1 ) = cs
294  rwork( i*2 ) = sn
295  END IF
296  10 continue
297  IF( nrhs.GT.1 ) THEN
298  DO 30 i = 1, nrhs
299  DO 20 j = 1, n - 1
300  cs = rwork( j*2-1 )
301  sn = rwork( j*2 )
302  CALL zdrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
303  20 continue
304  30 continue
305  END IF
306  END IF
307 *
308 * Scale.
309 *
310  nm1 = n - 1
311  orgnrm = dlanst( 'M', n, d, e )
312  IF( orgnrm.EQ.zero ) THEN
313  CALL zlaset( 'A', n, nrhs, czero, czero, b, ldb )
314  return
315  END IF
316 *
317  CALL dlascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
318  CALL dlascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
319 *
320 * If N is smaller than the minimum divide size SMLSIZ, then solve
321 * the problem with another solver.
322 *
323  IF( n.LE.smlsiz ) THEN
324  irwu = 1
325  irwvt = irwu + n*n
326  irwwrk = irwvt + n*n
327  irwrb = irwwrk
328  irwib = irwrb + n*nrhs
329  irwb = irwib + n*nrhs
330  CALL dlaset( 'A', n, n, zero, one, rwork( irwu ), n )
331  CALL dlaset( 'A', n, n, zero, one, rwork( irwvt ), n )
332  CALL dlasdq( 'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
333  $ rwork( irwu ), n, rwork( irwwrk ), 1,
334  $ rwork( irwwrk ), info )
335  IF( info.NE.0 ) THEN
336  return
337  END IF
338 *
339 * In the real version, B is passed to DLASDQ and multiplied
340 * internally by Q**H. Here B is complex and that product is
341 * computed below in two steps (real and imaginary parts).
342 *
343  j = irwb - 1
344  DO 50 jcol = 1, nrhs
345  DO 40 jrow = 1, n
346  j = j + 1
347  rwork( j ) = dble( b( jrow, jcol ) )
348  40 continue
349  50 continue
350  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
351  $ rwork( irwb ), n, zero, rwork( irwrb ), n )
352  j = irwb - 1
353  DO 70 jcol = 1, nrhs
354  DO 60 jrow = 1, n
355  j = j + 1
356  rwork( j ) = dimag( b( jrow, jcol ) )
357  60 continue
358  70 continue
359  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
360  $ rwork( irwb ), n, zero, rwork( irwib ), n )
361  jreal = irwrb - 1
362  jimag = irwib - 1
363  DO 90 jcol = 1, nrhs
364  DO 80 jrow = 1, n
365  jreal = jreal + 1
366  jimag = jimag + 1
367  b( jrow, jcol ) = dcmplx( rwork( jreal ),
368  $ rwork( jimag ) )
369  80 continue
370  90 continue
371 *
372  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
373  DO 100 i = 1, n
374  IF( d( i ).LE.tol ) THEN
375  CALL zlaset( 'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
376  ELSE
377  CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
378  $ ldb, info )
379  rank = rank + 1
380  END IF
381  100 continue
382 *
383 * Since B is complex, the following call to DGEMM is performed
384 * in two steps (real and imaginary parts). That is for V * B
385 * (in the real version of the code V**H is stored in WORK).
386 *
387 * CALL DGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
388 * $ WORK( NWORK ), N )
389 *
390  j = irwb - 1
391  DO 120 jcol = 1, nrhs
392  DO 110 jrow = 1, n
393  j = j + 1
394  rwork( j ) = dble( b( jrow, jcol ) )
395  110 continue
396  120 continue
397  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
398  $ rwork( irwb ), n, zero, rwork( irwrb ), n )
399  j = irwb - 1
400  DO 140 jcol = 1, nrhs
401  DO 130 jrow = 1, n
402  j = j + 1
403  rwork( j ) = dimag( b( jrow, jcol ) )
404  130 continue
405  140 continue
406  CALL dgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
407  $ rwork( irwb ), n, zero, rwork( irwib ), n )
408  jreal = irwrb - 1
409  jimag = irwib - 1
410  DO 160 jcol = 1, nrhs
411  DO 150 jrow = 1, n
412  jreal = jreal + 1
413  jimag = jimag + 1
414  b( jrow, jcol ) = dcmplx( rwork( jreal ),
415  $ rwork( jimag ) )
416  150 continue
417  160 continue
418 *
419 * Unscale.
420 *
421  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
422  CALL dlasrt( 'D', n, d, info )
423  CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
424 *
425  return
426  END IF
427 *
428 * Book-keeping and setting up some constants.
429 *
430  nlvl = int( log( dble( n ) / dble( smlsiz+1 ) ) / log( two ) ) + 1
431 *
432  smlszp = smlsiz + 1
433 *
434  u = 1
435  vt = 1 + smlsiz*n
436  difl = vt + smlszp*n
437  difr = difl + nlvl*n
438  z = difr + nlvl*n*2
439  c = z + nlvl*n
440  s = c + n
441  poles = s + n
442  givnum = poles + 2*nlvl*n
443  nrwork = givnum + 2*nlvl*n
444  bx = 1
445 *
446  irwrb = nrwork
447  irwib = irwrb + smlsiz*nrhs
448  irwb = irwib + smlsiz*nrhs
449 *
450  sizei = 1 + n
451  k = sizei + n
452  givptr = k + n
453  perm = givptr + n
454  givcol = perm + nlvl*n
455  iwk = givcol + nlvl*n*2
456 *
457  st = 1
458  sqre = 0
459  icmpq1 = 1
460  icmpq2 = 0
461  nsub = 0
462 *
463  DO 170 i = 1, n
464  IF( abs( d( i ) ).LT.eps ) THEN
465  d( i ) = sign( eps, d( i ) )
466  END IF
467  170 continue
468 *
469  DO 240 i = 1, nm1
470  IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
471  nsub = nsub + 1
472  iwork( nsub ) = st
473 *
474 * Subproblem found. First determine its size and then
475 * apply divide and conquer on it.
476 *
477  IF( i.LT.nm1 ) THEN
478 *
479 * A subproblem with E(I) small for I < NM1.
480 *
481  nsize = i - st + 1
482  iwork( sizei+nsub-1 ) = nsize
483  ELSE IF( abs( e( i ) ).GE.eps ) THEN
484 *
485 * A subproblem with E(NM1) not too small but I = NM1.
486 *
487  nsize = n - st + 1
488  iwork( sizei+nsub-1 ) = nsize
489  ELSE
490 *
491 * A subproblem with E(NM1) small. This implies an
492 * 1-by-1 subproblem at D(N), which is not solved
493 * explicitly.
494 *
495  nsize = i - st + 1
496  iwork( sizei+nsub-1 ) = nsize
497  nsub = nsub + 1
498  iwork( nsub ) = n
499  iwork( sizei+nsub-1 ) = 1
500  CALL zcopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
501  END IF
502  st1 = st - 1
503  IF( nsize.EQ.1 ) THEN
504 *
505 * This is a 1-by-1 subproblem and is not solved
506 * explicitly.
507 *
508  CALL zcopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
509  ELSE IF( nsize.LE.smlsiz ) THEN
510 *
511 * This is a small subproblem and is solved by DLASDQ.
512 *
513  CALL dlaset( 'A', nsize, nsize, zero, one,
514  $ rwork( vt+st1 ), n )
515  CALL dlaset( 'A', nsize, nsize, zero, one,
516  $ rwork( u+st1 ), n )
517  CALL dlasdq( 'U', 0, nsize, nsize, nsize, 0, d( st ),
518  $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
519  $ n, rwork( nrwork ), 1, rwork( nrwork ),
520  $ info )
521  IF( info.NE.0 ) THEN
522  return
523  END IF
524 *
525 * In the real version, B is passed to DLASDQ and multiplied
526 * internally by Q**H. Here B is complex and that product is
527 * computed below in two steps (real and imaginary parts).
528 *
529  j = irwb - 1
530  DO 190 jcol = 1, nrhs
531  DO 180 jrow = st, st + nsize - 1
532  j = j + 1
533  rwork( j ) = dble( b( jrow, jcol ) )
534  180 continue
535  190 continue
536  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
537  $ rwork( u+st1 ), n, rwork( irwb ), nsize,
538  $ zero, rwork( irwrb ), nsize )
539  j = irwb - 1
540  DO 210 jcol = 1, nrhs
541  DO 200 jrow = st, st + nsize - 1
542  j = j + 1
543  rwork( j ) = dimag( b( jrow, jcol ) )
544  200 continue
545  210 continue
546  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
547  $ rwork( u+st1 ), n, rwork( irwb ), nsize,
548  $ zero, rwork( irwib ), nsize )
549  jreal = irwrb - 1
550  jimag = irwib - 1
551  DO 230 jcol = 1, nrhs
552  DO 220 jrow = st, st + nsize - 1
553  jreal = jreal + 1
554  jimag = jimag + 1
555  b( jrow, jcol ) = dcmplx( rwork( jreal ),
556  $ rwork( jimag ) )
557  220 continue
558  230 continue
559 *
560  CALL zlacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
561  $ work( bx+st1 ), n )
562  ELSE
563 *
564 * A large problem. Solve it using divide and conquer.
565 *
566  CALL dlasda( icmpq1, smlsiz, nsize, sqre, d( st ),
567  $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
568  $ iwork( k+st1 ), rwork( difl+st1 ),
569  $ rwork( difr+st1 ), rwork( z+st1 ),
570  $ rwork( poles+st1 ), iwork( givptr+st1 ),
571  $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
572  $ rwork( givnum+st1 ), rwork( c+st1 ),
573  $ rwork( s+st1 ), rwork( nrwork ),
574  $ iwork( iwk ), info )
575  IF( info.NE.0 ) THEN
576  return
577  END IF
578  bxst = bx + st1
579  CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
580  $ ldb, work( bxst ), n, rwork( u+st1 ), n,
581  $ rwork( vt+st1 ), iwork( k+st1 ),
582  $ rwork( difl+st1 ), rwork( difr+st1 ),
583  $ rwork( z+st1 ), rwork( poles+st1 ),
584  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
585  $ iwork( perm+st1 ), rwork( givnum+st1 ),
586  $ rwork( c+st1 ), rwork( s+st1 ),
587  $ rwork( nrwork ), iwork( iwk ), info )
588  IF( info.NE.0 ) THEN
589  return
590  END IF
591  END IF
592  st = i + 1
593  END IF
594  240 continue
595 *
596 * Apply the singular values and treat the tiny ones as zero.
597 *
598  tol = rcnd*abs( d( idamax( n, d, 1 ) ) )
599 *
600  DO 250 i = 1, n
601 *
602 * Some of the elements in D can be negative because 1-by-1
603 * subproblems were not solved explicitly.
604 *
605  IF( abs( d( i ) ).LE.tol ) THEN
606  CALL zlaset( 'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
607  ELSE
608  rank = rank + 1
609  CALL zlascl( 'G', 0, 0, d( i ), one, 1, nrhs,
610  $ work( bx+i-1 ), n, info )
611  END IF
612  d( i ) = abs( d( i ) )
613  250 continue
614 *
615 * Now apply back the right singular vectors.
616 *
617  icmpq2 = 1
618  DO 320 i = 1, nsub
619  st = iwork( i )
620  st1 = st - 1
621  nsize = iwork( sizei+i-1 )
622  bxst = bx + st1
623  IF( nsize.EQ.1 ) THEN
624  CALL zcopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
625  ELSE IF( nsize.LE.smlsiz ) THEN
626 *
627 * Since B and BX are complex, the following call to DGEMM
628 * is performed in two steps (real and imaginary parts).
629 *
630 * CALL DGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
631 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
632 * $ B( ST, 1 ), LDB )
633 *
634  j = bxst - n - 1
635  jreal = irwb - 1
636  DO 270 jcol = 1, nrhs
637  j = j + n
638  DO 260 jrow = 1, nsize
639  jreal = jreal + 1
640  rwork( jreal ) = dble( work( j+jrow ) )
641  260 continue
642  270 continue
643  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
644  $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
645  $ rwork( irwrb ), nsize )
646  j = bxst - n - 1
647  jimag = irwb - 1
648  DO 290 jcol = 1, nrhs
649  j = j + n
650  DO 280 jrow = 1, nsize
651  jimag = jimag + 1
652  rwork( jimag ) = dimag( work( j+jrow ) )
653  280 continue
654  290 continue
655  CALL dgemm( 'T', 'N', nsize, nrhs, nsize, one,
656  $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
657  $ rwork( irwib ), nsize )
658  jreal = irwrb - 1
659  jimag = irwib - 1
660  DO 310 jcol = 1, nrhs
661  DO 300 jrow = st, st + nsize - 1
662  jreal = jreal + 1
663  jimag = jimag + 1
664  b( jrow, jcol ) = dcmplx( rwork( jreal ),
665  $ rwork( jimag ) )
666  300 continue
667  310 continue
668  ELSE
669  CALL zlalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
670  $ b( st, 1 ), ldb, rwork( u+st1 ), n,
671  $ rwork( vt+st1 ), iwork( k+st1 ),
672  $ rwork( difl+st1 ), rwork( difr+st1 ),
673  $ rwork( z+st1 ), rwork( poles+st1 ),
674  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
675  $ iwork( perm+st1 ), rwork( givnum+st1 ),
676  $ rwork( c+st1 ), rwork( s+st1 ),
677  $ rwork( nrwork ), iwork( iwk ), info )
678  IF( info.NE.0 ) THEN
679  return
680  END IF
681  END IF
682  320 continue
683 *
684 * Unscale and sort the singular values.
685 *
686  CALL dlascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
687  CALL dlasrt( 'D', n, d, info )
688  CALL zlascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
689 *
690  return
691 *
692 * End of ZLALSD
693 *
694  END