LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
clalsd.f
Go to the documentation of this file.
1 *> \brief \b CLALSD 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 CLALSD + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/clalsd.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/clalsd.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/clalsd.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE CLALSD( 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 * REAL RCOND
28 * ..
29 * .. Array Arguments ..
30 * INTEGER IWORK( * )
31 * REAL D( * ), E( * ), RWORK( * )
32 * COMPLEX B( LDB, * ), WORK( * )
33 * ..
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> CLALSD 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 REAL 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 REAL 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 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 REAL
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 array, dimension (N * NRHS).
140 *> \endverbatim
141 *>
142 *> \param[out] RWORK
143 *> \verbatim
144 *> RWORK is REAL array, dimension at least
145 *> (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS +
146 *> MAX( (SMLSIZ+1)**2, N*(1+NRHS) + 2*NRHS ),
147 *> where
148 *> NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
149 *> \endverbatim
150 *>
151 *> \param[out] IWORK
152 *> \verbatim
153 *> IWORK is INTEGER array, dimension (3*N*NLVL + 11*N).
154 *> \endverbatim
155 *>
156 *> \param[out] INFO
157 *> \verbatim
158 *> INFO is INTEGER
159 *> = 0: successful exit.
160 *> < 0: if INFO = -i, the i-th argument had an illegal value.
161 *> > 0: The algorithm failed to compute a singular value while
162 *> working on the submatrix lying in rows and columns
163 *> INFO/(N+1) through MOD(INFO,N+1).
164 *> \endverbatim
165 *
166 * Authors:
167 * ========
168 *
169 *> \author Univ. of Tennessee
170 *> \author Univ. of California Berkeley
171 *> \author Univ. of Colorado Denver
172 *> \author NAG Ltd.
173 *
174 *> \date September 2012
175 *
176 *> \ingroup complexOTHERcomputational
177 *
178 *> \par Contributors:
179 * ==================
180 *>
181 *> Ming Gu and Ren-Cang Li, Computer Science Division, University of
182 *> California at Berkeley, USA \n
183 *> Osni Marques, LBNL/NERSC, USA \n
184 *
185 * =====================================================================
186  SUBROUTINE clalsd( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
187  $ rank, work, rwork, iwork, info )
188 *
189 * -- LAPACK computational routine (version 3.4.2) --
190 * -- LAPACK is a software package provided by Univ. of Tennessee, --
191 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
192 * September 2012
193 *
194 * .. Scalar Arguments ..
195  CHARACTER uplo
196  INTEGER info, ldb, n, nrhs, rank, smlsiz
197  REAL rcond
198 * ..
199 * .. Array Arguments ..
200  INTEGER iwork( * )
201  REAL d( * ), e( * ), rwork( * )
202  COMPLEX b( ldb, * ), work( * )
203 * ..
204 *
205 * =====================================================================
206 *
207 * .. Parameters ..
208  REAL zero, one, two
209  parameter( zero = 0.0e0, one = 1.0e0, two = 2.0e0 )
210  COMPLEX czero
211  parameter( czero = ( 0.0e0, 0.0e0 ) )
212 * ..
213 * .. Local Scalars ..
214  INTEGER bx, bxst, c, difl, difr, givcol, givnum,
215  $ givptr, i, icmpq1, icmpq2, irwb, irwib, irwrb,
216  $ irwu, irwvt, irwwrk, iwk, j, jcol, jimag,
217  $ jreal, jrow, k, nlvl, nm1, nrwork, nsize, nsub,
218  $ perm, poles, s, sizei, smlszp, sqre, st, st1,
219  $ u, vt, z
220  REAL cs, eps, orgnrm, r, rcnd, sn, tol
221 * ..
222 * .. External Functions ..
223  INTEGER isamax
224  REAL slamch, slanst
225  EXTERNAL isamax, slamch, slanst
226 * ..
227 * .. External Subroutines ..
228  EXTERNAL ccopy, clacpy, clalsa, clascl, claset, csrot,
229  $ sgemm, slartg, slascl, slasda, slasdq, slaset,
230  $ slasrt, xerbla
231 * ..
232 * .. Intrinsic Functions ..
233  INTRINSIC abs, aimag, cmplx, int, log, REAL, sign
234 * ..
235 * .. Executable Statements ..
236 *
237 * Test the input parameters.
238 *
239  info = 0
240 *
241  IF( n.LT.0 ) THEN
242  info = -3
243  ELSE IF( nrhs.LT.1 ) THEN
244  info = -4
245  ELSE IF( ( ldb.LT.1 ) .OR. ( ldb.LT.n ) ) THEN
246  info = -8
247  END IF
248  IF( info.NE.0 ) THEN
249  CALL xerbla( 'CLALSD', -info )
250  return
251  END IF
252 *
253  eps = slamch( 'Epsilon' )
254 *
255 * Set up the tolerance.
256 *
257  IF( ( rcond.LE.zero ) .OR. ( rcond.GE.one ) ) THEN
258  rcnd = eps
259  ELSE
260  rcnd = rcond
261  END IF
262 *
263  rank = 0
264 *
265 * Quick return if possible.
266 *
267  IF( n.EQ.0 ) THEN
268  return
269  ELSE IF( n.EQ.1 ) THEN
270  IF( d( 1 ).EQ.zero ) THEN
271  CALL claset( 'A', 1, nrhs, czero, czero, b, ldb )
272  ELSE
273  rank = 1
274  CALL clascl( 'G', 0, 0, d( 1 ), one, 1, nrhs, b, ldb, info )
275  d( 1 ) = abs( d( 1 ) )
276  END IF
277  return
278  END IF
279 *
280 * Rotate the matrix if it is lower bidiagonal.
281 *
282  IF( uplo.EQ.'L' ) THEN
283  DO 10 i = 1, n - 1
284  CALL slartg( d( i ), e( i ), cs, sn, r )
285  d( i ) = r
286  e( i ) = sn*d( i+1 )
287  d( i+1 ) = cs*d( i+1 )
288  IF( nrhs.EQ.1 ) THEN
289  CALL csrot( 1, b( i, 1 ), 1, b( i+1, 1 ), 1, cs, sn )
290  ELSE
291  rwork( i*2-1 ) = cs
292  rwork( i*2 ) = sn
293  END IF
294  10 continue
295  IF( nrhs.GT.1 ) THEN
296  DO 30 i = 1, nrhs
297  DO 20 j = 1, n - 1
298  cs = rwork( j*2-1 )
299  sn = rwork( j*2 )
300  CALL csrot( 1, b( j, i ), 1, b( j+1, i ), 1, cs, sn )
301  20 continue
302  30 continue
303  END IF
304  END IF
305 *
306 * Scale.
307 *
308  nm1 = n - 1
309  orgnrm = slanst( 'M', n, d, e )
310  IF( orgnrm.EQ.zero ) THEN
311  CALL claset( 'A', n, nrhs, czero, czero, b, ldb )
312  return
313  END IF
314 *
315  CALL slascl( 'G', 0, 0, orgnrm, one, n, 1, d, n, info )
316  CALL slascl( 'G', 0, 0, orgnrm, one, nm1, 1, e, nm1, info )
317 *
318 * If N is smaller than the minimum divide size SMLSIZ, then solve
319 * the problem with another solver.
320 *
321  IF( n.LE.smlsiz ) THEN
322  irwu = 1
323  irwvt = irwu + n*n
324  irwwrk = irwvt + n*n
325  irwrb = irwwrk
326  irwib = irwrb + n*nrhs
327  irwb = irwib + n*nrhs
328  CALL slaset( 'A', n, n, zero, one, rwork( irwu ), n )
329  CALL slaset( 'A', n, n, zero, one, rwork( irwvt ), n )
330  CALL slasdq( 'U', 0, n, n, n, 0, d, e, rwork( irwvt ), n,
331  $ rwork( irwu ), n, rwork( irwwrk ), 1,
332  $ rwork( irwwrk ), info )
333  IF( info.NE.0 ) THEN
334  return
335  END IF
336 *
337 * In the real version, B is passed to SLASDQ and multiplied
338 * internally by Q**H. Here B is complex and that product is
339 * computed below in two steps (real and imaginary parts).
340 *
341  j = irwb - 1
342  DO 50 jcol = 1, nrhs
343  DO 40 jrow = 1, n
344  j = j + 1
345  rwork( j ) = REAL( B( JROW, JCOL ) )
346  40 continue
347  50 continue
348  CALL sgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
349  $ rwork( irwb ), n, zero, rwork( irwrb ), n )
350  j = irwb - 1
351  DO 70 jcol = 1, nrhs
352  DO 60 jrow = 1, n
353  j = j + 1
354  rwork( j ) = aimag( b( jrow, jcol ) )
355  60 continue
356  70 continue
357  CALL sgemm( 'T', 'N', n, nrhs, n, one, rwork( irwu ), n,
358  $ rwork( irwb ), n, zero, rwork( irwib ), n )
359  jreal = irwrb - 1
360  jimag = irwib - 1
361  DO 90 jcol = 1, nrhs
362  DO 80 jrow = 1, n
363  jreal = jreal + 1
364  jimag = jimag + 1
365  b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
366  80 continue
367  90 continue
368 *
369  tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
370  DO 100 i = 1, n
371  IF( d( i ).LE.tol ) THEN
372  CALL claset( 'A', 1, nrhs, czero, czero, b( i, 1 ), ldb )
373  ELSE
374  CALL clascl( 'G', 0, 0, d( i ), one, 1, nrhs, b( i, 1 ),
375  $ ldb, info )
376  rank = rank + 1
377  END IF
378  100 continue
379 *
380 * Since B is complex, the following call to SGEMM is performed
381 * in two steps (real and imaginary parts). That is for V * B
382 * (in the real version of the code V**H is stored in WORK).
383 *
384 * CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
385 * $ WORK( NWORK ), N )
386 *
387  j = irwb - 1
388  DO 120 jcol = 1, nrhs
389  DO 110 jrow = 1, n
390  j = j + 1
391  rwork( j ) = REAL( B( JROW, JCOL ) )
392  110 continue
393  120 continue
394  CALL sgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
395  $ rwork( irwb ), n, zero, rwork( irwrb ), n )
396  j = irwb - 1
397  DO 140 jcol = 1, nrhs
398  DO 130 jrow = 1, n
399  j = j + 1
400  rwork( j ) = aimag( b( jrow, jcol ) )
401  130 continue
402  140 continue
403  CALL sgemm( 'T', 'N', n, nrhs, n, one, rwork( irwvt ), n,
404  $ rwork( irwb ), n, zero, rwork( irwib ), n )
405  jreal = irwrb - 1
406  jimag = irwib - 1
407  DO 160 jcol = 1, nrhs
408  DO 150 jrow = 1, n
409  jreal = jreal + 1
410  jimag = jimag + 1
411  b( jrow, jcol ) = cmplx( rwork( jreal ), rwork( jimag ) )
412  150 continue
413  160 continue
414 *
415 * Unscale.
416 *
417  CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
418  CALL slasrt( 'D', n, d, info )
419  CALL clascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
420 *
421  return
422  END IF
423 *
424 * Book-keeping and setting up some constants.
425 *
426  nlvl = int( log( REAL( N ) / REAL( SMLSIZ+1 ) ) / log( two ) ) + 1
427 *
428  smlszp = smlsiz + 1
429 *
430  u = 1
431  vt = 1 + smlsiz*n
432  difl = vt + smlszp*n
433  difr = difl + nlvl*n
434  z = difr + nlvl*n*2
435  c = z + nlvl*n
436  s = c + n
437  poles = s + n
438  givnum = poles + 2*nlvl*n
439  nrwork = givnum + 2*nlvl*n
440  bx = 1
441 *
442  irwrb = nrwork
443  irwib = irwrb + smlsiz*nrhs
444  irwb = irwib + smlsiz*nrhs
445 *
446  sizei = 1 + n
447  k = sizei + n
448  givptr = k + n
449  perm = givptr + n
450  givcol = perm + nlvl*n
451  iwk = givcol + nlvl*n*2
452 *
453  st = 1
454  sqre = 0
455  icmpq1 = 1
456  icmpq2 = 0
457  nsub = 0
458 *
459  DO 170 i = 1, n
460  IF( abs( d( i ) ).LT.eps ) THEN
461  d( i ) = sign( eps, d( i ) )
462  END IF
463  170 continue
464 *
465  DO 240 i = 1, nm1
466  IF( ( abs( e( i ) ).LT.eps ) .OR. ( i.EQ.nm1 ) ) THEN
467  nsub = nsub + 1
468  iwork( nsub ) = st
469 *
470 * Subproblem found. First determine its size and then
471 * apply divide and conquer on it.
472 *
473  IF( i.LT.nm1 ) THEN
474 *
475 * A subproblem with E(I) small for I < NM1.
476 *
477  nsize = i - st + 1
478  iwork( sizei+nsub-1 ) = nsize
479  ELSE IF( abs( e( i ) ).GE.eps ) THEN
480 *
481 * A subproblem with E(NM1) not too small but I = NM1.
482 *
483  nsize = n - st + 1
484  iwork( sizei+nsub-1 ) = nsize
485  ELSE
486 *
487 * A subproblem with E(NM1) small. This implies an
488 * 1-by-1 subproblem at D(N), which is not solved
489 * explicitly.
490 *
491  nsize = i - st + 1
492  iwork( sizei+nsub-1 ) = nsize
493  nsub = nsub + 1
494  iwork( nsub ) = n
495  iwork( sizei+nsub-1 ) = 1
496  CALL ccopy( nrhs, b( n, 1 ), ldb, work( bx+nm1 ), n )
497  END IF
498  st1 = st - 1
499  IF( nsize.EQ.1 ) THEN
500 *
501 * This is a 1-by-1 subproblem and is not solved
502 * explicitly.
503 *
504  CALL ccopy( nrhs, b( st, 1 ), ldb, work( bx+st1 ), n )
505  ELSE IF( nsize.LE.smlsiz ) THEN
506 *
507 * This is a small subproblem and is solved by SLASDQ.
508 *
509  CALL slaset( 'A', nsize, nsize, zero, one,
510  $ rwork( vt+st1 ), n )
511  CALL slaset( 'A', nsize, nsize, zero, one,
512  $ rwork( u+st1 ), n )
513  CALL slasdq( 'U', 0, nsize, nsize, nsize, 0, d( st ),
514  $ e( st ), rwork( vt+st1 ), n, rwork( u+st1 ),
515  $ n, rwork( nrwork ), 1, rwork( nrwork ),
516  $ info )
517  IF( info.NE.0 ) THEN
518  return
519  END IF
520 *
521 * In the real version, B is passed to SLASDQ and multiplied
522 * internally by Q**H. Here B is complex and that product is
523 * computed below in two steps (real and imaginary parts).
524 *
525  j = irwb - 1
526  DO 190 jcol = 1, nrhs
527  DO 180 jrow = st, st + nsize - 1
528  j = j + 1
529  rwork( j ) = REAL( B( JROW, JCOL ) )
530  180 continue
531  190 continue
532  CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
533  $ rwork( u+st1 ), n, rwork( irwb ), nsize,
534  $ zero, rwork( irwrb ), nsize )
535  j = irwb - 1
536  DO 210 jcol = 1, nrhs
537  DO 200 jrow = st, st + nsize - 1
538  j = j + 1
539  rwork( j ) = aimag( b( jrow, jcol ) )
540  200 continue
541  210 continue
542  CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
543  $ rwork( u+st1 ), n, rwork( irwb ), nsize,
544  $ zero, rwork( irwib ), nsize )
545  jreal = irwrb - 1
546  jimag = irwib - 1
547  DO 230 jcol = 1, nrhs
548  DO 220 jrow = st, st + nsize - 1
549  jreal = jreal + 1
550  jimag = jimag + 1
551  b( jrow, jcol ) = cmplx( rwork( jreal ),
552  $ rwork( jimag ) )
553  220 continue
554  230 continue
555 *
556  CALL clacpy( 'A', nsize, nrhs, b( st, 1 ), ldb,
557  $ work( bx+st1 ), n )
558  ELSE
559 *
560 * A large problem. Solve it using divide and conquer.
561 *
562  CALL slasda( icmpq1, smlsiz, nsize, sqre, d( st ),
563  $ e( st ), rwork( u+st1 ), n, rwork( vt+st1 ),
564  $ iwork( k+st1 ), rwork( difl+st1 ),
565  $ rwork( difr+st1 ), rwork( z+st1 ),
566  $ rwork( poles+st1 ), iwork( givptr+st1 ),
567  $ iwork( givcol+st1 ), n, iwork( perm+st1 ),
568  $ rwork( givnum+st1 ), rwork( c+st1 ),
569  $ rwork( s+st1 ), rwork( nrwork ),
570  $ iwork( iwk ), info )
571  IF( info.NE.0 ) THEN
572  return
573  END IF
574  bxst = bx + st1
575  CALL clalsa( icmpq2, smlsiz, nsize, nrhs, b( st, 1 ),
576  $ ldb, work( bxst ), n, rwork( u+st1 ), n,
577  $ rwork( vt+st1 ), iwork( k+st1 ),
578  $ rwork( difl+st1 ), rwork( difr+st1 ),
579  $ rwork( z+st1 ), rwork( poles+st1 ),
580  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
581  $ iwork( perm+st1 ), rwork( givnum+st1 ),
582  $ rwork( c+st1 ), rwork( s+st1 ),
583  $ rwork( nrwork ), iwork( iwk ), info )
584  IF( info.NE.0 ) THEN
585  return
586  END IF
587  END IF
588  st = i + 1
589  END IF
590  240 continue
591 *
592 * Apply the singular values and treat the tiny ones as zero.
593 *
594  tol = rcnd*abs( d( isamax( n, d, 1 ) ) )
595 *
596  DO 250 i = 1, n
597 *
598 * Some of the elements in D can be negative because 1-by-1
599 * subproblems were not solved explicitly.
600 *
601  IF( abs( d( i ) ).LE.tol ) THEN
602  CALL claset( 'A', 1, nrhs, czero, czero, work( bx+i-1 ), n )
603  ELSE
604  rank = rank + 1
605  CALL clascl( 'G', 0, 0, d( i ), one, 1, nrhs,
606  $ work( bx+i-1 ), n, info )
607  END IF
608  d( i ) = abs( d( i ) )
609  250 continue
610 *
611 * Now apply back the right singular vectors.
612 *
613  icmpq2 = 1
614  DO 320 i = 1, nsub
615  st = iwork( i )
616  st1 = st - 1
617  nsize = iwork( sizei+i-1 )
618  bxst = bx + st1
619  IF( nsize.EQ.1 ) THEN
620  CALL ccopy( nrhs, work( bxst ), n, b( st, 1 ), ldb )
621  ELSE IF( nsize.LE.smlsiz ) THEN
622 *
623 * Since B and BX are complex, the following call to SGEMM
624 * is performed in two steps (real and imaginary parts).
625 *
626 * CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
627 * $ RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
628 * $ B( ST, 1 ), LDB )
629 *
630  j = bxst - n - 1
631  jreal = irwb - 1
632  DO 270 jcol = 1, nrhs
633  j = j + n
634  DO 260 jrow = 1, nsize
635  jreal = jreal + 1
636  rwork( jreal ) = REAL( WORK( J+JROW ) )
637  260 continue
638  270 continue
639  CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
640  $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
641  $ rwork( irwrb ), nsize )
642  j = bxst - n - 1
643  jimag = irwb - 1
644  DO 290 jcol = 1, nrhs
645  j = j + n
646  DO 280 jrow = 1, nsize
647  jimag = jimag + 1
648  rwork( jimag ) = aimag( work( j+jrow ) )
649  280 continue
650  290 continue
651  CALL sgemm( 'T', 'N', nsize, nrhs, nsize, one,
652  $ rwork( vt+st1 ), n, rwork( irwb ), nsize, zero,
653  $ rwork( irwib ), nsize )
654  jreal = irwrb - 1
655  jimag = irwib - 1
656  DO 310 jcol = 1, nrhs
657  DO 300 jrow = st, st + nsize - 1
658  jreal = jreal + 1
659  jimag = jimag + 1
660  b( jrow, jcol ) = cmplx( rwork( jreal ),
661  $ rwork( jimag ) )
662  300 continue
663  310 continue
664  ELSE
665  CALL clalsa( icmpq2, smlsiz, nsize, nrhs, work( bxst ), n,
666  $ b( st, 1 ), ldb, rwork( u+st1 ), n,
667  $ rwork( vt+st1 ), iwork( k+st1 ),
668  $ rwork( difl+st1 ), rwork( difr+st1 ),
669  $ rwork( z+st1 ), rwork( poles+st1 ),
670  $ iwork( givptr+st1 ), iwork( givcol+st1 ), n,
671  $ iwork( perm+st1 ), rwork( givnum+st1 ),
672  $ rwork( c+st1 ), rwork( s+st1 ),
673  $ rwork( nrwork ), iwork( iwk ), info )
674  IF( info.NE.0 ) THEN
675  return
676  END IF
677  END IF
678  320 continue
679 *
680 * Unscale and sort the singular values.
681 *
682  CALL slascl( 'G', 0, 0, one, orgnrm, n, 1, d, n, info )
683  CALL slasrt( 'D', n, d, info )
684  CALL clascl( 'G', 0, 0, orgnrm, one, n, nrhs, b, ldb, info )
685 *
686  return
687 *
688 * End of CLALSD
689 *
690  END