001:       SUBROUTINE CLALSD( UPLO, SMLSIZ, N, NRHS, D, E, B, LDB, RCOND,
002:      $                   RANK, WORK, RWORK, IWORK, INFO )
003: *
004: *  -- LAPACK routine (version 3.2) --
005: *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
006: *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
007: *     November 2006
008: *
009: *     .. Scalar Arguments ..
010:       CHARACTER          UPLO
011:       INTEGER            INFO, LDB, N, NRHS, RANK, SMLSIZ
012:       REAL               RCOND
013: *     ..
014: *     .. Array Arguments ..
015:       INTEGER            IWORK( * )
016:       REAL               D( * ), E( * ), RWORK( * )
017:       COMPLEX            B( LDB, * ), WORK( * )
018: *     ..
019: *
020: *  Purpose
021: *  =======
022: *
023: *  CLALSD uses the singular value decomposition of A to solve the least
024: *  squares problem of finding X to minimize the Euclidean norm of each
025: *  column of A*X-B, where A is N-by-N upper bidiagonal, and X and B
026: *  are N-by-NRHS. The solution X overwrites B.
027: *
028: *  The singular values of A smaller than RCOND times the largest
029: *  singular value are treated as zero in solving the least squares
030: *  problem; in this case a minimum norm solution is returned.
031: *  The actual singular values are returned in D in ascending order.
032: *
033: *  This code makes very mild assumptions about floating point
034: *  arithmetic. It will work on machines with a guard digit in
035: *  add/subtract, or on those binary machines without guard digits
036: *  which subtract like the Cray XMP, Cray YMP, Cray C 90, or Cray 2.
037: *  It could conceivably fail on hexadecimal or decimal machines
038: *  without guard digits, but we know of none.
039: *
040: *  Arguments
041: *  =========
042: *
043: *  UPLO   (input) CHARACTER*1
044: *         = 'U': D and E define an upper bidiagonal matrix.
045: *         = 'L': D and E define a  lower bidiagonal matrix.
046: *
047: *  SMLSIZ (input) INTEGER
048: *         The maximum size of the subproblems at the bottom of the
049: *         computation tree.
050: *
051: *  N      (input) INTEGER
052: *         The dimension of the  bidiagonal matrix.  N >= 0.
053: *
054: *  NRHS   (input) INTEGER
055: *         The number of columns of B. NRHS must be at least 1.
056: *
057: *  D      (input/output) REAL array, dimension (N)
058: *         On entry D contains the main diagonal of the bidiagonal
059: *         matrix. On exit, if INFO = 0, D contains its singular values.
060: *
061: *  E      (input/output) REAL array, dimension (N-1)
062: *         Contains the super-diagonal entries of the bidiagonal matrix.
063: *         On exit, E has been destroyed.
064: *
065: *  B      (input/output) COMPLEX array, dimension (LDB,NRHS)
066: *         On input, B contains the right hand sides of the least
067: *         squares problem. On output, B contains the solution X.
068: *
069: *  LDB    (input) INTEGER
070: *         The leading dimension of B in the calling subprogram.
071: *         LDB must be at least max(1,N).
072: *
073: *  RCOND  (input) REAL
074: *         The singular values of A less than or equal to RCOND times
075: *         the largest singular value are treated as zero in solving
076: *         the least squares problem. If RCOND is negative,
077: *         machine precision is used instead.
078: *         For example, if diag(S)*X=B were the least squares problem,
079: *         where diag(S) is a diagonal matrix of singular values, the
080: *         solution would be X(i) = B(i) / S(i) if S(i) is greater than
081: *         RCOND*max(S), and X(i) = 0 if S(i) is less than or equal to
082: *         RCOND*max(S).
083: *
084: *  RANK   (output) INTEGER
085: *         The number of singular values of A greater than RCOND times
086: *         the largest singular value.
087: *
088: *  WORK   (workspace) COMPLEX array, dimension (N * NRHS).
089: *
090: *  RWORK  (workspace) REAL array, dimension at least
091: *         (9*N + 2*N*SMLSIZ + 8*N*NLVL + 3*SMLSIZ*NRHS + (SMLSIZ+1)**2),
092: *         where
093: *         NLVL = MAX( 0, INT( LOG_2( MIN( M,N )/(SMLSIZ+1) ) ) + 1 )
094: *
095: *  IWORK  (workspace) INTEGER array, dimension (3*N*NLVL + 11*N).
096: *
097: *  INFO   (output) INTEGER
098: *         = 0:  successful exit.
099: *         < 0:  if INFO = -i, the i-th argument had an illegal value.
100: *         > 0:  The algorithm failed to compute an singular value while
101: *               working on the submatrix lying in rows and columns
102: *               INFO/(N+1) through MOD(INFO,N+1).
103: *
104: *  Further Details
105: *  ===============
106: *
107: *  Based on contributions by
108: *     Ming Gu and Ren-Cang Li, Computer Science Division, University of
109: *       California at Berkeley, USA
110: *     Osni Marques, LBNL/NERSC, USA
111: *
112: *  =====================================================================
113: *
114: *     .. Parameters ..
115:       REAL               ZERO, ONE, TWO
116:       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0, TWO = 2.0E0 )
117:       COMPLEX            CZERO
118:       PARAMETER          ( CZERO = ( 0.0E0, 0.0E0 ) )
119: *     ..
120: *     .. Local Scalars ..
121:       INTEGER            BX, BXST, C, DIFL, DIFR, GIVCOL, GIVNUM,
122:      $                   GIVPTR, I, ICMPQ1, ICMPQ2, IRWB, IRWIB, IRWRB,
123:      $                   IRWU, IRWVT, IRWWRK, IWK, J, JCOL, JIMAG,
124:      $                   JREAL, JROW, K, NLVL, NM1, NRWORK, NSIZE, NSUB,
125:      $                   PERM, POLES, S, SIZEI, SMLSZP, SQRE, ST, ST1,
126:      $                   U, VT, Z
127:       REAL               CS, EPS, ORGNRM, R, RCND, SN, TOL
128: *     ..
129: *     .. External Functions ..
130:       INTEGER            ISAMAX
131:       REAL               SLAMCH, SLANST
132:       EXTERNAL           ISAMAX, SLAMCH, SLANST
133: *     ..
134: *     .. External Subroutines ..
135:       EXTERNAL           CCOPY, CLACPY, CLALSA, CLASCL, CLASET, CSROT,
136:      $                   SGEMM, SLARTG, SLASCL, SLASDA, SLASDQ, SLASET,
137:      $                   SLASRT, XERBLA
138: *     ..
139: *     .. Intrinsic Functions ..
140:       INTRINSIC          ABS, AIMAG, CMPLX, INT, LOG, REAL, SIGN
141: *     ..
142: *     .. Executable Statements ..
143: *
144: *     Test the input parameters.
145: *
146:       INFO = 0
147: *
148:       IF( N.LT.0 ) THEN
149:          INFO = -3
150:       ELSE IF( NRHS.LT.1 ) THEN
151:          INFO = -4
152:       ELSE IF( ( LDB.LT.1 ) .OR. ( LDB.LT.N ) ) THEN
153:          INFO = -8
154:       END IF
155:       IF( INFO.NE.0 ) THEN
156:          CALL XERBLA( 'CLALSD', -INFO )
157:          RETURN
158:       END IF
159: *
160:       EPS = SLAMCH( 'Epsilon' )
161: *
162: *     Set up the tolerance.
163: *
164:       IF( ( RCOND.LE.ZERO ) .OR. ( RCOND.GE.ONE ) ) THEN
165:          RCND = EPS
166:       ELSE
167:          RCND = RCOND
168:       END IF
169: *
170:       RANK = 0
171: *
172: *     Quick return if possible.
173: *
174:       IF( N.EQ.0 ) THEN
175:          RETURN
176:       ELSE IF( N.EQ.1 ) THEN
177:          IF( D( 1 ).EQ.ZERO ) THEN
178:             CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B, LDB )
179:          ELSE
180:             RANK = 1
181:             CALL CLASCL( 'G', 0, 0, D( 1 ), ONE, 1, NRHS, B, LDB, INFO )
182:             D( 1 ) = ABS( D( 1 ) )
183:          END IF
184:          RETURN
185:       END IF
186: *
187: *     Rotate the matrix if it is lower bidiagonal.
188: *
189:       IF( UPLO.EQ.'L' ) THEN
190:          DO 10 I = 1, N - 1
191:             CALL SLARTG( D( I ), E( I ), CS, SN, R )
192:             D( I ) = R
193:             E( I ) = SN*D( I+1 )
194:             D( I+1 ) = CS*D( I+1 )
195:             IF( NRHS.EQ.1 ) THEN
196:                CALL CSROT( 1, B( I, 1 ), 1, B( I+1, 1 ), 1, CS, SN )
197:             ELSE
198:                RWORK( I*2-1 ) = CS
199:                RWORK( I*2 ) = SN
200:             END IF
201:    10    CONTINUE
202:          IF( NRHS.GT.1 ) THEN
203:             DO 30 I = 1, NRHS
204:                DO 20 J = 1, N - 1
205:                   CS = RWORK( J*2-1 )
206:                   SN = RWORK( J*2 )
207:                   CALL CSROT( 1, B( J, I ), 1, B( J+1, I ), 1, CS, SN )
208:    20          CONTINUE
209:    30       CONTINUE
210:          END IF
211:       END IF
212: *
213: *     Scale.
214: *
215:       NM1 = N - 1
216:       ORGNRM = SLANST( 'M', N, D, E )
217:       IF( ORGNRM.EQ.ZERO ) THEN
218:          CALL CLASET( 'A', N, NRHS, CZERO, CZERO, B, LDB )
219:          RETURN
220:       END IF
221: *
222:       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, N, 1, D, N, INFO )
223:       CALL SLASCL( 'G', 0, 0, ORGNRM, ONE, NM1, 1, E, NM1, INFO )
224: *
225: *     If N is smaller than the minimum divide size SMLSIZ, then solve
226: *     the problem with another solver.
227: *
228:       IF( N.LE.SMLSIZ ) THEN
229:          IRWU = 1
230:          IRWVT = IRWU + N*N
231:          IRWWRK = IRWVT + N*N
232:          IRWRB = IRWWRK
233:          IRWIB = IRWRB + N*NRHS
234:          IRWB = IRWIB + N*NRHS
235:          CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWU ), N )
236:          CALL SLASET( 'A', N, N, ZERO, ONE, RWORK( IRWVT ), N )
237:          CALL SLASDQ( 'U', 0, N, N, N, 0, D, E, RWORK( IRWVT ), N,
238:      $                RWORK( IRWU ), N, RWORK( IRWWRK ), 1,
239:      $                RWORK( IRWWRK ), INFO )
240:          IF( INFO.NE.0 ) THEN
241:             RETURN
242:          END IF
243: *
244: *        In the real version, B is passed to SLASDQ and multiplied
245: *        internally by Q'. Here B is complex and that product is
246: *        computed below in two steps (real and imaginary parts).
247: *
248:          J = IRWB - 1
249:          DO 50 JCOL = 1, NRHS
250:             DO 40 JROW = 1, N
251:                J = J + 1
252:                RWORK( J ) = REAL( B( JROW, JCOL ) )
253:    40       CONTINUE
254:    50    CONTINUE
255:          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
256:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
257:          J = IRWB - 1
258:          DO 70 JCOL = 1, NRHS
259:             DO 60 JROW = 1, N
260:                J = J + 1
261:                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
262:    60       CONTINUE
263:    70    CONTINUE
264:          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWU ), N,
265:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
266:          JREAL = IRWRB - 1
267:          JIMAG = IRWIB - 1
268:          DO 90 JCOL = 1, NRHS
269:             DO 80 JROW = 1, N
270:                JREAL = JREAL + 1
271:                JIMAG = JIMAG + 1
272:                B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
273:    80       CONTINUE
274:    90    CONTINUE
275: *
276:          TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
277:          DO 100 I = 1, N
278:             IF( D( I ).LE.TOL ) THEN
279:                CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, B( I, 1 ), LDB )
280:             ELSE
281:                CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS, B( I, 1 ),
282:      $                      LDB, INFO )
283:                RANK = RANK + 1
284:             END IF
285:   100    CONTINUE
286: *
287: *        Since B is complex, the following call to SGEMM is performed
288: *        in two steps (real and imaginary parts). That is for V * B
289: *        (in the real version of the code V' is stored in WORK).
290: *
291: *        CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, WORK, N, B, LDB, ZERO,
292: *    $               WORK( NWORK ), N )
293: *
294:          J = IRWB - 1
295:          DO 120 JCOL = 1, NRHS
296:             DO 110 JROW = 1, N
297:                J = J + 1
298:                RWORK( J ) = REAL( B( JROW, JCOL ) )
299:   110       CONTINUE
300:   120    CONTINUE
301:          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
302:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWRB ), N )
303:          J = IRWB - 1
304:          DO 140 JCOL = 1, NRHS
305:             DO 130 JROW = 1, N
306:                J = J + 1
307:                RWORK( J ) = AIMAG( B( JROW, JCOL ) )
308:   130       CONTINUE
309:   140    CONTINUE
310:          CALL SGEMM( 'T', 'N', N, NRHS, N, ONE, RWORK( IRWVT ), N,
311:      $               RWORK( IRWB ), N, ZERO, RWORK( IRWIB ), N )
312:          JREAL = IRWRB - 1
313:          JIMAG = IRWIB - 1
314:          DO 160 JCOL = 1, NRHS
315:             DO 150 JROW = 1, N
316:                JREAL = JREAL + 1
317:                JIMAG = JIMAG + 1
318:                B( JROW, JCOL ) = CMPLX( RWORK( JREAL ), RWORK( JIMAG ) )
319:   150       CONTINUE
320:   160    CONTINUE
321: *
322: *        Unscale.
323: *
324:          CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
325:          CALL SLASRT( 'D', N, D, INFO )
326:          CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
327: *
328:          RETURN
329:       END IF
330: *
331: *     Book-keeping and setting up some constants.
332: *
333:       NLVL = INT( LOG( REAL( N ) / REAL( SMLSIZ+1 ) ) / LOG( TWO ) ) + 1
334: *
335:       SMLSZP = SMLSIZ + 1
336: *
337:       U = 1
338:       VT = 1 + SMLSIZ*N
339:       DIFL = VT + SMLSZP*N
340:       DIFR = DIFL + NLVL*N
341:       Z = DIFR + NLVL*N*2
342:       C = Z + NLVL*N
343:       S = C + N
344:       POLES = S + N
345:       GIVNUM = POLES + 2*NLVL*N
346:       NRWORK = GIVNUM + 2*NLVL*N
347:       BX = 1
348: *
349:       IRWRB = NRWORK
350:       IRWIB = IRWRB + SMLSIZ*NRHS
351:       IRWB = IRWIB + SMLSIZ*NRHS
352: *
353:       SIZEI = 1 + N
354:       K = SIZEI + N
355:       GIVPTR = K + N
356:       PERM = GIVPTR + N
357:       GIVCOL = PERM + NLVL*N
358:       IWK = GIVCOL + NLVL*N*2
359: *
360:       ST = 1
361:       SQRE = 0
362:       ICMPQ1 = 1
363:       ICMPQ2 = 0
364:       NSUB = 0
365: *
366:       DO 170 I = 1, N
367:          IF( ABS( D( I ) ).LT.EPS ) THEN
368:             D( I ) = SIGN( EPS, D( I ) )
369:          END IF
370:   170 CONTINUE
371: *
372:       DO 240 I = 1, NM1
373:          IF( ( ABS( E( I ) ).LT.EPS ) .OR. ( I.EQ.NM1 ) ) THEN
374:             NSUB = NSUB + 1
375:             IWORK( NSUB ) = ST
376: *
377: *           Subproblem found. First determine its size and then
378: *           apply divide and conquer on it.
379: *
380:             IF( I.LT.NM1 ) THEN
381: *
382: *              A subproblem with E(I) small for I < NM1.
383: *
384:                NSIZE = I - ST + 1
385:                IWORK( SIZEI+NSUB-1 ) = NSIZE
386:             ELSE IF( ABS( E( I ) ).GE.EPS ) THEN
387: *
388: *              A subproblem with E(NM1) not too small but I = NM1.
389: *
390:                NSIZE = N - ST + 1
391:                IWORK( SIZEI+NSUB-1 ) = NSIZE
392:             ELSE
393: *
394: *              A subproblem with E(NM1) small. This implies an
395: *              1-by-1 subproblem at D(N), which is not solved
396: *              explicitly.
397: *
398:                NSIZE = I - ST + 1
399:                IWORK( SIZEI+NSUB-1 ) = NSIZE
400:                NSUB = NSUB + 1
401:                IWORK( NSUB ) = N
402:                IWORK( SIZEI+NSUB-1 ) = 1
403:                CALL CCOPY( NRHS, B( N, 1 ), LDB, WORK( BX+NM1 ), N )
404:             END IF
405:             ST1 = ST - 1
406:             IF( NSIZE.EQ.1 ) THEN
407: *
408: *              This is a 1-by-1 subproblem and is not solved
409: *              explicitly.
410: *
411:                CALL CCOPY( NRHS, B( ST, 1 ), LDB, WORK( BX+ST1 ), N )
412:             ELSE IF( NSIZE.LE.SMLSIZ ) THEN
413: *
414: *              This is a small subproblem and is solved by SLASDQ.
415: *
416:                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
417:      $                      RWORK( VT+ST1 ), N )
418:                CALL SLASET( 'A', NSIZE, NSIZE, ZERO, ONE,
419:      $                      RWORK( U+ST1 ), N )
420:                CALL SLASDQ( 'U', 0, NSIZE, NSIZE, NSIZE, 0, D( ST ),
421:      $                      E( ST ), RWORK( VT+ST1 ), N, RWORK( U+ST1 ),
422:      $                      N, RWORK( NRWORK ), 1, RWORK( NRWORK ),
423:      $                      INFO )
424:                IF( INFO.NE.0 ) THEN
425:                   RETURN
426:                END IF
427: *
428: *              In the real version, B is passed to SLASDQ and multiplied
429: *              internally by Q'. Here B is complex and that product is
430: *              computed below in two steps (real and imaginary parts).
431: *
432:                J = IRWB - 1
433:                DO 190 JCOL = 1, NRHS
434:                   DO 180 JROW = ST, ST + NSIZE - 1
435:                      J = J + 1
436:                      RWORK( J ) = REAL( B( JROW, JCOL ) )
437:   180             CONTINUE
438:   190          CONTINUE
439:                CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
440:      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
441:      $                     ZERO, RWORK( IRWRB ), NSIZE )
442:                J = IRWB - 1
443:                DO 210 JCOL = 1, NRHS
444:                   DO 200 JROW = ST, ST + NSIZE - 1
445:                      J = J + 1
446:                      RWORK( J ) = AIMAG( B( JROW, JCOL ) )
447:   200             CONTINUE
448:   210          CONTINUE
449:                CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
450:      $                     RWORK( U+ST1 ), N, RWORK( IRWB ), NSIZE,
451:      $                     ZERO, RWORK( IRWIB ), NSIZE )
452:                JREAL = IRWRB - 1
453:                JIMAG = IRWIB - 1
454:                DO 230 JCOL = 1, NRHS
455:                   DO 220 JROW = ST, ST + NSIZE - 1
456:                      JREAL = JREAL + 1
457:                      JIMAG = JIMAG + 1
458:                      B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
459:      $                                 RWORK( JIMAG ) )
460:   220             CONTINUE
461:   230          CONTINUE
462: *
463:                CALL CLACPY( 'A', NSIZE, NRHS, B( ST, 1 ), LDB,
464:      $                      WORK( BX+ST1 ), N )
465:             ELSE
466: *
467: *              A large problem. Solve it using divide and conquer.
468: *
469:                CALL SLASDA( ICMPQ1, SMLSIZ, NSIZE, SQRE, D( ST ),
470:      $                      E( ST ), RWORK( U+ST1 ), N, RWORK( VT+ST1 ),
471:      $                      IWORK( K+ST1 ), RWORK( DIFL+ST1 ),
472:      $                      RWORK( DIFR+ST1 ), RWORK( Z+ST1 ),
473:      $                      RWORK( POLES+ST1 ), IWORK( GIVPTR+ST1 ),
474:      $                      IWORK( GIVCOL+ST1 ), N, IWORK( PERM+ST1 ),
475:      $                      RWORK( GIVNUM+ST1 ), RWORK( C+ST1 ),
476:      $                      RWORK( S+ST1 ), RWORK( NRWORK ),
477:      $                      IWORK( IWK ), INFO )
478:                IF( INFO.NE.0 ) THEN
479:                   RETURN
480:                END IF
481:                BXST = BX + ST1
482:                CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, B( ST, 1 ),
483:      $                      LDB, WORK( BXST ), N, RWORK( U+ST1 ), N,
484:      $                      RWORK( VT+ST1 ), IWORK( K+ST1 ),
485:      $                      RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
486:      $                      RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
487:      $                      IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
488:      $                      IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
489:      $                      RWORK( C+ST1 ), RWORK( S+ST1 ),
490:      $                      RWORK( NRWORK ), IWORK( IWK ), INFO )
491:                IF( INFO.NE.0 ) THEN
492:                   RETURN
493:                END IF
494:             END IF
495:             ST = I + 1
496:          END IF
497:   240 CONTINUE
498: *
499: *     Apply the singular values and treat the tiny ones as zero.
500: *
501:       TOL = RCND*ABS( D( ISAMAX( N, D, 1 ) ) )
502: *
503:       DO 250 I = 1, N
504: *
505: *        Some of the elements in D can be negative because 1-by-1
506: *        subproblems were not solved explicitly.
507: *
508:          IF( ABS( D( I ) ).LE.TOL ) THEN
509:             CALL CLASET( 'A', 1, NRHS, CZERO, CZERO, WORK( BX+I-1 ), N )
510:          ELSE
511:             RANK = RANK + 1
512:             CALL CLASCL( 'G', 0, 0, D( I ), ONE, 1, NRHS,
513:      $                   WORK( BX+I-1 ), N, INFO )
514:          END IF
515:          D( I ) = ABS( D( I ) )
516:   250 CONTINUE
517: *
518: *     Now apply back the right singular vectors.
519: *
520:       ICMPQ2 = 1
521:       DO 320 I = 1, NSUB
522:          ST = IWORK( I )
523:          ST1 = ST - 1
524:          NSIZE = IWORK( SIZEI+I-1 )
525:          BXST = BX + ST1
526:          IF( NSIZE.EQ.1 ) THEN
527:             CALL CCOPY( NRHS, WORK( BXST ), N, B( ST, 1 ), LDB )
528:          ELSE IF( NSIZE.LE.SMLSIZ ) THEN
529: *
530: *           Since B and BX are complex, the following call to SGEMM
531: *           is performed in two steps (real and imaginary parts).
532: *
533: *           CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
534: *    $                  RWORK( VT+ST1 ), N, RWORK( BXST ), N, ZERO,
535: *    $                  B( ST, 1 ), LDB )
536: *
537:             J = BXST - N - 1
538:             JREAL = IRWB - 1
539:             DO 270 JCOL = 1, NRHS
540:                J = J + N
541:                DO 260 JROW = 1, NSIZE
542:                   JREAL = JREAL + 1
543:                   RWORK( JREAL ) = REAL( WORK( J+JROW ) )
544:   260          CONTINUE
545:   270       CONTINUE
546:             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
547:      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
548:      $                  RWORK( IRWRB ), NSIZE )
549:             J = BXST - N - 1
550:             JIMAG = IRWB - 1
551:             DO 290 JCOL = 1, NRHS
552:                J = J + N
553:                DO 280 JROW = 1, NSIZE
554:                   JIMAG = JIMAG + 1
555:                   RWORK( JIMAG ) = AIMAG( WORK( J+JROW ) )
556:   280          CONTINUE
557:   290       CONTINUE
558:             CALL SGEMM( 'T', 'N', NSIZE, NRHS, NSIZE, ONE,
559:      $                  RWORK( VT+ST1 ), N, RWORK( IRWB ), NSIZE, ZERO,
560:      $                  RWORK( IRWIB ), NSIZE )
561:             JREAL = IRWRB - 1
562:             JIMAG = IRWIB - 1
563:             DO 310 JCOL = 1, NRHS
564:                DO 300 JROW = ST, ST + NSIZE - 1
565:                   JREAL = JREAL + 1
566:                   JIMAG = JIMAG + 1
567:                   B( JROW, JCOL ) = CMPLX( RWORK( JREAL ),
568:      $                              RWORK( JIMAG ) )
569:   300          CONTINUE
570:   310       CONTINUE
571:          ELSE
572:             CALL CLALSA( ICMPQ2, SMLSIZ, NSIZE, NRHS, WORK( BXST ), N,
573:      $                   B( ST, 1 ), LDB, RWORK( U+ST1 ), N,
574:      $                   RWORK( VT+ST1 ), IWORK( K+ST1 ),
575:      $                   RWORK( DIFL+ST1 ), RWORK( DIFR+ST1 ),
576:      $                   RWORK( Z+ST1 ), RWORK( POLES+ST1 ),
577:      $                   IWORK( GIVPTR+ST1 ), IWORK( GIVCOL+ST1 ), N,
578:      $                   IWORK( PERM+ST1 ), RWORK( GIVNUM+ST1 ),
579:      $                   RWORK( C+ST1 ), RWORK( S+ST1 ),
580:      $                   RWORK( NRWORK ), IWORK( IWK ), INFO )
581:             IF( INFO.NE.0 ) THEN
582:                RETURN
583:             END IF
584:          END IF
585:   320 CONTINUE
586: *
587: *     Unscale and sort the singular values.
588: *
589:       CALL SLASCL( 'G', 0, 0, ONE, ORGNRM, N, 1, D, N, INFO )
590:       CALL SLASRT( 'D', N, D, INFO )
591:       CALL CLASCL( 'G', 0, 0, ORGNRM, ONE, N, NRHS, B, LDB, INFO )
592: *
593:       RETURN
594: *
595: *     End of CLALSD
596: *
597:       END
598: