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